Message Boards

WOLFRAM COMMUNITY

10817 Views

2 Replies

1 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Image Processing Mathematics Optics Physics Wolfram Language Modeling

fitting the Airy diffraction pattern with FindFit

Alex Mikh

Posted 11 years ago

Hello! I have 2D data (images) which are results of Gaussian type beam diffraction. According to the theory of diffraction the pattern can be discribed by function ~ J1(r)/r. I tried to fit my data with this function using the FIndFit routine, but had an error : Blockquote Power::infy: "Infinite expression 1/0. encountered" Blockquote data01 = Import["366.0000cm_660.0000nm.png"]; data02 = ImageData[data01]; min02 = Min[data02]; max02 = Max[data02]; (Dimensionalities of images) dim01a = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[1]]; dim01b = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[2]]; dim01 = Min[dim01a, dim01b]; ReducPix = 5; data03 = Table[ (data02[[i]][[j]][[1]] - min02)/max02, {i, 1, dim01b, ReducPix}, {j, 1, dim01a, ReducPix}]; (Dimensionalities) dim03a = Dimensions[data03][[1]]; dim03b = Dimensions[data03][[2]]; dim03 = Min[dim03a, dim03b] ListPlot3D[data03, PlotRange -> All] s[x_, y_] := Sqrt[a2 (x - a1)^2 + b2 (y - b1)^2] Bessel3D = Piecewise[ {{1/4, x == 0 && y == 0}, {(a0 BesselJ[1, s[x, 0]]/s[x, 0]), y == 0}, {(a0 BesselJ[1, s[0, y]]/s[0, y]), x == 0}}, {(a0 BesselJ[1, s[x, y]]/s[x, y])}] data3D = Flatten[MapIndexed[{#2[[1]], #2[[2]], #1} &, data03, {2}], 1]; fit = NonlinearModelFit[ data3D, (Bessel3D)^2, {{a0, 1}, {a1, 50}, {a2, 10}, {b1, 50}, {b2, 10}}, {x, y}]; I tried to use different method in the function (NMinimize, Gradient, so on ) but it did not help. It is really puzzling for me. Will appreciate any suggestions!!

Hello! I have 2D data (images) which are results of Gaussian type beam diffraction. enter image description here According to the theory of diffraction the pattern can be discribed by function ~ J1(r)/r. I tried to fit my data with this function using the FIndFit routine, but had an error :

Blockquote Power::infy: "Infinite expression 1/0. encountered" Blockquote

data01 = Import["366.0000cm_660.0000nm.png"];
data02 = ImageData[data01];
min02 = Min[data02];
max02 = Max[data02];

(*Dimensionalities of images*)
dim01a = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[1]];
dim01b = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[2]];
dim01 = Min[dim01a, dim01b];

ReducPix = 5;
data03 = Table[ (data02[[i]][[j]][[1]] - min02)/max02, {i, 1, dim01b, 
    ReducPix}, {j, 1, dim01a, ReducPix}];

(*Dimensionalities*)
dim03a = Dimensions[data03][[1]];
dim03b = Dimensions[data03][[2]];
dim03 = Min[dim03a, dim03b]

ListPlot3D[data03, PlotRange -> All]
s[x_, y_] := Sqrt[a2 (x - a1)^2 + b2 (y - b1)^2]
Bessel3D = Piecewise[
  {{1/4, x == 0 && y == 0},
   {(a0 BesselJ[1, s[x, 0]]/s[x, 0]), y == 0},
   {(a0 BesselJ[1, s[0, y]]/s[0, y]), x == 0}},
  {(a0 BesselJ[1, s[x, y]]/s[x, y])}]

data3D = Flatten[MapIndexed[{#2[[1]], #2[[2]], #1} &, data03, {2}], 1];
fit = NonlinearModelFit[
   data3D, (Bessel3D)^2, {{a0, 1}, {a1, 50}, {a2, 10}, {b1, 50}, {b2, 
     10}}, {x, y}];

I tried to use different method in the function (NMinimize, Gradient, so on ) but it did not help. It is really puzzling for me. Will appreciate any suggestions!!

POSTED BY: Alex Mikh

2 Replies

Sort By:

David Keith

Posted 11 years ago

Hi Alex, This is an interesting problem. The airy function has a singularity at 0 which can be removed by defining the value there as 1, which is the limit. I do that in the code below, and also make use of the radial symmetry to fit for only the central amplitude and the center. (Which is not all that is needed.) I have found that the code below works intermittently. It seems still dependent on choices of values made during the fit function, even though I think I have removed the singularity. I Hopefully someone with more knowledge will speak up. In[1]:= (2 BesselJ[1, r]/r)^2 /. r -> 0 During evaluation of In[1]:= Power::infy: Infinite expression 1/0^2 encountered. >> During evaluation of In[1]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >> Out[1]= Indeterminate In[2]:= Limit[(2 BesselJ[1, r]/r)^2, r -> 0] Out[2]= 1 In[3]:= airy[r_] := (2 BesselJ[1, r]/r)^2 /; r != 0. In[4]:= airy[0.] = 1 Out[4]= 1 In[5]:= airyData = Flatten[Table[{x, y, 5 airy[Sqrt[(x - 1)^2 + (y - 2)^2]]}, {x, -10, 10, .1}, {y, -10, 10, .1}], 1]; In[6]:= ListPlot3D[airyData]; In[7]:= NonlinearModelFit[airyData, a airy[Sqrt[(x - h)^2 + (y - k)^2]], {a, h, k}, {x, y}]

Hi Alex,

This is an interesting problem. The airy function has a singularity at 0 which can be removed by defining the value there as 1, which is the limit. I do that in the code below, and also make use of the radial symmetry to fit for only the central amplitude and the center. (Which is not all that is needed.)

I have found that the code below works intermittently. It seems still dependent on choices of values made during the fit function, even though I think I have removed the singularity. I Hopefully someone with more knowledge will speak up.

In[1]:= (2 BesselJ[1, r]/r)^2 /. r -> 0

During evaluation of In[1]:= Power::infy: Infinite expression 1/0^2 encountered. >>

During evaluation of In[1]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

Out[1]= Indeterminate

In[2]:= Limit[(2 BesselJ[1, r]/r)^2, r -> 0]

Out[2]= 1

In[3]:= airy[r_] := (2 BesselJ[1, r]/r)^2 /; r != 0.

In[4]:= airy[0.] = 1

Out[4]= 1

In[5]:= airyData = 
  Flatten[Table[{x, y, 5 airy[Sqrt[(x - 1)^2 + (y - 2)^2]]}, {x, -10, 
     10, .1}, {y, -10, 10, .1}], 1];

In[6]:= ListPlot3D[airyData];

In[7]:= NonlinearModelFit[airyData, 
 a airy[Sqrt[(x - h)^2 + (y - k)^2]], {a, h, k}, {x, y}]

POSTED BY: David Keith

Alex Mikh

Posted 11 years ago

Hi David! thanks for your interest to this problem! You are right, there is the singularity at zero point. I tried to eliminate it using a Picewise definiton of the model function in my code but for some reason it did no help and I still had the error. Looks like your definition of the function `airy[r_] := (2 BesselJ[1, r]/r)^2 /; r != 0. airy[0.] = 1;` is able to overcome that issue. Great! An addition specification of the range for parameters in NonlinearModelFit helps as well. However when I'm trying to fit my data the result is not accurate. A new error pops up: NonlinearModelFit::nrlnum: The function value {2.3360210^119-3.4376510^103 I,3.2420410^118+0. I,4.6252810^117+6.1297110^101 I,6.7870110^116-1.18210^101 I,<<44>>,1.1806110^96+0. I,1.15086*10^96+0. I,<<12238>>} is not a list of real numbers with dimensions {12288} at {a0,a1,a2,b1,b2} = {0.206555,69.9985,-2.9634,50.0133,-2.99061}. >> not sure what does it mean.. The 3D plots with parameters from the fitting are not very similiar to the original data at all and as you mentioned very sensitive to the initial values (guess) in the fitting. Thank you anyway! Thanks for your suggestions!

POSTED BY: Alex Mikh

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback