Faster method for polyfit along 3rd dimension of a large 3D matrix?
9 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Justin
el 9 de Ag. de 2012
Comentada: Jan
el 11 de Feb. de 2022
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
2 comentarios
Walter Roberson
el 9 de Ag. de 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.
Respuesta aceptada
Teja Muppirala
el 10 de Ag. de 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]);
4 comentarios
Pavel Psota
el 10 de Feb. de 2022
+1, cool, indeed! Could the fitting coefficients be obtained from your solution?
Jan
el 11 de Feb. de 2022
A small improvement is to avoid the expensive power operation:
% Replace:
V = bsxfun(@power,x,0:z);
% by:
V = [ones(300, 1), cumprod(repmat(x, 1, z), 2)];
Más respuestas (1)
Martin Offterdinger
el 9 de Abr. de 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
1 comentario
Jan
el 9 de Abr. de 2019
Editada: Jan
el 9 de Abr. de 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...)
Ver también
Categorías
Más información sobre Logical en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!