## Faster method for polyfit along 3rd dimension of a large 3D matrix?

### Justin (view profile)

on 9 Aug 2012
Latest activity Edited by Jan

on 9 Apr 2019

### Teja Muppirala (view profile)

I am currently working with large data cubes in the form of an MxNxP matrix. The P-dimension represents the signal at each (m,n) pixel. I must obtain a Z-order polynomial fit(where Z varies depending on the situation) for each signal in the data cube. Currently, I utilize a "for" loop to obtain the signal at each pixel, obtain the polynomial coefficients, then calculate the fitted curve. The code I use is fundamentally similar to the following:
dataCube = rand(1000,1000,300);
x = rand(300,1);
sizeCube = size(dataCube);
polyCube = zeros(sizeCube);
for ii = 1:sizeCube(1);
for iii = 1:sizeCube(2);
signal = squeeze(dataCube(ii,iii,:));
a = polyfit(x,signal,z)
y = polyval(a,x);
polyCube(ii,iii,:) = y;
end
end
Because of the quantity of iterations in the for loop, this operation takes a considerable amount of time for each data cube. Is there a faster way to obtain the polynomial fitting, without having to resort to the iterative process I use here. Perhaps, something similar to the filter function where you can apply the filter to a specific dimension of a matrix, rather than extracting each signal?
filteredCube = filter(b,a,dataCube,[],3)
Thanks, Justin

Walter Roberson

### Walter Roberson (view profile)

on 9 Aug 2012
Have you considered re-ordering your data, at least during the processing, so that the dimension you are fitting over is the first dimension? Access (and assignment) over the first dimension is faster.
Justin

### Justin (view profile)

on 14 Aug 2012
Walter,
I attempted reordering the matrix and performing the process as above, over the 1st dimension. There was an incremental improvement in speed, but not nearly the type of advance I was hoping I could achieve.
Check out the answer by Teja below. Its exactly what I was hoping for and works perfectly.
Thanks for the input! Justin

### Teja Muppirala (view profile)

on 10 Aug 2012

This can be accomplished in a fraction of the time with some matrix operations.
dataCube = rand(100,100,300);
sizeCube = size(dataCube);
x = rand(300,1);
z = 3;
V = bsxfun(@power,x,0:z);
M = V*pinv(V);
polyCube = M*reshape(permute(dataCube,[3 1 2]),sizeCube(3),[]);
polyCube = reshape(polyCube,[sizeCube(3) sizeCube(1) sizeCube(2)]);
polyCube = permute(polyCube,[2 3 1]);

Justin

### Justin (view profile)

on 14 Aug 2012
Teja,
Thank you! This is a slick solution that drastically reduces the processing time. It went from ~50 minutes to process a 1920x1040x301 image_cube (@ 0.0015 seconds per loop), to about 15 seconds!
Thank you for the help! Justin
Jan

### Jan (view profile)

on 18 Dec 2017
+1, Teja, what a powerful solution!

### Martin Offterdinger (view profile)

on 9 Apr 2019

Dear Teja,
I am having a similar problem- actually a simpler one even. I have the same array, but I always need to fit a first-order polynom (linear, z=1 in your code). Is it possible to get the coefficients of the linear fit (p1,p2) from your solution as well?
Thanks,
Martin

Jan

on 9 Apr 2019