# Polynomial fitting of each pixel in an image without looping

5 views (last 30 days)

Show older comments

##### 0 Comments

### Accepted Answer

ChristianW
on 11 Feb 2013

Edited: ChristianW
on 13 Feb 2013

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

### More Answers (4)

Image Analyst
on 10 Feb 2013

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

ChristianW
on 11 Feb 2013

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 Comments

ChristianW
on 11 Feb 2013

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
on 12 Feb 2013

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

### See Also

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!