interp2 with huge matrix

Hello

Problem: I have a matrix (M) with 6000 by 11000 elements of complex or double type. The task is to preform a second order polynomial transformation (shift and warp) with interp2().

[X Y] = meshgrid(1:11000,1:6000); % very fast
% A and B are the coefficients
transY = A(1) + A(2)*Y + A(3)*X + A(4)*X.*Y + A(5)*X.^2 + A(6)*Y.^2; % very fast
transX = B(1) + B(2)*Y + B(3)*X + B(4)*X.*Y + B(5)*X.^2 + B(6)*Y.^2; % very fast
ZI = interp2(X,Y,M,transX,transY); % out of memory

The interp2() is very slow and as you can see there are four very large matrices and ZI would be the 5th, so running out of memory is always the case.

How to solve this problem?

So far I have tried to cut the matrix into tiles, but this does not work because of the NaN values for each tile.

Thanks in advance!

Answer 1

[the cyclist]

By "cut the matrix into tiles" do you mean doing the interpolation by sections (such as X=1:1000, then X=1001:1100, etc)? I would expect that to work, if handled properly and being careful about edge cases. Maybe you should post the code you used to try to do that part.

[Michael]

Thanks!

My code is:

function [OUT] = RESAMPLE(DATA)

% this function cuts a 2D Matrix into tiles and applies a +200 shift in Y

% and +300 in X

[imagewidth, imagelength] = size(DATA)

tilewidth = 800;

tilelength = 800;

OUT = zeros([imagewidth, imagelength]);

num_tiles_width = int32(fix(imagewidth/tilewidth));

num_tiles_length = int32(fix(imagelength/tilelength));

restwidth = mod(imagewidth, tilewidth);

restlength = mod(imagelength, tilelength);

if restwidth ~= 0

num_tiles_width + 1

end

if restlength ~= 0

num_tiles_length + 1

end

for i = 0:num_tiles_width

%fprintf('ROW: %i\n',i);

if i == num_tiles_width

%fprintf('STOP\n')

left = (i * tilewidth) + 1;

right= (i * tilewidth) + restwidth;

else

left = (i * tilewidth) + 1;

right= (i+1) * tilewidth;

end

for j = 0:num_tiles_length

if j == num_tiles_length

%fprintf('STOP\n')

top = (j * tilelength) + 1;

bottom = (j * tilelength) + restlength;

else

top = (j * tilelength) + 1;

bottom = (j + 1) * tilelength;

end

fprintf('left: %i \t right: %i \t top: %i \t bottom: %i \n',left,right,top,bottom);

OUT(left:right,top:bottom) = interpolatetile(left,right,top,bottom);

end

end

function [ZI] = interpolatetile(l,r,t,b);

M = (DATA(l:r,t:b));

[meshX meshY] = meshgrid(t:b,l:r);

meshX = double(meshX);

meshY = double(meshY);

M = DATA(l:r,t:b);

%transY = A(1) + A(2)*meshY + A(3)*meshX + A(4)*meshX.*meshY + A(5)*meshX.^2 + A(6)*meshY.^2;

%transX = B(1) + B(2)*meshY + B(3)*meshX + B(4)*meshX.*meshY + B(5)*meshX.^2 + B(6)*meshY.^2;

transY = meshY + 200;

transX = meshX + 300;

M = M + double(r)*double(b);

ZI = interp2(meshX,meshY,M,transX,transY);

end

end

The NaN values of each tile (because of the shift every tile has NaNs) are the problem, and i do not know how to handle that.

Michael

[Ashish Uthama]

Michael, consider editing your question to add this code with proper formatting. Aside: blkproc or blockproc might be useful(?).