Copyright (c) [2019] Madeline Hamilton, Pavel Nadolsky (nadolsky@smu.edu)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the “Software”), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

1. The above copyright notice, citation policy, permission notice, and program description 
shall be included in all copies or substantial portions of the Software.

2. You shall cite the following research article on which this program is based in the publications 
that use any results obtained with this program:

Fast Ellipsoidal Fitting of Discrete Multi-Dimensional Data,  
by Rafey Anwar, Madeline Hamilton, Pavel Nadolsky
E-print: arXiv:1901.05511 [physics.data-an] (https://arxiv.org/abs/1901.05511)
SMU Journal of Undergraduate Research: Vol. 5 , Article 4 
https://doi.org/10.25172/jour.5.1.4


Permission notice: THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
 
Program description: Reconstruction of a 3-dimensional ellipsoid from quasi-Gaussian data
Version 2019-01-10
This notebook helps to visualize the ellipsoid reconstruction algorithm developed by R. Anwar, M. Hamilton and P. Nadolsky. 
It generates 3-dimensional quasi-ellipsoidal data and reconstructs an ellipsoid describing the data using the algorithm.

Visual and Convex Hull Setup

This cell sets up the visuals for the rest of the notebook and initializes necessary matrices for the Convex Hull method. This cell must be executed before the fit.

(* Visual setup *)

Switch[StringTake[$Version,{1,3}]//ToExpression//Round,
6,<<"BarCharts`";<<"Histograms`";<<"PieCharts`";<<"ErrorBarPlots`";SetOptions[{Plot,ListPlot,ErrorListPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
7,SetOptions[{Plot,ListPlot,ErrorListPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
8,SetOptions[{Plot,ListPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
9,SetOptions[{Plot,ListPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
10,SetOptions[{Plot,ListPlot,ListLogPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
11,SetOptions[{Plot,ListPlot,ListLogPlot,ListLogLinearPlot,Graphics,ParametricPlot},BaseStyle->{FontFamily->"Helvetica",16}],
_,<< Graphics`Graphics`;<<Graphics`FilledPlot`;
];(* Switch -> *)

(* Initializing matrices necessary for the Convex Hull method *)

CC=(({{0, 0, 2, 0, 0, 0}, {0, -1, 0, 0, 0, 0}, {2, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}}));
{U,L,V}=SingularValueDecomposition[CC];
U2={U[[1]],U[[2]],U[[3]],U[[6]],U[[5]],U[[4]]};
V2={V[[1]],V[[2]],V[[3]],V[[6]],V[[5]],V[[4]]};
CC-U.L.Transpose[V];
CC-U2.L.Transpose[V2];
U.Transpose[U];
V.Transpose[V];
Transpose[U].U;
Transpose[V].V;
cc1=CC[[1;;3,1;;3]];
u1=U[[1;;3,1;;3]];
v1=V[[1;;3,1;;3]];
l=L[[1;;3,1;;3]];
linv=Inverse[l];

Generate and Examine the Input Data

Here, a dataset of 1000 observations from a 3-dimensional Gaussian distribution is produced and visualized. The “true” matrix representing the ellipsoid is calculated. plist contains the data points of the input ellipsoid. MeanX contains the coordinates of the ellipsoid’s center. {a,b,c} are the principal axes of the ellipsoid. The rotation of the ellipsoid with respect to the {x,y,z} is given by the orthogonal matrix OM.

(* Generate the quasi-ellipsoidal data *)
d=3;
Recompute=True;  (* Set Recompute=True to recompute the points of the ellipsoid *)
If[Recompute,
{a,b,c}={1.5,0.5,3}; Npt=1000;
plist0=Table[{φ,θ}={2π RandomReal[],π RandomReal[]};
{a Cos[φ] Sin[θ],b Sin[φ] Sin[θ],c Cos[θ]},{i,Npt}]];

{w1, w2}={{-0.3,0.4,0.35},{0.3,-0.5,0.2}};
OM=RotationMatrix[1.17,w2].RotationMatrix[0.712,w1];
MeanX={0.9,4.3,5.2};
plist=Table[(MeanX+OM.plist0[[ip]])/.List:>LF,{ip,Length[plist0]}];

(* Create a 3D scatter plot of the data *)
ListPointPlot3D[plist/.LF:>List]

(* Calculate the "answer"; the matrix characterizing the ellipsoid *)
aOriginal=OM.DiagonalMatrix[{(1/(a^2)),(1/(b^2)),(1/(c^2))}].Transpose[OM]
(* Visualize the data (the blue points) overlaid on the ellipsoid's true equation (the green mesh) *)
p1=plist/.LF[a__]:>{a}-MeanX//ListPointPlot3D;
p2=ContourPlot3D[1-(({{x1, x2, x3}})).aOriginal.(({{x1}, {x2}, {x3}}))==0,{x1,-3,3},{x2,-3,3},{x3,-3,3},Contours->10,ContourStyle->Directive[Opacity[0],Orange],MeshStyle->RGBColor[0.3,0.9,0],AxesLabel->{"x1","x2","x3"},PlotRange->{{-3,3},{-3,3},{-3,3}}];
Show[p2,p1]

Projection Reconstruction

This cell iterates through the 3 unique orthogonal two-dimensional projections of the data and fits ellipses using both the Covariance Matrix method and the Convex Hull method. Each projection of the data is visualized with their fitted ellipses.

mean=MeanX; 
Print["Reconstruct the ellipsoid's projections:\naCM (covariance matrix, solid red); aHull (convex hull, black dashed)"];
{aCM,aHull}={ConstantArray[0,d*(d-1)/2],ConstantArray[0,d*(d-1)/2]};
Module[{x1,y1,x0,y0,dx,dy,xmin,xmax,ymin,ymax,δx,δy,cosphi,phi,s02,xy,oneZScore,a,id},
id=0;
For[ix=1,ix≤d-1,ix++,
For[iy=ix+1,iy≤d,iy++,
Print["Plotting plane x"<>ToString[ix]<>"x"<>ToString[iy]];
id++;
x1=plist/.LF[a__]:>{a}[[ix]]; (* x values for p1 *)
y1=plist/.LF[a__]:>{a}[[iy]];  (* y values for p1 *)
(* Average x1 and average y1 *)
(*{x0,y0}={Mean[x1],Mean[y1]} ;*)
{x0,y0}={mean[[ix]],mean[[iy]]};
Print["The center of the ellipse is at {"<>ToString[x0]<>","<>ToString[y0]<>"}"];
{dx, dy} = {x1-x0, y1 - y0}; (* Deviations from the average values *);
{xmin,xmax}={Min[#],Max[#]}&@ x1;
{ymin,ymax}={Min[#],Max[#]}&@ y1;
{δx, δy} ={((xmax-xmin)/2),((ymax-ymin)/2)} ;
(* -> correlation ellipse *)
cosphi=Mean[dx dy]/δx/δy;   
phi=ArcCos[cosphi];
{z1x,z1y} = {1* δx, 1* δy};
(* a correlation ellipse takes the form 0 =1- a(x^2) - b x y - c(y^2) - dx - ey  *)
s02 =((Sin[phi])^2) ;(* constant = 1/f *)
aCM[[id]]= {{(1/(((δx)^2) s02)), -(cosphi/(δx δy s02))},{-(cosphi/(δx δy s02)),(1/(((δy)^2) s02))}}; (* multipling by 1/f*)

(* -> Convex hull ellipse *)
Print["-----------------------"];
xy = plist/.LF[a__]:>{{a}[[ix]],{a}[[iy]]}; (* x and y values *)
oneZScore= Select[xy,(((#[[1]]>(x0-z1x)) && (#[[1]]<(x0+z1x))) && ((#[[2]]>(y0-z1y)) && (#[[2]]<(y0+z1y)) ))&]; 

(* picking points within one zscore *)
cH=ConvexHullMesh[oneZScore]; (* finding convex hull of points *)
convexHullPoints=Cases[Show[cH],GraphicsComplex[pts_,___]:>pts][[1]]; (* isolating the vertices of the convex hull *)
pointsToFit= convexHullPoints/. {x_, y_} :> {x-x0,y-y0};
{d1,d2}={pointsToFit/.{x_,y_}:>{(x^2),x y, (y^2)},pointsToFit/.{x_,y_}:>{1}};
{s11,s12,s21,s22}={Transpose[d1].d1,Transpose[d1].d2,Transpose[d2].d1,Transpose[d2].d2};
{sv11,sv12,sv21,sv22}={linv.Transpose[u1].s11.v1,linv.Transpose[u1].s12,s21.v1,s22};

If[ Det[sv22]==0,
Print["Error: Singular sv22"],
ES=Sort[Eigensystem[sv11-sv12.Inverse[sv22].sv21]//Transpose,#1[[1]]>#2[[1]]&];
NumPositive=Select[Transpose[ES][[1]],#>0&]//Length;
If[NumPositive≠1,
Print["Error: number of positive eigenvalues is ",NumPositive],
av1=Transpose[{ES[[1]][[2]]}];
av2=-Inverse[sv22].sv21.av1;
(* Preliminary estimate for a1 and a2 (with a wrong normalization factor) *)
{a1,a2}={v1.av1,av2};
{X0, Y0} = {0,0};
NormFac = 1;
](*If[NumPositive≠1 *)
];(*If[Det[sv22]==0...*)
a=Transpose[{Join[Transpose[a1][[1]],{0,0,a2[[1]][[1]]}]}]//Flatten;
constant2 = -(((a[[6]]))^(-1)); (* dividing the coefficients by the constant f*)
aHull[[id]] ={{constant2 a[[1]], (1/2)constant2 a[[2]]},{(1/2)constant2 a[[2]],constant2 a[[3]]}}; 

(* Convex hull ->  *)


parametricPlot1 = ParametricPlot[{x0+ δx Cos[u],y0+ δy Cos[u+phi]},{u,0,2π},PlotStyle->RGBColor[1,0,0]];
parametricPlotA=ContourPlot[1-((({{x-x0, y-y0}})).aHull[[id]].(({{x-x0}, {y-y0}})))[[1]][[1]]==0,{x,xmin,xmax},{y,ymin,ymax},ContourStyle->{RGBColor[0.0,0,0],Dashing[{0.06,0.02}],Thickness[0.006]},PlotPoints->100];
parametricPlotB=ContourPlot[1+((({{x-x0, y-y0}})).aCM[[id]].(({{x-x0}, {y-y0}})))[[1]][[1]]==0,{x,xmin,xmax},{y,ymin,ymax},ContourStyle->{RGBColor[0.0,0.5,0.5],Dashing[{0.002,0.02}],Thickness[0.01]},PlotPoints->100];
pp2=ListPlot[plist/.LF[a__]:>{{a}[[ix]],{a}[[iy]]}, PlotMarkers -> {"●",4},PlotStyle->Black,AspectRatio->1,ImageSize->{500,Automatic},Frame->True,AxesOrigin->{0,0},Axes->True,FrameLabel->{"x"<>ToString[ix],"x"<>ToString[iy]}]; 
Print[Show[pp2,parametricPlot1,parametricPlotA,(* parametricPlotB, *)
Graphics[{
RGBColor[0,0,1],Text["★",{x0,y0},{0,0}],
Dashing[{0.02,0.02}],
Line[{{x0-δx,y0-δy},{x0+δx,y0-δy},{x0+δx,y0+δy},{x0-δx,y0+δy},{x0-δx,y0-δy}}],
}], PlotRange->{{x0-1.2δx,x0+1.2δx},{y0-1.2δy,y0+1.2δy}} ,AxesOrigin->{x0,y0}]];
](* For iy -> *)
](* For ix -> *)
](*Module -> *)

Visualize the Original Ellipsoid and the Fitted Ellipses

Here, the ellipses fitted in the previous cell are visualized together with the original 3D ellipsoid

p12=ContourPlot3D[1-((({{x1, x2}})).aHull[[1]].(({{x1}, {x2}})))[[1]][[1]]==0,{x1,-3,3},{x2,-3,3},{x3,-3-(10^(-8)),3+(10^(-8))},Contours->10,ContourStyle->Directive[Opacity[0.2],Orange],Mesh->None];
p13=ContourPlot3D[1-((({{x1, x3}})).aHull[[2]].(({{x1}, {x3}})))[[1]][[1]]==0,{x1,-3,3},{x2,-3-(10^(-8)),3+(10^(-8))},{x3,-3,3},Contours->10,ContourStyle->Directive[Opacity[0.2],Magenta],Mesh->None];
p23=ContourPlot3D[1-((({{x2, x3}})).aHull[[3]].(({{x2}, {x3}})))[[1]][[1]]==0,{x1,-3-(10^(-8)),3+(10^(-8))},{x2,-3,3},{x3,-3,3},Contours->10,ContourStyle->Directive[Opacity[0.2],Cyan],Mesh->None];

Show[p2,p1,p12,p13,p23]

Examine the Input Ellipsoid and the Matrices of the Individual Projections

The original input ellipsoid and the three matrices determining the projections’ ellipses are as follows:

Print["The matrix determining the original ellipsoid is:"];
Print[aOriginal];
Print["The matrix determining the x1 vs x2 projection's ellipse is:"];
Print[aHull[[1]]];
Print["The matrix determining the x1 vs x3 projection's ellipse is:"];
Print[aHull[[2]]];
Print["The matrix determining the x2 vs x3 projection's ellipse is:"];
Print[aHull[[3]]];

If we take the inverses of these matrices, the relationship between the coefficients of the projections’ ellipses and the coefficients of the full ellipsoid becomes clear:

Print["The inverse of the matrix determining the original ellipsoid is:"];
Print[Inverse[aOriginal]];
Print["The inverse of the matrix determining the x1 vs x2 projection's ellipse is:"];
Print[Inverse[aHull[[1]]]];
Print["The inverse of the matrix determining the x1 vs x3 projection's ellipse is:"];
Print[Inverse[aHull[[2]]]];
Print["The inverse of the matrix determining the x2 vs x3 projection's ellipse is:"];
Print[Inverse[aHull[[3]]]];

Reconstruct and Evaluate the 3D Ellipsoid

Reconstruct the matrix aOutFinal of the 3 D ellipsoid from the matrix aIn of three 2 D projections of the ellipsoid. If the projections belong to the same ellipsoid, their diagonal coefficients can be arranged into a “consistency coefficient” acon computed below that must be close to zero.

i=1;
{aOutCM,aOutHull}={ConstantArray[0,{d,d}],ConstantArray[0,{d,d}]};
{adiagCM,adiagHull}={ConstantArray[{},d],ConstantArray[{},d]};
{daOutCM,daOutHull}={ConstantArray[0,d],ConstantArray[0,d]};
{aconCM,aconHull}={1,1};
For[ix=1,ix≤d-1,ix++,
For[iy=ix+1,iy≤d,iy++,
{aInHull,aInCM}={Inverse[aHull[[i]]],Inverse[aCM[[i]]]};
aOutCM[[ix,iy]]=aOutCM[[iy,ix]]=aInCM[[1,2]];
aOutHull[[ix,iy]]=aOutHull[[iy,ix]]=aInHull[[1,2]];
adiagCM[[ix]]=Append[adiagCM[[ix]],aInCM[[1,1]]];
adiagCM[[iy]]=Append[adiagCM[[iy]],aInCM[[2,2]]];
adiagHull[[ix]]=Append[adiagHull[[ix]],aInHull[[1,1]]];
adiagHull[[iy]]=Append[adiagHull[[iy]],aInHull[[2,2]]];
i++
]
]

For[ix=1,ix≤d,ix++,
{aOutCM[[ix,ix]],aOutHull[[ix,ix]]}={Mean[adiagCM[[ix]]],Mean[adiagHull[[ix]]]};
{daOutCM[[ix]],daOutHull[[ix]]}={StandardDeviation[adiagCM[[ix]]],StandardDeviation[adiagHull[[ix]]]};
aconCM*=daOutCM[[ix]]/aOutCM[[ix]][[ix]];
aconHull*=daOutHull[[ix]]/aOutHull[[ix]][[ix]];
]
{aconCM,aconHull}={(((Abs[aconCM]))^(1/d)),(((Abs[aconHull]))^(1/d))}//Chop;
{aFinalCM,aFinalHull}={Inverse[aOutCM],Inverse[aOutHull]};
Print["The input ellipsoid matrix is (Subscript[A, original])=",aOriginal];
Print["The reconstructed matrices are (Subscript[A, Hull])=",aFinalHull," and  (Subscript[A, cov . mat . ])=",aFinalCM];
Print["The consistency coefficients (average differences in diagonal elements) are (Subscript[acon, Hull])=",aconHull," and (Subscript[acon, CM])=",aconCM];