Thread Subject: 3d surface fitting

Hi,

I have a 320x240 pixel grayscale image, which captures the
intensity profile information of a beam. The beam is a
nondiffracting bessel beam... and hence, features a bessel
function distribution in 3D.

I'm trying to fit a 5 parameter bessel function to the data,
so that I can extract the bessel parameters. While I've done
some 2d curve fitting before, I'm not sure how to do 3D
surface fitting.

I tried using lsqcurvefit and some other functions, but
nothing seems to work. Mainly, I'm not sure how to pass the
data and the objective function correctly.

My 5 parameter objective function is:

a * [besselj(0,b * sqrt( (x-c)^2 + (y-d)^2 ) )]^2 + e

My image data is a 320x240 array of type 'double'.

Can somebody please guide me about fitting my data to the
equation, and producing a mesh of it?

Thanks in advance!
Arjun

"Arjun Chennu" <arjun.chennu@gmail.com> wrote in message
<fum3b2$nc9$1@fred.mathworks.com>...
> Hi,
>
> I have a 320x240 pixel grayscale image, which captures the
> intensity profile information of a beam. The beam is a
> nondiffracting bessel beam... and hence, features a bessel
> function distribution in 3D.
>
> I'm trying to fit a 5 parameter bessel function to the data,
> so that I can extract the bessel parameters. While I've done
> some 2d curve fitting before, I'm not sure how to do 3D
> surface fitting.
>
> I tried using lsqcurvefit and some other functions, but
> nothing seems to work. Mainly, I'm not sure how to pass the
> data and the objective function correctly.
>
> My 5 parameter objective function is:
>
> a * [besselj(0,b * sqrt( (x-c)^2 + (y-d)^2 ) )]^2 + e
>
> My image data is a 320x240 array of type 'double'.
>
> Can somebody please guide me about fitting my data to the
> equation, and producing a mesh of it?

[x,y] = meshgrid(1:240,1:320);
fun = @(c) c(!)*besselj(0,c(2)*sqrt((x-c(3)).^2+(y-c(4)).^2)).^2+c(5);

John

> [x,y] = meshgrid(1:240,1:320);
> fun = @(c)
c(!)*besselj(0,c(2)*sqrt((x-c(3)).^2+(y-c(4)).^2)).^2+c(5);
>
> John


Hi John,

Thanks a lot for your quick reply!

Indeed, I have defined my function as:

function err = myfit(p, realdata)
[x,y] = meshgrid(1:240,1:320);
calcval = p(1)*besselj(0,
p(2)*sqrt((x-p(3)).^2+(y-p(4)).^2)).^2 +p(5);
size(calcval)
err = (calcval - realdata).^2;
return

However, when I use lsqcurvefit:

[p, resnorm] = lsqcurvefit(@myfit,p_start,dd)

