Polynomial fitting of each pixel in an image without looping

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?

Answer 1

[ChristianW]

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

[charlotte]

Thank you, modifirf it a bit and it works well!

Answer 2

[Matt J]

I think looping is your best bet.

Answer 3

[Image Analyst]

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};

[charlotte]

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

Answer 4

[ChristianW]

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

[charlotte]

The over all result is 5 scalar values for each pixels, for 256x256 pixels

[charlotte]

The dimension of each C-element is 256x256x7, your code gives error: index exceeds matrix dimensions

[ChristianW]

Ok, I misunderstood alot here. I tought your C = cell(3,7) with 256x256 matrix each in it.

And my result are 5 scalar k values.

[charlotte]

Maybe I was unclear.. of course i can rearrange my data so I have cell(3,7) with 256x256, but i the result im looking for is 5 scalar k values for each pixel (256x256 pixels)

[Matt J]

Even when correctly formulated as a simultaneous system of equations across all pixels, I'm pessimistic that this will give you better performance than a for-loop. By combining into a single system, you lose simplifying information about how the system can be decomposed into smaller problems. I could be wrong, though.

[charlotte]

ok, thank you. The problem that is that it takes about 256 seconds.

[ChristianW]

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.

Answer 5

[Teja Muppirala]

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

[charlotte]

Thank you!