Asked by Justin
on 9 Aug 2012

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

Answer by Teja Muppirala
on 10 Aug 2012

Accepted Answer

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
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
on 18 Dec 2017

+1, Teja, what a powerful solution!

Sign in to comment.

Answer by Martin Offterdinger
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

@Martin: Please do not attach a new question in the section for answer. Open a new thread instead and remove this pseudo-answer. Including a link to this thread is a good idea. Thanks.

What's wrong with setting z=1? Which array is "the same" and why do you need to determine the fit multiple times for the same array? (Please explain this in your new question...)

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Walter Roberson (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/45680-faster-method-for-polyfit-along-3rd-dimension-of-a-large-3d-matrix#comment_93812

## Justin (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/45680-faster-method-for-polyfit-along-3rd-dimension-of-a-large-3d-matrix#comment_94457

Sign in to comment.