2
$\begingroup$

I have six 1D continuous functions that I want to use in order to find a 2D function. Each of these 1D functions has to be placed at a unique angle $\phi$ w.r.t the $x$-axis as given in the image below. In the image, for example, a 1D function in black color is at an angle $\phi = -1.54337$ rad and the red colored 1D function is at an angle $\phi = 0.5$ rad.

enter image description here

This way $n$ 1D functions can be used. Here, I have used six different 1D functions and have placed them at different angles. I'm quite new to Mathematica, therefore, my code is a little messy. The Edited code to construct the data is here:

OneDimFunc1 = 0.03136086274687064` (30.953065429366962`E^(-25.926092562282516` \ (-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.8693935369154184` (0.195211521347951` + xdash)^2));
OneDimFunc2 = 0.03136086274687064` (30.953065429366962`E^(-22 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2 (0.195211521347951` + xdash)^2));
OneDimFunc3 = 0.03136086274687064` (30.953065429366962` E^(-15 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.5 (0.195211521347951` + xdash)^2));
OneDimFunc4 = 0.03136086274687064` (30.953065429366962` E^(-27 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-2.9 (0.195211521347951` + xdash)^2));
OneDimFunc5 = 0.03136086274687064` (30.953065429366962` E^(-17 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-1.8 (0.195211521347951` + xdash)^2));
OneDimFunc6 = 0.03136086274687064` (30.953065429366962` E^(-23 \(-0.0017564449181118602` + xdash)^2) + 0.9338124308343011` E^(-3.2 (0.195211521347951` + xdash)^2));

NumData = 100;         (* Number of data points *)
xD = Array[# &, NumData, {-1, 1}];

zDt1 = Flatten[Abs[OneDimFunc1] /. {xdash -> xD}]; 
x2Ddata1 = Array[# &, NumData, {-1, 1}] Cos[-1.54337];  (* \[Phi] = -1.54337*)
y2Ddata1 = Array[# &, NumData, {-1, 1}] Sin[-1.54337];
Zdt1 = Transpose[{x2Ddata1, y2Ddata1, zDt1}]; 

zDt2 = Flatten[Abs[OneDimFunc2] /. {xdash -> xD}];  
x2Ddata2 = Array[# &, NumData, {-1, 1}] Cos[-1];   (* \[Phi] = 1*)
y2Ddata2 = Array[# &, NumData, {-1, 1}] Sin[-1];
Zdt2 = Transpose[{x2Ddata2, y2Ddata2, zDt2}];

zDt3 = Flatten[Abs[OneDimFunc3] /. {xdash -> xD}]; 
x2Ddata3 = Array[# &, NumData, {-1, 1}] Cos[-0.5];  (* \[Phi] = -0.5*)
y2Ddata3 = Array[# &, NumData, {-1, 1}] Sin[-0.5];
Zdt3 = Transpose[{x2Ddata3, y2Ddata3, zDt3}];

zDt4 = Flatten[Abs[OneDimFunc4] /. {xdash -> xD}];          
x2Ddata4 = Array[# &, NumData, {-1, 1}] Cos[-0];   (* \[Phi] = 0*)
y2Ddata4 = Array[# &, NumData, {-1, 1}] Sin[-0]; 
Zdt4 = Transpose[{x2Ddata4, y2Ddata4, zDt4}];

zDt5 = Flatten[Abs[OneDimFunc5] /. {xdash -> xD}];     
x2Ddata5 = Array[# &, NumData, {-1, 1}] Cos[0.5];  (* \[Phi] = 0.5*)
y2Ddata5 = Array[# &, NumData, {-1, 1}] Sin[0.5];
Zdt5 = Transpose[{x2Ddata5, y2Ddata5, zDt5}];

zDt6 = Flatten[Abs[OneDimFunc6] /. {xdash -> xD}];    
x2Ddata6 = Array[# &, NumData, {-1, 1}] Cos[1];  (* \[Phi] = 1*)
y2Ddata6 = Array[# &, NumData, {-1, 1}] Sin[1]; 
Zdt6 = Transpose[{x2Ddata6, y2Ddata6, zDt6}];

TwoDimdata = Join[Zdt1, Zdt2, Zdt3, Zdt4, Zdt5, Zdt6];  (* Combined data points *)
ListPointPlot3D[TwoDimdata]

After combining the data I fit 2D Gaussian functions to the data using nonlinear fitting as follows:

NBasisFunc2D = 3;
TwoDimModel = Sum[a[i] Exp[-( ((xout - center1[i])/b[i])^2 + ((yout - center2[i])/b[i])^2)], {i, 1, NBasisFunc2D, 1}];  (* Bivariate Gaussians *)
TwoDimParameters = Flatten[Transpose[{Join[Array[a, NBasisFunc2D][[;; ;; 1]], Array[b, NBasisFunc2D][[;; ;; 1]], Array[center1, NBasisFunc2D][[;; ;; 1]], Array[center2, NBasisFunc2D][[;; ;; 1]]]}]];
nlm2D = NonlinearModelFit[TwoDimdata, TwoDimModel, TwoDimParameters, {xout, yout}, MaxIterations -> 1000];

FittedPrams = nlm2D["BestFitParameters"]      

Show[Plot3D[TwoDimModel /. FittedPrams, {xout, -0.45, 0.45}, {yout, -0.45, 0.45}], Graphics3D[{Red, Point /@ TwoDimdata}]]

(Edit: The 2D Gaussian fit is poor after using six different 1D funtions)

I am not sure whether this is a correct way to find out the 2D function. Can I directly use the 1D continuous functions to find out the 2D function?. If yes then how do I do it?. Or should I use cubic spline interpolation instead of least squares fit?.

UPDATE:

I tried the Matlab griddata to interpolate my data. From Matlab documentation: vq = griddata(x,y,v,xq,yq) fits a surface of the form $v = f(x,y)$ to the scattered data in the vectors $(x,y,v)$. The griddata function interpolates the surface at the query points specified by $(xq,yq)$ and returns the interpolated values, $vq$. The Matlab code is:

[xq,yq] = meshgrid(-1:0.01:1);
vq = griddata(x,y,Fxy,xq,yq,'cubic');
surface(xq,yq,vq)
hold on;
plot3(x,y,Fxy,'or','MarkerFaceColor', 'r')
hold off;

where x,y and Fxy is taken from TwoDimdata in Mathematica:

x = TwoDimdata[[All, 1]] ;
y = TwoDimdata[[All, 2]];
Fxy = TwoDimdata[[All, 3]];  

enter image description here

When I try the same in Mathematica I recieve an error:

Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All

How do I do this Matlab equivalent interpolation in Mathematica?

$\endgroup$
  • $\begingroup$ Since your original 1D slice is symmetric and you just want to revolve it around the z-axis in you original plot you just have an expression in polar coordinates where the slice from $(0, \infty)$ of your 1D is $r$ and the $\phi$ is as ususal. $\endgroup$ – b3m2a1 Mar 19 '18 at 17:42
  • $\begingroup$ The code that I have posted here uses just one 1D function (for the sake of simplicity) so maybe you understood that my 1D slice is symmetric. But there are six different 1D functions that I would like to use and those are not symmetric. Therefore, I cannot just revolve it around the $z$-axis. $\endgroup$ – dykes Mar 19 '18 at 19:03
  • $\begingroup$ Do you have any idea how they depend parametrically on $\varphi$? If so you can include that in an analytic expression and then convert from polar to cartesian. Otherwise I think you'll have to fake it via interpolation of some form. $\endgroup$ – b3m2a1 Mar 19 '18 at 19:06
  • 1
    $\begingroup$ You might want to look up transfinite interpolation. $\endgroup$ – J. M. will be back soon Mar 19 '18 at 19:36
  • $\begingroup$ @b3m2a1 I just have these equations on their own terms. I do not have any idea about their parametric dependence on the rotation angle $\phi $. $\endgroup$ – dykes Mar 20 '18 at 0:13
2
$\begingroup$

This is an alternative way to construct your data:

ϕ = {-1.54337, -1, -0.5, 0, .5, 1.};
aa = {-25.926092562282516`, -22, -15, -27, -17, -23};
bb = {-2.8693935369154184`, -2, -2.5, -2.9, -1.8, -3.2};
{α, β, γ, δ, ρ} =  {0.03136086274687064, 30.953065429366962`, -0.0017564449181118602`, 
   0.9338124308343011`, 0.195211521347951`};
cc = Transpose[{ϕ, aa, bb}];
NumData = 100;
ClearAll[oneDfuncs]
oneDfuncs[x_] := {x Cos@#, x Sin@#, Abs[α(β E^(#2 (γ + x)^2) + δ E^(#3 (ρ + x)^2))]}&@@@cc;
twoDdata = Transpose[oneDfuncs /@ Range[-1, 1, 2/99]];

Show[ParametricPlot3D[Evaluate@oneDfuncs[x], {x, -1, 1}, 
  PlotLegends -> ("OneDimFunc" <> ToString[#] & /@ Range[6])], 
 ListPointPlot3D[twoDdata]]

enter image description here

The data generated is the same as your TwoDimdata

Join @@ twoDdata == TwoDimdata

True

$\endgroup$
  • $\begingroup$ How do I interpolate twoDdata with InterpolationOrder-> 3?. It shows an error that " Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All" $\endgroup$ – dykes Mar 22 '18 at 11:43
  • $\begingroup$ @dykes, i don't know how to do interpolation on unstructured grids. Maybe some of these Q/As may be useful. $\endgroup$ – kglr Mar 22 '18 at 20:34
  • $\begingroup$ It's almost working : Show[ParametricPlot3D[Evaluate@oneDfuncs[x], {x, -1, 1}], ListSurfacePlot3D[ Flatten[Transpose[(oneDfuncs /@ Range[-1, 1, 2/99])], 1], MaxPlotPoints -> 80, BoxRatios -> {1, 1, 1/2}, ColorFunction -> Function[{x, y, z}, Hue[z]]]] $\endgroup$ – Mariusz Iwaniuk Mar 22 '18 at 23:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.