it says it needs a fourth argument. (dd is the experimental
data to which I'm trying to fit my function myfit).

The help file says:
x = lsqcurvefit(fun,x0,xdata,ydata) starts at x0 and finds
coefficients x to best fit the nonlinear function
fun(x,xdata) to the data ydata (in the least-squares sense).
ydata must be the same size as the vector (or matrix) F
returned by fun.

so here ydata is dd.

But what is xdata? I'm not sure what I'm supposed to pass as
xdata.

Help muchos appreciated! Thanks again.

Arjun

"Arjun Chennu" <arjun.chennu@gmail.com> wrote in message
<fun86d$o8a$1@fred.mathworks.com>...
>
> > [x,y] = meshgrid(1:240,1:320);
> > fun = @(c)
> c(!)*besselj(0,c(2)*sqrt((x-c(3)).^2+(y-c(4)).^2)).^2+c(5);
> >
> > John
>
>
> Hi John,
>
> Thanks a lot for your quick reply!
>
> Indeed, I have defined my function as:
>
> function err = myfit(p, realdata)
> [x,y] = meshgrid(1:240,1:320);
> calcval = p(1)*besselj(0,
> p(2)*sqrt((x-p(3)).^2+(y-p(4)).^2)).^2 +p(5);
> size(calcval)
> err = (calcval - realdata).^2;
> return
>
> However, when I use lsqcurvefit:
>
> [p, resnorm] = lsqcurvefit(@myfit,p_start,dd)
>
> it says it needs a fourth argument. (dd is the experimental
> data to which I'm trying to fit my function myfit).
>
> The help file says:
> x = lsqcurvefit(fun,x0,xdata,ydata) starts at x0 and finds
> coefficients x to best fit the nonlinear function
> fun(x,xdata) to the data ydata (in the least-squares sense).
> ydata must be the same size as the vector (or matrix) F
> returned by fun.
>
> so here ydata is dd.
>
> But what is xdata? I'm not sure what I'm supposed to pass as
> xdata.
>
> Help muchos appreciated! Thanks again.
>
> Arjun
>

There are two things to do here. First, Do NOT
square the residuals. lsqcurvefit will subtract
the aim (ydata) and square the residuals itself.

Next, you can probably save some small amount
of time by not generating the meshgrid for every
call into your function. Since you must pass in
xdata anyway, this should already contain x
and y.

Using an anonymous function, the call might
look like this:

[x,y] = meshgrid(1:240,1:320);
xdata = {x,y};
fun = @(c,xdata) c(1)*besselj(0, ...
     c(2)*sqrt((xdata{1}-c(3)).^2+(xdata{2}-c(4)).^2)).^2+c(5);

[p, resnorm] = lsqcurvefit(fun,p_start,xdata,dd);

John

> Using an anonymous function, the call might
> look like this:
>
> [x,y] = meshgrid(1:240,1:320);
> xdata = {x,y};
> fun = @(c,xdata) c(1)*besselj(0, ...
>
c(2)*sqrt((xdata{1}-c(3)).^2+(xdata{2}-c(4)).^2)).^2+c(5);
>
> [p, resnorm] = lsqcurvefit(fun,p_start,xdata,dd);
>
> John


Ah, I see. I did not know of this "cell array" option with
the curly braces. Thanks for that tip. :)

Can you tell me if there is some way I can set the lb and ub
(bounds) so that I can constrain them to real numbers?
Complex values in the parameters are troubling.

I've got my lsqcurvefit line trying to come up with a set of
5 parameters now, however the final result is rather far
from the data that I have. (The bessel is too thin). I'm
trying reducing the tolerance to 1e-9, but don't think it
has helped. Anything I might be missing?

Arjun

"Arjun Chennu" <arjun.chennu@gmail.com> wrote in message
<funs5q$jaf$1@fred.mathworks.com>...
> > Using an anonymous function, the call might
> > look like this:
> >
> > [x,y] = meshgrid(1:240,1:320);
> > xdata = {x,y};
> > fun = @(c,xdata) c(1)*besselj(0, ...
> >
> c(2)*sqrt((xdata{1}-c(3)).^2+(xdata{2}-c(4)).^2)).^2+c(5);
> >
> > [p, resnorm] = lsqcurvefit(fun,p_start,xdata,dd);
> >
> > John
>
>
> Ah, I see. I did not know of this "cell array" option with
> the curly braces. Thanks for that tip. :)
>
> Can you tell me if there is some way I can set the lb and ub
> (bounds) so that I can constrain them to real numbers?
> Complex values in the parameters are troubling.

No. You cannot constrain it in this way.

But what you can do is to recognize that
besselj(0,x) has essentially zero imaginary
part over the entire real line. But, when x
is negative, it generates a tiny imaginary
part.

besselj(0,5)
ans =
      -0.1776

besselj(0,-5)
ans =
      -0.1776 - 4.5591e-17i

besselj(0,-25)
ans =
     0.096267 + 1.0399e-17i
besselj(0,0)

So you can ignore that imaginary part, by
inserting a call to real on the result.

fun2 = @(c,xdata) real(fun(c,xdata));
[p, resnorm] = lsqcurvefit(fun2,p_start,xdata,dd);


> I've got my lsqcurvefit line trying to come up with a set of
> 5 parameters now, however the final result is rather far
> from the data that I have. (The bessel is too thin). I'm
> trying reducing the tolerance to 1e-9, but don't think it
> has helped. Anything I might be missing?

If the curve has the wrong fundamental shape,
then insisting that it take on the shape you
want will not help.

Effectively, this is all you are doing when you
tighten the tolerance. You are demanding
more and more fervently that the curve takes
on a different shape. Your demands will fall
on deaf ears I'm afraid.

You might want to take a read through this
document. It talks about the idea of a
fundamental shape of a curve, and what
you can do.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=10864&objectType=FILE

John

> If the curve has the wrong fundamental shape,
> then insisting that it take on the shape you
> want will not help.
>
> Effectively, this is all you are doing when you
> tighten the tolerance. You are demanding
> more and more fervently that the curve takes
> on a different shape. Your demands will fall
> on deaf ears I'm afraid.
>
> You might want to take a read through this
> document. It talks about the idea of a
> fundamental shape of a curve, and what
> you can do.

Thanks! That was indeed an interesting read on data
modelling... cleared up a few concepts better! :-)

However, from my theoretical predictions the intensity
profile must be the square of a bessel function (of the
first kind). That is the model I'm applying...

... and the fit does find a bessel, but not centered at the
same spot as my empirical curve, and neither the same
'thickness' of the central core.

Thus: a huge resnorm (of 1e7 or 1e8 order), and a residual
with bumps and humps from the zero plane.


While I realize I'm trying to coerce the minimization to
render a bessel form, there ARE some problems with my data:

1. In some of the frames, the central core has saturated the
ccd, and so the "heads" of the bessel are chopped off at the
ends.

2. Due to a slight amount of astigmatism, the annular rings
are not perfectly symmetrical but have a higher intensity in
one direction, or so.


I thought these limitations would not change severly affect
the fit, because the smallest residual should be the one
with the bessel's "head" filled out, with the annular rings
matching up.

But that doesn't seem to be the case. Is there an
intelligent way to edit my model to incorporate the planar
top of the bessel into the fit?

This has been a journey into modelling for me. Thanks for
the guidance! :)

Arjun

"Arjun Chennu" <arjun.chennu@gmail.com> wrote in message
<fur8d1$8ja$1@fred.mathworks.com>...
>
> > If the curve has the wrong fundamental shape,
> > then insisting that it take on the shape you
> > want will not help.
> >
> > Effectively, this is all you are doing when you
> > tighten the tolerance. You are demanding
> > more and more fervently that the curve takes
> > on a different shape. Your demands will fall
> > on deaf ears I'm afraid.
> >
> > You might want to take a read through this
> > document. It talks about the idea of a
> > fundamental shape of a curve, and what
> > you can do.
>
> Thanks! That was indeed an interesting read on data
> modelling... cleared up a few concepts better! :-)
>
> However, from my theoretical predictions the intensity
> profile must be the square of a bessel function (of the
> first kind). That is the model I'm applying...
>
> ... and the fit does find a bessel, but not centered at the
> same spot as my empirical curve, and neither the same
> 'thickness' of the central core.
>
> Thus: a huge resnorm (of 1e7 or 1e8 order), and a residual
> with bumps and humps from the zero plane.
>
>
> While I realize I'm trying to coerce the minimization to
> render a bessel form, there ARE some problems with my data:
>
> 1. In some of the frames, the central core has saturated the
> ccd, and so the "heads" of the bessel are chopped off at the
> ends.
>
> 2. Due to a slight amount of astigmatism, the annular rings
> are not perfectly symmetrical but have a higher intensity in
> one direction, or so.
>
>
> I thought these limitations would not change severly affect
> the fit, because the smallest residual should be the one
> with the bessel's "head" filled out, with the annular rings
> matching up.
>
> But that doesn't seem to be the case. Is there an
> intelligent way to edit my model to incorporate the planar
> top of the bessel into the fit?
>
> This has been a journey into modelling for me. Thanks for
> the guidance! :)
>
> Arjun

Modeling is a sometimes tricky thing, as you
are finding. It sounds like your data has some
problems. Why does data always come out that
way? ;-)

One option is to down-weight data points that
you don't trust. Its not that hard to build a
weighted fit.

Do you have good enough starting values?
Sometimes if the starting values are poor on
problematic data, the whole thing diverges on
you. One way to improve things is to use a
partitioned fit. The idea is to use a tool like
my fminspleas, which reduces the problem
to a 3 variable one from 5 variables. So you
only need to choose starting values for those
three parameters.

data = {x,y};
funlist = {@(params,data) besselj(0, ...
   params(1)*sqrt( (data{1}-params(2)).^2 + ...
   (data{2}-params(3)).^2 )).^2,1};

[INLP,ILP] = fminspleas(funlist,NLPstart,data,z)

The nonlinear parameters are [b,c,d] in your
model, and the linear parameters are [a,e].

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
objectId=10093&objectType=FILE

I'm not certain about what the problem is with
your data. Is it the asymmetry? Or is it that
the central peak has the wrong shape?

You might try down weighting the outer lobes,
or vice-versa. Or perhaps build an asymmetric
bessel, so that the radial behavior is a function
of angle. Or, you might add a secondary term
to tweak the shape of that central peak. There
are lots of things you might try here, so don't
give in.

John

> Modeling is a sometimes tricky thing, as you
> are finding. It sounds like your data has some
> problems. Why does data always come out that
> way? ;-)

I don't know, but there should be a statistical law about
the errors of experimentation.... as oppossed to just "being
careful". :)

I'm really thankful for your (yet another) tip about
fminspleas. Reducing it to a 3 NLP fit already starting
giving me the right bessel form, but still couldn't do it
well for half the datasets.

After that I used the built-in weighting options to focus
mainly on the central lobe, and now the fits work for
(almost) all the datasets! :-)

I'm not sure I really understand though what fminspleas does
different from lsqcurvefit. They both seem to solve problems
in a least square sense, but the main difference seems to be
in the way the parameters are handled/prioritized.

However, I couldn't really understand why it is faster or
more effective?

Thanks again! :) A real pleasure to talk with an 'expert'.

"Arjun Chennu" <arjun.chennu@gmail.com> wrote in message
<fv29eq$38p$1@fred.mathworks.com>...
> > Modeling is a sometimes tricky thing, as you
> > are finding. It sounds like your data has some
> > problems. Why does data always come out that
> > way? ;-)
>
> I don't know, but there should be a statistical law about
> the errors of experimentation.... as oppossed to just "being
> careful". :)
>
> I'm really thankful for your (yet another) tip about
> fminspleas. Reducing it to a 3 NLP fit already starting
> giving me the right bessel form, but still couldn't do it
> well for half the datasets.
>
> After that I used the built-in weighting options to focus
> mainly on the central lobe, and now the fits work for
> (almost) all the datasets! :-)

(Great!)

 
> I'm not sure I really understand though what fminspleas does
> different from lsqcurvefit. They both seem to solve problems
> in a least square sense, but the main difference seems to be
> in the way the parameters are handled/prioritized.
>
> However, I couldn't really understand why it is faster or
> more effective?

lsqcurvefit must search a parameter space of
five dimensions on your problem.

When fminspleas does the optimization, it
breaks the problem into two different classes
of parameters. Two of the parameters can be
estimated using a simple linear regression,
so fminspleas does not need to iterate over
that portion of the parameter space. These
are the "linear" parameters. In fact, a better
description is "conditionally linear". The
reduction of the problem dimensionality from
five to three dimensions is far more significant
than is the fact that fminsearch is not really a
terribly efficient optimizer.

John