001/*
002 * Copyright 2010-2015 Institut Pasteur.
003 * 
004 * This file is part of Icy.
005 * 
006 * Icy is free software: you can redistribute it and/or modify
007 * it under the terms of the GNU General Public License as published by
008 * the Free Software Foundation, either version 3 of the License, or
009 * (at your option) any later version.
010 * 
011 * Icy is distributed in the hope that it will be useful,
012 * but WITHOUT ANY WARRANTY; without even the implied warranty of
013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
014 * GNU General Public License for more details.
015 * 
016 * You should have received a copy of the GNU General Public License
017 * along with Icy. If not, see <http://www.gnu.org/licenses/>.
018 */
019package plugins.kernel.roi.roi2d;
020
021import icy.math.ArrayMath;
022import icy.resource.ResourceUtil;
023import icy.sequence.Sequence;
024import icy.type.point.Point5D;
025
026import java.awt.geom.Ellipse2D;
027import java.awt.geom.Point2D;
028import java.awt.geom.Rectangle2D;
029import java.util.Collection;
030
031/**
032 * @author Stephane
033 */
034public class ROI2DEllipse extends ROI2DRectShape
035{
036    /**
037     * @deprecated
038     */
039    @Deprecated
040    public ROI2DEllipse(Point2D topLeft, Point2D bottomRight, boolean cm)
041    {
042        this(topLeft, bottomRight);
043    }
044
045    public ROI2DEllipse(Point2D topLeft, Point2D bottomRight)
046    {
047        super(new Ellipse2D.Double(), topLeft, bottomRight);
048
049        // set icon (default name is defined by getDefaultName()) 
050        setIcon(ResourceUtil.ICON_ROI_OVAL);
051    }
052
053    /**
054     * Create a ROI ellipse from its rectangular bounds.
055     */
056    public ROI2DEllipse(double xmin, double ymin, double xmax, double ymax)
057    {
058        this(new Point2D.Double(xmin, ymin), new Point2D.Double(xmax, ymax));
059    }
060
061    /**
062     * @deprecated
063     */
064    @Deprecated
065    public ROI2DEllipse(Rectangle2D rectangle, boolean cm)
066    {
067        this(rectangle);
068    }
069
070    public ROI2DEllipse(Rectangle2D rectangle)
071    {
072        this(rectangle.getMinX(), rectangle.getMinY(), rectangle.getMaxX(), rectangle.getMaxY());
073    }
074
075    public ROI2DEllipse(Ellipse2D ellipse)
076    {
077        this(ellipse.getBounds2D());
078    }
079
080    /**
081     * @deprecated
082     */
083    @Deprecated
084    public ROI2DEllipse(Point2D pt, boolean cm)
085    {
086        this(pt);
087    }
088
089    public ROI2DEllipse(Point2D pt)
090    {
091        this(new Point2D.Double(pt.getX(), pt.getY()), pt);
092    }
093
094    /**
095     * Generic constructor for interactive mode
096     */
097    public ROI2DEllipse(Point5D pt)
098    {
099        this(pt.toPoint2D());
100    }
101
102    public ROI2DEllipse()
103    {
104        this(new Point2D.Double(), new Point2D.Double());
105    }
106    
107    @Override
108    public String getDefaultName()
109    {
110        return "Ellipse2D";
111    }
112
113    public Ellipse2D getEllipse()
114    {
115        return (Ellipse2D) shape;
116    }
117
118    public void setEllipse(Ellipse2D ellipse)
119    {
120        setBounds2D(ellipse.getBounds2D());
121    }
122
123    @Override
124    public double getLength(Sequence sequence) throws UnsupportedOperationException
125    {
126        final Ellipse2D ellipse = getEllipse();
127        return computeEllipsePerimeter(ellipse.getWidth() * 0.5d * sequence.getPixelSizeX(), ellipse.getHeight() * 0.5d
128                * sequence.getPixelSizeY());
129    }
130
131    @Override
132    public double computeNumberOfContourPoints()
133    {
134        final Ellipse2D ellipse = getEllipse();
135        return computeEllipsePerimeter(ellipse.getWidth() * 0.5d, ellipse.getHeight() * 0.5d);
136    }
137
138    /**
139     * Calculating the perimeter of an ellipse is non-trivial. Here we follow the approximation
140     * proposed in:<br/>
141     * Ramanujan, S., "Modular Equations and Approximations to Pi," Quart. J. Pure. Appl. Math.,
142     * v.45 (1913-1914), 350-372
143     * 
144     * @since Icy 1.5.3.2
145     */
146    public static double computeEllipsePerimeter(double w, double h)
147    {
148        double result = (w - h) / (w + h);
149        result *= result;
150
151        return Math.PI * (w + h) * (1 + (result / 4) + ((result * result) / 64) + ((result * result * result) / 256));
152
153    }
154
155    @Override
156    public double computeNumberOfPoints()
157    {
158        final Ellipse2D ellipse = getEllipse();
159        return Math.PI * ellipse.getWidth() * ellipse.getHeight() / 4d;
160    }
161
162    /**
163     * Adjust the ROI to fit the specified list of coordinates with a circle
164     * 
165     * @param points
166     *        the list of points to fit
167     */
168    public void setToFitCircle(Collection<? extends Point2D> points)
169    {
170        int nbPoints = points.size();
171
172        double[] xCoords = new double[nbPoints];
173        double[] yCoords = new double[nbPoints];
174
175        for (Point2D point : points)
176        {
177            nbPoints--;
178            xCoords[nbPoints] = point.getX();
179            yCoords[nbPoints] = point.getY();
180        }
181
182        setToFitCircle(xCoords, yCoords);
183    }
184
185    /**
186     * Circle fit by Taubin. <br/>
187     * Reference: G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar Space Curves
188     * Defined By Implicit
189     * Equations, With Applications To Edge And Range Image Segmentation", IEEE Trans. PAMI, Vol.
190     * 13, pages 1115-1138,
191     * (1991).<br/>
192     * This method is a port to Icy from the original <a
193     * href="http://www.mathworks.com/matlabcentral/fileexchange/22678">Matlab code from Nikolai
194     * Chernov (2009)</a>
195     * 
196     * @param xCoords
197     *        the X coordinates of the points to fit
198     * @param yCoords
199     *        the Y coordinates of the points to fit
200     */
201    private void setToFitCircle(double[] xCoords, double[] yCoords)
202    {
203        int n = xCoords.length;
204
205        if (n != yCoords.length)
206            throw new IllegalArgumentException("Coordinate arrays must have the same size");
207
208        Point2D centroid = new Point2D.Double(ArrayMath.mean(xCoords), ArrayMath.mean(yCoords));
209
210        // temporary buffer to save memory
211        double[] buffer = new double[n];
212
213        double[] X = ArrayMath.subtract(xCoords, centroid.getX());
214        double[] Y = ArrayMath.subtract(yCoords, centroid.getY());
215        double[] XX = ArrayMath.multiply(X, X);
216        double[] YY = ArrayMath.multiply(Y, Y);
217        double[] Z = ArrayMath.add(XX, YY);
218
219        // compute normalized moments
220        double Mxx = ArrayMath.sum(XX) / n;
221        double Mxy = ArrayMath.sum(ArrayMath.multiply(X, Y, buffer)) / n;
222        double Myy = ArrayMath.sum(YY) / n;
223        double Mxz = ArrayMath.sum(ArrayMath.multiply(X, Z, buffer)) / n;
224        double Myz = ArrayMath.sum(ArrayMath.multiply(Y, Z, buffer)) / n;
225        double Mzz = ArrayMath.sum(ArrayMath.multiply(Z, Z, buffer)) / n;
226
227        // computing the coefficients of the characteristic polynomial
228
229        double Mz = Mxx + Myy;
230        double Cov_xy = Mxx * Myy - Mxy * Mxy;
231        double A3 = 4 * Mz;
232        double A2 = -3 * Mz * Mz - Mzz;
233        double A1 = Mzz * Mz + 4 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz - Mz * Mz * Mz;
234        double A0 = Mxz * Mxz * Myy + Myz * Myz * Mxx - Mzz * Cov_xy - 2 * Mxz * Myz * Mxy + Mz * Mz * Cov_xy;
235        double A22 = A2 + A2;
236        double A33 = A3 + A3 + A3;
237
238        double xold, xnew = 0;
239        double yold, ynew = 1e+20;
240        double epsilon = 1e-12;
241        int IterMax = 20;
242
243        // Newton's method starting at x=0
244
245        for (int iter = 0; iter < IterMax; iter++)
246        {
247            yold = ynew;
248            ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
249            if (Math.abs(ynew) > Math.abs(yold))
250            {
251                System.err.println("Circle fitting error: Newton-Taubin goes wrong direction: |ynew| > |yold|");
252                xnew = 0;
253                break;
254            }
255
256            double Dy = A1 + xnew * (A22 + xnew * A33);
257            xold = xnew;
258            xnew = xold - ynew / Dy;
259
260            if (Math.abs((xnew - xold) / xnew) < epsilon)
261                break;
262
263            if (iter >= IterMax)
264            {
265                System.err.println("Circle fitting error: Newton-Taubin will not converge");
266                xnew = 0;
267            }
268
269            if (xnew < 0)
270            {
271                System.out.println("Newton-Taubin negative root: x=" + xnew);
272                xnew = 0;
273            }
274        }
275
276        // computing the circle parameters
277
278        double DET = xnew * xnew - xnew * Mz + Cov_xy;
279        double xCenter = (Mxz * (Myy - xnew) - Myz * Mxy) / DET / 2;
280        double yCenter = (Myz * (Mxx - xnew) - Mxz * Mxy) / DET / 2;
281        double radius = Math.sqrt(xCenter * xCenter + yCenter * yCenter + Mz);
282        xCenter += centroid.getX();
283        yCenter += centroid.getY();
284
285        setEllipse(new Ellipse2D.Double(xCenter - radius, yCenter - radius, 2 * radius, 2 * radius));
286    }
287}