k-Wave User Forum » Topic: error about 'CartInterp' and 'DataCast'

k-Wave User Forum » Topic: error about 'CartInterp' and 'DataCast' http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:08:03 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "error about 'CartInterp' and 'DataCast'" http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast#post-449 Thu, 12 Apr 2012 00:19:15 +0000 Bradley Treeby 449@http://www.k-wave.org/forum/ 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 <code>tsearch</code> function used by the original <code>gridDataFast</code> (which is slowly being phased out of MATLAB), and seem to use a lot less memory than <code>TriScatteredInterp</code>. 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 <code>t</code> 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. huangchao on "error about 'CartInterp' and 'DataCast'" http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast#post-448 Wed, 11 Apr 2012 21:27:45 +0000 huangchao 448@http://www.k-wave.org/forum/ 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 huangchao on "error about 'CartInterp' and 'DataCast'" http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast#post-392 Fri, 16 Mar 2012 23:08:14 +0000 huangchao 392@http://www.k-wave.org/forum/ 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 Bradley Treeby on "error about 'CartInterp' and 'DataCast'" http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast#post-187 Thu, 30 Jun 2011 23:49:41 +0000 Bradley Treeby 187@http://www.k-wave.org/forum/ Hi Chao, Unfortunately this is a limitation of the MATLAB function <code>TriScatteredInterp</code> which is used for linear interpolation within k-Wave. This function only supports doubles. For example, try <pre><code>>> 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.</code></pre> In the 2D case, the k-Wave code checks if the parameter <code>'DataCast'</code> is used, and if so it switches to using <code>gridDataFast</code> instead of <code>TriScatteredInterp</code>. However, <code>gridDataFast</code> 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 <a href="http://www.k-wave.org/documentation/example_ivp_3D_simulation.php#heading6">here</a>. Kind regards, Brad. huangchao on "error about 'CartInterp' and 'DataCast'" http://www.k-wave.org/forum/topic/error-about-cartinterp-and-datacast#post-186 Thu, 30 Jun 2011 21:13:41 +0000 huangchao 186@http://www.k-wave.org/forum/ 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