How to interpolate between matrices over time?

Hi all, I have 20 matrices, each of the size 128X128X128. Each matrix contains the temperatures in space, and every matrix is in different time. I want to found out for each pixel on what time the temperature was 10, so I want to get 128X128X128 matrix contains the time when the temperature in that pixel was 10. Of course I can do 128^3 interpolations, but I want it to be efficient.

Answer 1

[SK]

To complete the last interpolation bit, it is as simple as:

R = (M(linhi) - 10)./(M(linhi) - M(linlo));
pixtimes = (1-R).*Time(IXlo) + R.*(Time(IXhi)-Time(IXlo));
pixtimes = reshape(pixtimes, [128, 128, 128]);

Takes under 1 second on my machine.

[Aviad]

Thank you very much! Now I need to add a compilaction: What can I do if for each pixel the temperature I want to interpolate is different? (i.e. for one pixel I want to know the time when the temperature was 11, for another pixel I want to know the time when the temperature was 30 etc.) Thanks again.

[Sean de Wolski]

Loop over all of them and solve the zero finding problem at each point.

[SK]

ALternatively, If T is your matrix of pixel temperatures (128 x 128 x 128), instead of doing

Blo = M < 10;

 you could do:

T = transpose(T(:));
T = ones(20,1)*T;
M = M(:,:)
Blo = (M < T);
Bhi = (M >= T);

Answer 2

[Sean de Wolski]

Well you could do 128^3 optimization problems too. I don't see any reason why you can't do this with interp1.

x = rand(128,128,128,10); % simulate your data stacked in the 4th dimension - time
v = permute(x,[4 1 2 3]); % first dimension is the one we interpolate over

tic
% Do the interpolation
vv = interp1(1:10,v,1:0.05:10);
toc

% permute back so time is 4th dim
vfinal = ipermute(vv,[4 1 2 3]);

On my laptop it's taking about 8s.

Of course this only interpolates to that resolution. If you need better, it might make sense to solve this with fzero for each 128^3 voxel in time. The objective function would be using interp1 on those 10 points. Let me know if this approach interests you.

Answer 3

[SK]

I misunderstood your problem the first time thinking that you wanted to find the temperature at time 10. But you want it the other way around - the time at which the temperature was 10. That is a little tricky

They may be more than one such time for each pixel, no? In that case do you want both times, only the first ...? Anyway supposing that there is exactly one and only one such time, you could do (Assuming that your matrices have been put into 4-D form (20x128x128x128):

M = M(:,:);
Blo = (M < 10);
Bhi = (M >= 10);

[~, J] = max(Blo(end:-1:1, :));
IXlo = 20 - J + 1;

[~, IXhi] = max(Bhi);

N = 128^3;
linlo = 20*(0 : N-1) + IXlo;
linhi = 20*(0 : N-1) + IXhi;

Now linlo contains the linear indices of M of the last temp that is < 10. linhi contains the linear indices in M of the first temp that is >= than 10. Then, assuming that you have the 20 times in a vector, say:

Time = 0.5 : 0.2 : 4.3;

you can do a simple linear interpolation of the 'Time' vector between M(linlo) and M(linhi) to get the N times in a long vector (call it say pixtimes) of length N. Then reshape pixtimes into 128x128x128.