Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

From this question I successfully made an elliptical fit for my data. However, when I try and collect the datapoints within the ellipse via this question the coordinates of my ellipse correspond to an ellipse of a rasterized image of my listplot data and not the original. How can I transform or scale the ellipse that fits the rasterized image to the original listplot.

Ok, with all that said heres an example.

data = RandomReal[NormalDistribution[], {100000, 2}]
p = ListPlot[data, ImageSize -> 4000];
f = FillingTransform@ColorNegate@Binarize@p // DeleteSmallComponents
{c, s, t} = 1 /. ComponentMeasurements[f, {"Centroid", "SemiAxes", "Orientation"}]
Show[Rasterize[p], Graphics[{Red, Rotate[Circle[c, s],t]}]]

where I get a nice image:

example data

however, the coordinates c, s, t of the ellipse are pixle coordiantes corresponding to the rasterized image rather than data coordinates.

So when I need the ellipses parameter to do any calculations I get bunk results.

The image processing approach would be the best as I am filtering out the most dense cluster of data.

Thanks so much.

share|improve this question

migrated from stackoverflow.com Aug 5 '12 at 13:26

This question came from our site for professional and enthusiast programmers.

1  
related stackoverflow.com/questions/11783394/… –  belisarius Aug 5 '12 at 0:05
    
I added another solution to your question, at your referenced question here: mathematica.stackexchange.com/questions/4519/… –  Chris Degnen Aug 6 '12 at 10:06
add comment

2 Answers

up vote 6 down vote accepted
(*Generate Data and fit*)
data1 = RandomReal[NormalDistribution[10, 1], {10^4}];(*test data*)
data2 = RandomReal[NormalDistribution[20, 5], {10^4}];(*test data*)
data  = Transpose@{data1, data2};
r     = RotationTransform[Pi/8];
data3 = r /@ data;

(*we need to specify PlotRange due to a kown bug in AbsoluteOptions[] *)
prange = {Min@#, Max@#} & /@ {First@#, Last@#} &@Transpose@data3;

p = ListPlot[data3, Axes -> None, PlotRange -> prange];
f = FillingTransform@ColorNegate@Binarize@p // DeleteSmallComponents;
{co, so, to} = 1 /. ComponentMeasurements[f, {"Centroid", "SemiAxes", "Orientation"}];

(*Transform Image to Graphic coordinates*)
c = Rescale[co[[#]], {1, ImageDimensions[f][[#]]}, prange[[#]]] & /@ {1, 2};
s = Rescale[so[[#]], {0, Norm@ImageDimensions@f}, {0, Norm@(Differences/@ prange)}] & /@ {1, 2};
t = -ArcTan@Rescale[Tan@to, {0, 1/Divide @@ ImageDimensions[f]}, 
                            {0, 1/First@(Divide @@ Differences /@ prange)}];
(* Replot graphic*)

{s1, s2} = s;
{cx, cy} = c;
f0 = Sqrt[s1 s1 - s2 s2];
f1 = {cx + f0 Cos[t], cy - f0 Sin[t]};
f2 = {cx - f0 Cos[t], cy + f0 Sin[t]};
r = 2 Sqrt[f0 f0 + s2 s2];

sd = Select[data3, EuclideanDistance[#, f1] + EuclideanDistance[#, f2] < r &];
Show[p, Graphics[{Red, PointSize[Large], Point@sd}], Axes -> True]

Mathematica graphics

share|improve this answer
1  
Thanks again bel. If there is anything more that I can do to thank you let me know. –  tquarton Aug 5 '12 at 16:52
1  
add comment

Here is repost of the code mention in my comment above, but using belisarius' data set since the data generator I used doesn't work in Mathematica 7. (No harm posting the code twice I guess).

It uses the EllipsoidQuantile function (which I think is pretty neat), with standard deviation set to 2, which you can change as you please.

(* belisarius' data set up code *)
data1 = RandomReal[NormalDistribution[10, 1], {1000}]; (* test data *)
data2 = RandomReal[NormalDistribution[20, 5], {1000}]; (* test data *)
data = Transpose@{data1, data2};
r = RotationTransform[Pi/8];
data3 = r /@ data;
prange = {Min@#, Max@#} & /@ {First@#, Last@#} &@Transpose@data3;

{{xmin, xmax}, {ymin, ymax}} = prange;
(* For values within two standard deviations, (approx 95.45% of values) *)
sd = 2;
cl = 2*(CDF[NormalDistribution[0, 1], sd] - 0.5);
Needs["MultivariateStatistics`"];
e = EllipsoidQuantile[data3, cl];
ctr = e[[1]];
{r1, r2} = e[[2]];
inc = ArcTan[e[[3, 1, 2]]/e[[3, 1, 1]]]*180/Pi;
Print["Ellipse center = " <> ToString@ctr];
Print["Ellipse radii (r1, r2) = " <> ToString@{r1, r2}]; Print[
 StringJoin["Ellipse inclination = ", ToString@inc, " degrees"]];

(* Find the foci of the ellipse *)
f = Sqrt[r1^2 - r2^2];
dx = f*Cos[inc Degree];
dy = f*Sin[inc Degree];
f1 = ctr - {dx, dy};
f2 = ctr + {dx, dy};

edge = ctr + r1*e[[3, 1]];
rlim = EuclideanDistance[edge, f1] + EuclideanDistance[edge, f2];
(* more of belisarius' code *)
inside[{x_, y_}, {f1_, f2_}] := 
  Sum[EuclideanDistance[{x, y}, i], {i, {f1, f2}}];
sd = Select[data3, inside[#, {f1, f2}] < rlim &];

Show[
 RegionPlot[
  inside[{x, y}, {f1, f2}] < rlim, {x, xmin, xmax}, {y, ymin, ymax}],
 ListPlot[data3],
 Graphics[{Green, Point@sd}],
 Graphics@e,
 Graphics[{Black, Thick, Dashing[0.05],
   Rotate[Circle[ctr, {r1, r2}], inc Degree]}],
 Graphics[{Red, Line[{ctr + r1*e[[3, 1]], ctr, ctr + r2*e[[3, 2]]}]}],
 Graphics[{PointSize[Large], Point[{f1, f2}]}],
 PlotRange -> {{xmin, xmax}, {ymin, ymax}}, 
 AspectRatio -> (ymax - ymin)/(xmax - xmin)]

enter image description here

share|improve this answer
    
Thanks Chris! That's a nice approach as well! Thank you for your time to give me more options. –  tquarton Aug 6 '12 at 20:39
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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