This software provides functions for the fast, hardware-accelerated simulation of therapeutic ultrasound beams. The package consists of C/C++ libraries with a MATLAB interface and interface functions to FOCUS.

The software is developed at the Ferrarra Laboratory in the Department of Biomedical Engineering at the University of California at Davis and is supported by NIH3R21EB009434 and by ARRA supplement 3R21EB009434-01S1.

Current Release:

- Mac OSX 10.6, Intel 32 Bit (dmg/mpkg) [abt. 30MB] [md5: 48234cce15d2d4f2436810a654a94c26 ]
- Windows, x84 32 Bit (ZIP) [abt. 30MB] [md5: 8569ebd424179f0181f919398949f7b3 ]

Version | Date | Comment |
---|---|---|

0.0.2 | 2011-01-09 | Added Windows and Linux versions. |

0.0.1 | 2010-11-09 | Initial Release of OSX version. |

- You will need an Intel Macintosh computer that has a reasonably current NVIDIA graphics processor installed. All modern MacBook Pro's should be sufficient.

- A current version of Matlab is sufficient for the non-GPU version of the software. To run GPU-based code, you need NVIDIA's libraries and a compatible graphics board installed.

- The binary requires a 64bit linux system and has been compiled on an Ubuntu Linux system.

We provide two user interfaces to our code: A native interface (indicated by the ML prefix) and an interface that tries to replicate the function calls of FOCUS. If you use the FOCUS interface, ensure that you have FOCUS installed and running on your machine. In case of OSX, where no FOCUS binaries are provided, our functions can be used to replace some of FOCUS' functionality. The examples below only use FOCUS' MATLAB files and do not rely on the binaries.

The code shown in blue boxes shows our native interface, the code shown in green boxes uses FOCUS' interface.

MLRectangularPistonGrid

??? Error using ==> MLRectangularPistonGrid

Must have five or six input arguments: ([4x4 transform], [sx, sy, sz], l, s, k, ['cuda'|'cudadp'|'cpu'|'cpusp'|'cpudp'], [nabscissas])

Compute the near-field in front of a piston:
??? Error using ==> MLRectangularPistonGrid

Must have five or six input arguments: ([4x4 transform], [sx, sy, sz], l, s, k, ['cuda'|'cudadp'|'cpu'|'cpusp'|'cpudp'], [nabscissas])

%% setup the basic parameters

f0=1540000; % excitation frequency (Hz)

c=1500; % speed of sound (m/s)

lambda = c/f0; % wavelength (m)

k=2*pi/lambda; % wavenumber (1/m)

rho=1000; % acoustic density

%% setup piston parameters and the computational grid

s=5.0*lambda; % piston width

l=7.5*lambda; % piston height

siz = [ 61 1 101 ]; % grid sizes

sampling = [s/2/(siz(1)-1)/2*3 1 s*s/4/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

%% do the calculation

pressure = MLRectangularPistonGrid( transform, siz, l, s, k, 'cuda', 40)*(1j*c*rho);

%% display and store the result

figure(1); surf(abs(squeeze(pressure))); axis tight

saveas(figure(1), 'figure1.png')

f0=1540000; % excitation frequency (Hz)

c=1500; % speed of sound (m/s)

lambda = c/f0; % wavelength (m)

k=2*pi/lambda; % wavenumber (1/m)

rho=1000; % acoustic density

%% setup piston parameters and the computational grid

s=5.0*lambda; % piston width

l=7.5*lambda; % piston height

siz = [ 61 1 101 ]; % grid sizes

sampling = [s/2/(siz(1)-1)/2*3 1 s*s/4/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

%% do the calculation

pressure = MLRectangularPistonGrid( transform, siz, l, s, k, 'cuda', 40)*(1j*c*rho);

%% display and store the result

figure(1); surf(abs(squeeze(pressure))); axis tight

saveas(figure(1), 'figure1.png')

%% setup the basic parameters

f0=1540000;

define_media;

lambda = lossless.soundspeed/f0;

%% setup piston parameters and the computational grid

s=5.0*lambda;

l=7.5*lambda;

siz = [61 1 101];

sampling = [s/2/(siz(1)-1)/2*3 1 s*s/4/lambda/(siz(3)-1)];

ps=set_coordinate_grid(sampling,0,1.5*s/2,0,0,0,s^2/4/lambda);

transducer=get_rect(s/2, l/2, [0 0 0], [0 0 0]);

%% do the calculation

pressure=gpu_fnm_call(transducer,ps,lossless,40,f0,0);

%% display and store the result

figure(1); surf(abs(squeeze(pressure))); axis tight

saveas(figure(1), 'figure1.png')

This should result in a picture file called "figure1.png" that looks as follows:f0=1540000;

define_media;

lambda = lossless.soundspeed/f0;

%% setup piston parameters and the computational grid

s=5.0*lambda;

l=7.5*lambda;

siz = [61 1 101];

sampling = [s/2/(siz(1)-1)/2*3 1 s*s/4/lambda/(siz(3)-1)];

ps=set_coordinate_grid(sampling,0,1.5*s/2,0,0,0,s^2/4/lambda);

transducer=get_rect(s/2, l/2, [0 0 0], [0 0 0]);

%% do the calculation

pressure=gpu_fnm_call(transducer,ps,lossless,40,f0,0);

%% display and store the result

figure(1); surf(abs(squeeze(pressure))); axis tight

saveas(figure(1), 'figure1.png')

MLCircularPistonGrid

??? Error using ==> MLCircularPistonGrid

Must have four or five input argument: (transform4x4, [sizex, sizey, sizez], radius, k, [, 'cuda' | 'gpu',[ abscissas]]

with the following parameters
??? Error using ==> MLCircularPistonGrid

Must have four or five input argument: (transform4x4, [sizex, sizey, sizez], radius, k, [, 'cuda' | 'gpu',[ abscissas]]

f0=1000000; % excitation frequency

c=1500; % speed of sound

lambda = c/f0; % wavelength

k=2*pi/lambda; % wavenumber

radius=5.0*lambda; % piston radius

siz = [ 61 1 101 ];

sampling = [1.5*radius/(siz(1)-1) 1 radius^2/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

rho=1000;

pressure = MLCircularPistonGrid( transform, siz, radius, k, 'cpu', 40)*(1j*rho*c);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

c=1500; % speed of sound

lambda = c/f0; % wavelength

k=2*pi/lambda; % wavenumber

radius=5.0*lambda; % piston radius

siz = [ 61 1 101 ];

sampling = [1.5*radius/(siz(1)-1) 1 radius^2/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

rho=1000;

pressure = MLCircularPistonGrid( transform, siz, radius, k, 'cpu', 40)*(1j*rho*c);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

f0=1540000;

define_media;

lambda = lossless.soundspeed/f0;

radius=5.0*lambda;

siz = [61 1 101];

sampling = [1.5*radius/(siz(1)-1) 1 radius^2/lambda/(siz(3)-1)];

ps=set_coordinate_grid(sampling,0,sampling(1)*(siz(1)-1),0,0,0,sampling(3)*(siz(3)-1));

transducer=get_circ(radius, [0 0 0], [0 0 0]);

pressure=cpu_fnm_call(transducer,ps,lossless,40,f0,0);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

which should result in a picture file called "figure2.png" that looks as follows:define_media;

lambda = lossless.soundspeed/f0;

radius=5.0*lambda;

siz = [61 1 101];

sampling = [1.5*radius/(siz(1)-1) 1 radius^2/lambda/(siz(3)-1)];

ps=set_coordinate_grid(sampling,0,sampling(1)*(siz(1)-1),0,0,0,sampling(3)*(siz(3)-1));

transducer=get_circ(radius, [0 0 0], [0 0 0]);

pressure=cpu_fnm_call(transducer,ps,lossless,40,f0,0);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

%% setup parameters

f0=1000000; % excitation frequency

c=1500; % speed of sound

lambda = c/f0; % wavelength

k=2*pi/lambda; % wavenumber

radius=5.0*lambda; % piston radius

siz = [ 256 1 101 ];

sampling = [1.5*radius/(61-1) 1 radius*radius/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

%% compute a single piston

pressure = MLCircularPistonGrid( transform, siz, radius, k, 'cpu', 40);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

%% compute the convolution of four elements, triggered synchronously

array_pressure = MLRectilinearArrayConvolution( pressure, [1; 1; 1; 1], 4, 1, 64, 1, 'cuda')

%% draw figure

figure(3); surf(abs(squeeze(array_pressure))); axis tight; view(225,45);

saveas(figure(3), 'figure3.png');

The result is a 64x1x101 grid that is aligned with the outer of the four pistons as shown here:f0=1000000; % excitation frequency

c=1500; % speed of sound

lambda = c/f0; % wavelength

k=2*pi/lambda; % wavenumber

radius=5.0*lambda; % piston radius

siz = [ 256 1 101 ];

sampling = [1.5*radius/(61-1) 1 radius*radius/lambda/(siz(3)-1)];

transform = [ sampling(1) 0 0 0; 0 sampling(2) 0 0; 0 0 sampling(3) 0; 0 0 0 1 ];

%% compute a single piston

pressure = MLCircularPistonGrid( transform, siz, radius, k, 'cpu', 40);

figure(2); surf(abs(squeeze(pressure))); axis tight; view(225,45);

saveas(figure(2), 'figure2.png')

%% compute the convolution of four elements, triggered synchronously

array_pressure = MLRectilinearArrayConvolution( pressure, [1; 1; 1; 1], 4, 1, 64, 1, 'cuda')

%% draw figure

figure(3); surf(abs(squeeze(array_pressure))); axis tight; view(225,45);

saveas(figure(3), 'figure3.png');

Keep in mind that any call to MATLAB involves copying data from MATLAB to C/C++ to the GPU. This may slow down some operations tremendously.

Fast ultrasound beam prediction for linear and regular two-dimensional arrays using a general-purpose graphics processing unit
Abstract Book of the 2010 IEEE International Ultrasonics Symposium (IUS), San Diego, CA , IEEE, UFFC pp. 619, P5-M4-9 October 11-14 2010 [conference] | |

Abstract: In our preliminary work, we have demonstrated that ultrasonic drug release using mild hyperthermia is feasible; temperature sensitive drug delivery particles are injected and allowed to accumulate within a region of interest, typically over 6-24 hours. Ultrasound is then scanned through the volume of interest to locally increase the temperature by 2-4 degree celsius for a short period, releasing a drug and increasing its efficacy. In order to translate these methods, ultrasound AND drug dose must each be predicted and validated. We have made great strides on drug dose—we now tackle ultrasound dose for this special problem of mild hyperthermia. Predicting the 3-D heating profile in tissue as a function of time requires knowledge of the acoustic beam intensity as a function of time as it is steered in 3-D space. Ultrasound beams generated from an aperture large enough to deliver the acoustic power necessary for tissue heating must be focused to some degree in order to work at reasonable depths from the aperture face and at depth-of-fields that are not much larger that the desired heating region. Focusing necessitates that the beam be scanned in order to treat larger volumes, such as the volume encompassing a tumor. The beam scanning is accomplished dynamically, in real time, and knowledge of the beam is necessary in order to make predictions on where energy is being deposited and how much is being deposited. This information is important for predicting thermal dose according to Arrhenius-based methods such as Cumulative Equivalent Minutes at 43°C CEM4). In summary, we present our initial results using GPUs to accelerate the prediction of large-scale acoustic fields relevant to ultrasound hyperthermia. The FNM was chosen over other methods due to its rapid convergence and accuracy in the near- field up to the surface of the radiator. Algorithms for regularly-sampled 1-D and 2-D RSM Times (in milliseconds) FNM Times (in milliseconds) piston-arrays were presented and benchmarked to demonstrate the power of GPU- MinAccuracy 1.0% 0.1% 0.01% 1.0% 0.1% 0.01% Abscissas based field prediction for scenarios commonly found in ultrasound therapy. Our goal of generating the complex acoustic field for a volume containing 10 million points within 1 second is easily achieved using the algorithms and hardware presented. Our future goal is to incorporate this new ability into our ultrasound hyperthermia system for purposes of treatment planning and thermal dose prediction and control. | |

Function approximations to accelerate 3-D beam predictions for thermal dose calculations.
Proceedings of the 2009 IEEE International Ultrasonics Symposium pp. 1777-1779 September 20-23 2009 [ISBN: 978-1-4244-4389-5] [conference] | |

Abstract: We demonstrate that ultrasound beam predictions are accelerated by 22x using a fast approximation to evaluate complex exponentials with insignificant errors compared to full precision calculations. Since TAP is the more relevant quantity for thermal diffusion-limited heating, the small error of less than 0.1% indicates a better measure of the fast approximation compared to maximum intensity error of less than 1%. It appears in these preliminary results that errors in the beam prediction can be reduced by selecting a value for L closer to that which minimizes error for sine/cosine pairs based on the Pythagorean trigonometric identity. The same methodology may be applied to other beam prediction algorithms to speed their execution at a relatively minor degradation in accuracy. |

- 2011-02-04: Fixed missing semicolon in example code and corrected angle of view
- 2010-10-08: initial release