## Polynomial fitting of each pixel in an image without looping

on 10 Feb 2013

### ChristianW (view profile)

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?

## Products

### ChristianW (view profile)

on 11 Feb 2013
Edited by ChristianW

### ChristianW (view profile)

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
```

charlotte

### charlotte (view profile)

on 13 Feb 2013

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

### Matt J (view profile)

on 10 Feb 2013

I think looping is your best bet.

### Image Analyst (view profile)

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

charlotte

### charlotte (view profile)

on 11 Feb 2013

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 (view profile)

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
```

Matt J

### Matt J (view profile)

on 11 Feb 2013

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

### charlotte (view profile)

on 11 Feb 2013

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

ChristianW

### ChristianW (view profile)

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 (view profile)

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
```

charlotte

### charlotte (view profile)

on 13 Feb 2013

Thank you!

#### Join the 15-year community celebration.

Play games and win prizes!

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi