4
$\begingroup$

I want to fit (the best/closest possible) parallelogram to a set of XY data points and obtain the co-ordinates of the vertices.

  1. The X and Y data are two independent sinusoidal type of waves, which on plotting against each other give a parallelogram-like Lissajous figure. Example data can be found in this link -- Example Data

  2. The parallelogram (almost) always has the 4 vertices in four different quadrants.

At the moment I am manually fitting a parallelogram to extract the co-ordinates. I have 60 such data sets. Any kind of help is greatly appreciated. (I am new to Mathematica).

EDIT: The result of

BoundingRegion[myData, "MinOrientedRectangle"]

in yellow in the linked image is not what I am trying to obtain in this case. The red lines in the image is what I am trying to achieve, and then extract the cyan coloured co-ordinates.

enter image description here

$\endgroup$
  • 2
    $\begingroup$ BoundingRegion[pts, "MinOrientedRectangle"] could be a useful starting point. $\endgroup$ – J. M. will be back soon Oct 5 '16 at 16:26
  • $\begingroup$ @jackryan - your image didn't come through, it seems like an empty jpeg $\endgroup$ – Jason B. Oct 5 '16 at 17:57
  • $\begingroup$ @JasonB Apologies; I just added a link to the image. $\endgroup$ – jackryan Oct 5 '16 at 18:02
10
$\begingroup$

Okay, the way we want to solve that is easy. We'll look at the distance to the center. This gives us a data function we can fit with a trigonometric function.

So lets start:

lData=Sqrt[Plus@@(#^2)]&/@data; (*compute distance data*)
maxAmplitude=Max[lData]; (*find maximum*)
lpData=Table[{i/Length[lData]*10Pi,lData[[i]]/maxAmplitude*2},{i,1,Length[lData]}]; (*rescale data to [0,10Pi] and [0,2]*)

fit=NonlinearModelFit[lpData,Cos[a*x+b]+1,{a,b},x]; (*fit with a simple cosine function*)

Show[ListPlot[lpData,PlotStyle->Black],Plot[fit["BestFit"],{x,0,10Pi},PlotStyle->Red],Graphics[{Blue,PointSize[Large],Point[{#,2}]&/@Table[(-b+2Pi*i)/a/.fit["BestFitParameters"],{i,1,5}]}]] (*display result*)

enter image description here

The next task is pretty starightforward. We search for the Max,Min-Points. We know that:

$$x_{max}=\frac{-b+2\pi\cdot n}{a}$$ $$x_{min}=\frac{-b+2\pi\cdot n+\pi}{a}$$

With $n\in\mathbb{N}$. So we can use that to rescale to our data-range. We'll get multiple points so we take the average from them:

mulCornerData=data[[Round[Length[data]*(Table[(-b+Pi*i)/a/.fit["BestFitParameters"],{i,1,10}]/(10Pi))]]];
cornerData={Mean[mulCornerData[[#]]&/@{1,5,9}],Mean[mulCornerData[[#]]&/@{2,6,10}],Mean[mulCornerData[[#]]&/@{3,7}],Mean[mulCornerData[[#]]&/@{4,8}]};
Show[ListPlot[data,PlotStyle->Black],ListLinePlot[cornerData[[Mod[#,4]+1]]&/@Range[0,4],PlotStyle->Red]]

enter image description here

$\endgroup$
  • $\begingroup$ You made it so simple; Very big thanks! Is there way to add a constraint to strictly get all co-ordinates in different quadrants? $\endgroup$ – jackryan Oct 5 '16 at 19:59
7
$\begingroup$

There is already an answer, but I think with less assumptions we can create a better fit.

After loading the data I transformed them in polar coordinates (angle, radius) and chopped off the first pair.

dataTransf = {ArcTan[First@#, Last@#], Norm[#]} & /@ data[[2 ;;]];

In polar coordinates you can represent a line with the $Csc(\phi-\alpha)$ where $\alpha$ is the angle of the line. With this knowledge we build a function of 5 (as there are some data points prior to the first peak). We got 2 pair of lines, were each has the same angle, but different signum.

fitfunc = Piecewise[{
{-c2 Csc[ϕ - ϕ2], ϕ < a - π},
{c1 Csc[ϕ - ϕ1], a - π <= ϕ <= -π + b},
{c2 Csc[ϕ - ϕ2], -π + b < ϕ < a},
{-c1 Csc[ϕ - ϕ1], a <= ϕ <= b},
{-c2 Csc[ϕ - ϕ2], a < ϕ < π + a}}]

To become an rectangle the lines have to be continuous. So we use continuity conditions at a and b:

-c2 Csc[a - ϕ2] == c1 Csc[a - ϕ1]
c1 Csc[b - ϕ1] == c2 Csc[b - ϕ2]

And solve those for a and b to get: Edit: Of coure Mathematica can do this

Solve[-c2 Csc[a - ϕ2] == c1 Csc[a - ϕ1], a]
Solve[c1 Csc[b - ϕ1] == c2 Csc[b - ϕ2], b]

You get 4 possible solutions each and I have to admitt: I dont know why but it works fine with the 2nd solution for a and the 3rd for b.

{a -> ArcCos[-((c2 Cos[ϕ1] + c1 Cos[ϕ2])/Sqrt[
c2^2 Cos[ϕ1]^2 + 2 c1 c2 Cos[ϕ1] Cos[ϕ2] + 
 c1^2 Cos[ϕ2]^2 + c2^2 Sin[ϕ1]^2 + 
 2 c1 c2 Sin[ϕ1] Sin[ϕ2] + c1^2 Sin[ϕ2]^2])],
b->ArcCos[(c2 Cos[ϕ1] - c1 Cos[ϕ2])/Sqrt[
c2^2 Cos[ϕ1]^2 - 2 c1 c2 Cos[ϕ1] Cos[ϕ2] + 
c1^2 Cos[ϕ2]^2 + c2^2 Sin[ϕ1]^2 - 
2 c1 c2 Sin[ϕ1] Sin[ϕ2] + c1^2 Sin[ϕ2]^2]]}

By replacing those in the equation we reduce it to 4 parameters $\phi 1$, $\phi 2$, $c1$, $c2$. Now after all the prework we can actually fit.

fitmod = NonlinearModelFit[dataTransf, 
fitfunc, {{ϕ1, -1}, ϕ2, c1, {c2, 100}}, ϕ];

Here we have to take care of the intial values. (-1 for \ϕ1 and 100 for c2 do the job.)

Now we plot it nicely

Show[ListPlot[data, PlotRange -> {{-3600, 3600}, {-800, 800}}],
PolarPlot[fitmod[ϕ], {ϕ, -π, π}, 
ColorFunction -> Function[{x, y, z}, Hue[0, 1, 1]]], 
ImageSize -> Large]

Show[ListPlot[dataTransf, PlotRange -> {-1000, 3600}], 
Plot[fitmod[ϕ], {ϕ, -π, π}, 
PlotRange -> {-1000, 3300}, 
ColorFunction -> Function[{x, y}, Hue[0.0, 1, 1]]], 
ImageSize -> Large]

and get:

DataFit

enter image description here

Edit2:

fitmod["BestFitParameters"]

gives the parameters. The parameter $a,b,a+\pi,b+\pi$ describe the angles of the vertices and evaluating fitfunc there gives the radius. Then we can use:

$$x=r \sin(\phi)$$ $$y=r\cos(\phi)$$

If there are still any questions concerning the math or the code let me know.

$\endgroup$
  • $\begingroup$ pretty cool. Well done. $\endgroup$ – Julien Kluge Oct 6 '16 at 16:46
  • 1
    $\begingroup$ @meneken17 Very elegant approach (took 3-4 hour to understand the maths). [Phi]1 = 0 worked on Mathematica v11. 1. Could you please put the code for solution to a and b. 2. Also, how can I extract the co-ordinates of the four vertices? Thank you. $\endgroup$ – jackryan Oct 6 '16 at 19:59

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.