k-Wave

1. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi, Dr. Cox and Dr. Treeby,

   Now I am trying to run a 3D simulation, and Cartesian 'sensor.mask' is used. The problem is that if 'CartInterp' is set to be 'linear' and 'DataCast' is set to be 'single', then an error will prompt, saying:

   ??? Error using ==> subsasgn
   The input points must be a double array.

   Error in ==> kspaceFirstOrder3D at 935
   F_interp.V = reshape(p, [], 1);

   However, if I used
   'input_args = {'CartInterp', 'linear', 'DataCast','off'}'
   or
   'input_args = {'CartInterp', 'nearest', 'DataCast','single'}',
   then there is no problem in both of these 2 cases.

   I think the culprit may be the 'F_interp.V', which perhaps needs to be cast to 'single', too.

   Chao

   Posted 13 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Chao,

   Unfortunately this is a limitation of the MATLAB function TriScatteredInterp which is used for linear interpolation within k-Wave. This function only supports doubles. For example, try

   >> F = TriScatteredInterp(single(rand(100,1)), single(rand(100,1)),
   single(rand(100,1)));
   ??? Error using ==> TriScatteredInterp
   The input points must be a double array.

   In the 2D case, the k-Wave code checks if the parameter 'DataCast' is used, and if so it switches to using gridDataFast instead of TriScatteredInterp. However, gridDataFast only works in 2D. At this stage, the two options you mention are the best way to avoid this error.

   There are some notes in the documentation here.

   Kind regards,

   Brad.

   Posted 13 years ago #

3. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi,

   I just found that there may be a mistake in the help document of the optimising k-Wave performance example.

   It says, quote, "Note, the interpolation function used within kspaceFirstOrder2D and kspaceFirstOrder3D does not currently support GPU usage, so the optional input parameter 'CartInterp' should be set to 'nearest' if using a Cartesian sensor mask." But actually, when calling 'kspaceFirstOrder2D', 'CartInterp' can be set to 'linear' when 'DataCast' is set to 'gsingle', because the 'gridDataFast' used for interpolation in 'kspaceFirstOrder2D' supports both 'single' and 'gsingle'.

   Best,
   Chao

   Posted 12 years ago #

4. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi, Dr. Treeby,

   I found that the same interpolation results given by 'gridDataFast' are slightly different from the ones from 'TriScatteredInterp'. And, as you mentioned in the previous post, 'TriScatteredInterp' cannot be used with 'single' or 'gsingle' data type.

   So I modified the 'gridDataFast' so that it will give the exact same interpolation results as 'TriScatteredInterp', and it supports 'DataCast', which means we can use GPUs to accelerate the computation in 3D case if 'CartInterp' is set to be 'linear'.

   In 2D case, the modified version is:

   function [tri, bc] = gridDataFast2D(x, y, xi, yi)

   % enforce x, y, xi, yi to be column vectors
   x = x(:);
   y = y(:);
   xi = xi(:);
   yi = yi(:);

   % triangulize the data
   tri = DelaunayTri(x, y);

   % catch trinagulation error
   if isempty(tri)
   error('Data cannot be triangulated.');
   end

   % find the nearest triangle and the corresponding Barycentric coordinates
   [t, bc] = pointLocation(tri,[xi yi]);

   tri = tri(t,:);

   After the 'tri' and 'bc' are computed in 'kspaceFirstOrder_createStorageVariables', the 'sensor_data' can be calculated by sensor_data(:, t_index) = sum(p(sensor.tri) .* sensor.bc, 2) in the kspaceFirstOrder2D.

   Likely, the 3D version is:

   function [tri, bc] = gridDataFast3D(x, y, z, xi, yi, zi)

   % enforce x, y, z, xi, yi, zi to be column vectors
   x = x(:);
   y = y(:);
   z = z(:);
   xi = xi(:);
   yi = yi(:);
   zi = zi(:);

   % triangulize the data
   tri = DelaunayTri(x, y, z);

   % catch trinagulation error
   if isempty(tri)
   error('Data cannot be triangulated.');
   end

   % find the nearest triangle and the corresponding Barycentric coordinates
   [t, bc] = pointLocation(tri,[xi yi zi]);

   tri = tri(t,:);

   Similarly, After the 'tri' and 'bc' are computed in 'kspaceFirstOrder_createStorageVariables', the 'sensor_data' can be calculated by sensor_data(:, t_index) = sum(p(sensor.tri) .* sensor.bc, 2) in the kspaceFirstOrder3D.

   Also, I found the modified 'gridDataFast2D' and 'gridDataFast3D' are faster than 'gridDataFast' and 'TriScatteredInterp'. For example, in 2D case, we had:

   tic;
   [tri, bc] = gridDataFast2D(kgrid.x, kgrid.y, sensor.mask(1,:), sensor.mask(2,:));
   i1 = sum(source.p0(tri) .* bc,2);
   toc
   Elapsed time is 1.072967 seconds.

   tic;
   [zi, del_tri] = gridDataFast(kgrid.x, kgrid.y, zeros(kgrid.Nx, kgrid.Ny), sensor.mask(1,:), sensor.mask(2,:));
   i2 = gridDataFast(kgrid.x, kgrid.y, source.p0, sensor.mask(1,:), sensor.mask(2,:), del_tri);
   toc
   Elapsed time is 5.227701 seconds.

   tic;
   F_interp = TriScatteredInterp(reshape(kgrid.x, [], 1), reshape(kgrid.y, [], 1), reshape(zeros(kgrid.Nx, kgrid.Ny), [], 1));
   F_interp.V = reshape(source.p0, [], 1);
   i3 = F_interp(sensor.mask(1, :), sensor.mask(2, :));
   toc
   Elapsed time is 2.090216 seconds.

   In 3D:

   tic;
   [tri, bc] = gridDataFast3D(kgrid.x, kgrid.y, kgrid.z, sensor.mask(1,:), sensor.mask(2,:), sensor.mask(3,:));
   i1 = sum(source.p0(tri) .* bc,2);
   toc
   Elapsed time is 98.273495 seconds.

   tic;
   F_interp = TriScatteredInterp(reshape(kgrid.x, [], 1), reshape(kgrid.y, [], 1), reshape(kgrid.z, [], 1), reshape(zeros(kgrid.Nx, kgrid.Ny, kgrid.Nz), [], 1));
   F_interp.V = reshape(source.p0, [], 1);
   i2 = F_interp(sensor.mask(1,:), sensor.mask(2,:), sensor.mask(3,:));
   toc
   Elapsed time is 135.914373 seconds.

   Best,
   Chao

   Posted 12 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Chao,

   Thanks for your functions, they look great! I did some quick testing and agree they seem to be significantly faster, both for precomputing the triangulation, and for evaluating the interpolation. They also avoid using the tsearch function used by the original gridDataFast (which is slowly being phased out of MATLAB), and seem to use a lot less memory than TriScatteredInterp. So all-in-all, I think they'd make a useful addition to the toolbox.

   A possible modification would be to catch any NANS inside t which will occur if the evaluation points are outside the grid dimensions. For example, something like

   if any(isnan(t))
   error('Cartesian sensor coordinates must lie within the grid dimensions');
   end

   Brad.

   Posted 12 years ago #

RSS feed for this topic

Reply

You must log in to post.