Polynomial fitting of each pixel in an image without looping

3 次查看(过去 30 天)
I have 3x7 images of 256x256 pixels stored in a cell array, i.e. for each pixel i have 7 x-values, 7 y-values and 7 z-values. I want to find the coefficients for z=k1*x + k2*y + k3*x^2 + k4*y^2 + k5*x*y in a least square sense for each pixel without looping over each pixel. Is there a more efficient way to do this?

采纳的回答

ChristianW
ChristianW 2013-2-11
编辑:ChristianW 2013-2-13
Referring to your 256 sec calulation time:
Got it to 1 sec and 0.8 with parfor. (on my cpu)
dim = 256;
C = mat2cell(randi(255,dim*3,dim*7), dim*ones(1,3), dim*ones(1,7));
tic
C0 = cellfun(@(x) reshape(x,1,[]),C,'UniformOutput',false);
X = cat(1,C0{1,:});
Y = cat(1,C0{2,:});
Z = cat(1,C0{3,:});
K = cell(size(C{1}));
for i = 1:size(X,2) % 1:NumberOfPixelsPerImage
K{i} = [X(:,i), Y(:,i), X(:,i).^2, Y(:,i).^2, X(:,i).*Y(:,i)]\Z(:,i);
end
toc

更多回答(4 个)

Matt J
Matt J 2013-2-10
I think looping is your best bet.

Image Analyst
Image Analyst 2013-2-10
I don't understand your data layout. So you have a cell array with 3 rows and 7 columns. What is inside each cell? Is each cell a 256 by 256 array of either x, y, or z values, like
{256x256 x1 values}, {256x256 x2 values},....{256x256 x7 values};
{256x256 y1 values}, {256x256 y2 values},....{256x256 y7 values};
{256x256 z1 values}, {256x256 z2 values},....{256x256 z7 values};
  1 个评论
charlotte
charlotte 2013-2-11
The data is stored as you described in your code: {256x256 x1 values}, {256x256 x2 values},....{256x256 x7 values}; {256x256 y1 values}, {256x256 y2 values},....{256x256 y7 values}; {256x256 z1 values}, {256x256 z2 values},....{256x256 z7 values};.
I want to find the least square solution for the equations for each pixel: z1=k1*x1 + k2*y1 + k3*x1^2 + k4*y1^2 +k5*x1*y1 z2=k1*x2 + k2*y2 + k3*x2^2 + k4*y2^2 +k5*x2*y1 ... z7=k1*x7 + k2*y7 + k3*x7^2 + k4*y7^2 +k5*x7*y7

请先登录,再进行评论。


ChristianW
ChristianW 2013-2-11
Is the overall result just 5 scalar k values?
X = cat(1,C{1,:});
Y = cat(1,C{2,:});
Z = cat(1,C{3,:});
M = [X(:), Y(:), X(:).^2, Y(:).^2, X(:).*Y(:)];
K = M\Z(:); % Z = M*K
  7 个评论
charlotte
charlotte 2013-2-11
ok, thank you. The problem that is that it takes about 256 seconds.
ChristianW
ChristianW 2013-2-11
I'll give it a second shot. I need some help with the math.
Lets talk about a single pixel only. With 7 xValues in X(:,1), each row one of the 7 pictures (analogously for Y and Z), like this:
X = [pixel1_image1
pixel1_image2
...
pixel1_image7];
With these inputs, does this function solve the equations for that pixel?
function K = fcn(X,Y,Z)
M = [X, Y, X.^2, Y.^2, X.*Y];
K = M\Z; % Z = M*K
I need a check for that function or an example input with correct output to varify.

请先登录,再进行评论。


Teja Muppirala
Teja Muppirala 2013-2-12
Here's an approach using sparse matrices to do it. this works in about 0.3 seconds for me, and gives the coefficients in a 5x65536 matrix K.
It should be noted that using a simple for-loop is much simpler to implement, and still works in about 0.6 seconds if you preallocate properly.
% Making random data
dim = 256;
C = mat2cell(randi(255,dim*3,dim*7), dim*ones(1,3), dim*ones(1,7));
tic
C0 = cellfun(@(x) reshape(x,1,[]),C,'UniformOutput',false);
X = cat(1,C0{1,:});
Y = cat(1,C0{2,:});
Z = cat(1,C0{3,:});
M = permute( cat(3,X,Y,X.^2,Y.^2,X.*Y), [1 3 2]);
% Generate the locations of the block-diagonal sparse entries
jloc = repmat(1:(dim^2*5),7,1);
iloc = bsxfun(@plus, repmat((1:7)',1,5) ,reshape( 7*(0:dim^2-1) , 1, 1, []));
SM = sparse(iloc(:),jloc(:),M(:));
% Do the pseudoinversion
K = (SM'*SM) \ (SM'*Z(:));
K = reshape(K,5,[]);
toc

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by