001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 * image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause).
006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009
010package imagingbook.common.geometry.fd;
011
012import imagingbook.common.geometry.basic.Pnt2d;
013import imagingbook.common.math.Complex;
014import org.apache.commons.math3.analysis.UnivariateFunction;
015import org.apache.commons.math3.optim.MaxEval;
016import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
017import org.apache.commons.math3.optim.univariate.BrentOptimizer;
018import org.apache.commons.math3.optim.univariate.SearchInterval;
019import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
020import org.apache.commons.math3.optim.univariate.UnivariateOptimizer;
021import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
022
023import java.awt.Shape;
024import java.awt.geom.AffineTransform;
025import java.awt.geom.Ellipse2D;
026import java.awt.geom.Path2D;
027import java.util.Locale;
028
029import static imagingbook.common.math.Arithmetic.mod;
030import static imagingbook.common.math.Arithmetic.sqr;
031
032/**
033 * <p>
034 * This class represents elliptic Fourier descriptors. See Ch. 26 of [1] for additional details including invariance
035 * calculations.
036 * </p>
037 * <p>
038 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction Using Java</em>, 2nd ed,
039 * Springer (2016).
040 * </p>
041 *
042 * @author WB
043 * @version 2022/10/24
044 * @see FourierDescriptorUniform
045 * @see FourierDescriptorTrigonometric
046 */
047public class FourierDescriptor {
048
049        public static final int MinReconstructionSamples = 50;
050
051        private final int mp;           // number of Fourier pairs
052        private final Complex[] G;      // complex-valued DFT spectrum (of length 2 * mp + 1)
053        private final double scale;     // scale factor to apply for reconstructing original (sample) size
054
055        /**
056         * Constructor using the default scale (1.0).
057         *
058         * @param G a complex-valued DFT spectrum
059         */
060        public FourierDescriptor(Complex[] G) {
061                this(G, 1.0);
062        }
063
064        /**
065         * Constructor using a specific scale.
066         *
067         * @param G a complex-valued DFT spectrum
068         * @param scale the reconstruction scale
069         */
070        public FourierDescriptor(Complex[] G, double scale) {
071                this.mp = (G.length - 1) / 2;
072                this.G = truncate(G, this.mp);  // always make G odd-sized
073                this.scale = scale;
074        }
075
076        /**
077         * Constructor for cloning Fourier descriptors.
078         *
079         * @param fd an existing instance
080         */
081        public FourierDescriptor(FourierDescriptor fd) {
082                this(fd.G.clone(), fd.scale);
083        }
084
085        // ----------------------------------------------------------------
086
087        /**
088         * <p>
089         * Truncate the given DFT spectrum to the specified number of coefficient pairs. Truncation removes the
090         * highest-frequency coefficients. The resulting spectrum is always odd-sized. If the number of coefficient pairs is
091         * zero, only coefficient 0 remains, i.e., the new spectrum has length 1. An exception is thrown if the original
092         * spectrum has fewer coefficient pairs than needed. For example, for an even-sized spectrum with 10 coefficients
093         * G[m] = a,b,...,j,
094         * </p>
095         * <pre>
096         * m    = 0 1 2 3 4 5 6 7 8 9
097         * G[m] = a b c d e f g h i j </pre>
098         * <p>
099         * the truncated spectrum for mp = 3 is
100         * </p>
101         * <pre>
102         * m'    = 0 1 2 3 4 5 6
103         * G[m'] = a b c d h i j </pre>
104         * <p>
105         * I.e., the highest frequency coefficients (e, f, g) of the original spectrum are removed.
106         *
107         * @param G the original DFT spectrum (with length greater than 2 * mp + 1)
108         * @param mp the number of remaining coefficient pairs.
109         * @return the truncated spectrum (always odd-sized)
110         */
111        static Complex[] truncate(Complex[] G, int mp) {
112                if (mp < 0) {
113                        throw new IllegalArgumentException("number of coefficient pairs must be >= 0 but is " + mp);
114                }
115                int M = G.length;
116                int Mnew = 2 * mp + 1;
117                if (Mnew > M) {
118                        throw new IllegalArgumentException("spectrum has fewer coefficient pairs than needed");
119                }
120                Complex[] Gnew = new Complex[Mnew];
121                // fill the new spectrum:
122                for (int m = 0; m < Mnew; m++) {
123                        if (m <= Mnew / 2)
124                                Gnew[m] = G[m];
125                        else
126                                Gnew[m] = G[M - Mnew + m];
127                }
128                return Gnew;
129        }
130
131        /**
132         * Truncates the given DFT spectrum to the maximum number of coefficient pairs. The resulting spectrum is always
133         * odd-sized. If the original spectrum is odd-sized, the same spectrum is return, otherwise the single
134         * highest-frequency coefficient is removed.
135         *
136         * @param G the original DFT spectrum
137         * @return the truncated spectrum (always odd-sized)
138         * @see #truncate(Complex[], int)
139         */
140        static Complex[] truncate(Complex[] G) {
141                return truncate(G, (G.length - 1) / 2);
142        }
143
144        // ----------------------------------------------------------------
145
146        /**
147         * Returns the size (length) of this Fourier descriptor's array of DFT coefficients G. The resulting value is always
148         * odd.
149         *
150         * @return the size of the DFT coefficient array (M)
151         */
152        public int size() {
153                return G.length;
154        }
155
156        /**
157         * Returns the scale of this Fourier descriptor, that is, the factor required to scale it to its original (sample)
158         * size. The scale factor is changed in the process of making descriptors scale-invariant.
159         *
160         * @return the scale factor
161         */
162        public double getScale() {
163                return scale;
164        }
165
166        /**
167         * Returns the number of Fourier coefficient pairs for this descriptor. The result is (M-1)/2, M being the size of
168         * the DFT coefficient array G.
169         *
170         * @return the number of Fourier coefficient pairs
171         */
172        public int getMp() {
173                return this.mp;
174        }
175
176        /**
177         * Returns the array G = {G[0], ..., G[M-1]} of complex-valued DFT coefficients.
178         *
179         * @return the array of DFT coefficients
180         */
181        public Complex[] getCoefficients() {
182                return G;
183        }
184
185        /**
186         * <p>
187         * Returns the complex-valued DFT coefficient with the specified frequency index m. The returned coefficient is G[k]
188         * with k = (m mod G.length). Unique coefficients are returned for m = 0, ..., M-1, where M is the size of the DFT
189         * coefficient array, or m = -mp, ..., +mp where mp is the number of coefficient pairs. For example, given a Fourier
190         * descriptor with a 9-element spectrum,
191         * </p>
192         * <pre>
193         * k    = 0 1 2 3 4 5 6 7 8
194         * G[k] = a b c d e f g h i </pre>
195         * <p>
196         * the following values are returned:
197         * </p>
198         * <pre>
199         * getCoefficient(0)  &rarr; G[0] = a
200         * getCoefficient(1)  &rarr; G[1] = b
201         * ...
202         * getCoefficient(4)  &rarr; G[4] = e
203         * getCoefficient(-1) &rarr; G[8] = i
204         * ...
205         * getCoefficient(-4) &rarr; G[5] = f
206         * getCoefficient(-5) &rarr; G[4] = e
207         * ...</pre>
208         *
209         * @param m frequency index (positive or negative)
210         * @return the associated DFT coefficient
211         * @see #getCoefficientIndex(int)
212         */
213        public Complex getCoefficient(int m) {
214                return G[getCoefficientIndex(m)];
215        }
216
217        /**
218         * Returns the Fourier coefficient pair (G[-m], G[+m]) as a {@link Complex} valued array.
219         *
220         * @param m frequency index
221         * @return the DFT coefficient pair
222         */
223        public Complex[] getCoefficientPair(int m) {
224                return new Complex[] {G[getCoefficientIndex(-m)], G[getCoefficientIndex(+m)]};
225        }
226
227        /**
228         * Returns the coefficient array index for the specified frequency index, which may be positive or negative.
229         *
230         * @param m frequency index (positive or negative)
231         * @return the coefficient array index
232         * @see #getCoefficient(int)
233         */
234        public int getCoefficientIndex(int m) {
235                return mod(m, G.length);
236        }
237
238        // ----------------------------------------------------------------
239        // Invariance-related methods
240        // ----------------------------------------------------------------
241
242        /**
243         * <p>
244         * Makes this Fourier descriptor invariant to scale, start-point and rotation. The descriptors center position
245         * (coefficient 0) is preserved. Performs the following normalization steps in sequence:
246         * </p>
247         * <ol>
248         * <li>scale invariance,</li>
249         * <li>start-point invariance,</li>
250         * <li>rotation invariance.</li>
251         * </ol>
252         * <p>
253         * Multiple candidate descriptors are returned, since start-point invariance is
254         * not unique. See Sec. 26.5 (Alg. 26.2) of [1] for additional details. The
255         * original (this) descriptor is not modified.
256         * </p>
257         * <p>
258         * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An
259         * Algorithmic Introduction Using Java</em>, 2nd ed, Springer (2016).
260         * </p>
261         *
262         * @return an array of modified Fourier descriptors
263         * @see #makeScaleInvariant()
264         * @see #makeStartPointInvariant()
265         * @see #makeRotationInvariant()
266         */
267        public FourierDescriptor[] makeInvariant() {
268                // Step 1: make scale invariant
269                FourierDescriptor fdS = this.makeScaleInvariant();
270
271                // Step 2: make start-point invariant (not unique, produces 2 versions)
272                FourierDescriptor[] fdAB = fdS.makeStartPointInvariant();
273
274                // Step 3: make all versions rotation invariant
275                for (int i = 0; i < fdAB.length; i++) {
276                        fdAB[i] = fdAB[i].makeRotationInvariant();
277                }
278
279                return fdAB;
280        }
281        
282        // ----------------------------------------------------------------
283
284        /**
285         * Returns a new scale invariant Fourier descriptor by normalizing the L2 norm of the sub-vector {G[-mp], ...,
286         * G[-1], G[1], ..., G[mp]}. Coefficient G_0 remains unmodified. The new Fourier descriptor carries the associated
287         * scale (see {@link #getScale()}). The original Fourier descriptor is not modified.
288         *
289         * @return a new scale-normalized Fourier descriptor
290         */
291        public FourierDescriptor makeScaleInvariant() {
292                double sum = 0;
293                for (int i = 1; i < G.length; i++) {
294                        sum = sum + G[i].abs2();
295                }
296                double scale = Math.sqrt(sum);  // = L2 norm
297                double s = 1 / scale;
298                
299                Complex[] G2 = new Complex[G.length];
300                G2[0] = G[0];
301                for (int i = 1; i < G2.length; i++) {
302                        G2[i] = G[i].multiply(s);
303                }
304                return new FourierDescriptor(G2, scale);
305        }
306
307        /**
308         * Returns a pair of start-point normalized Fourier descriptors. The original Fourier descriptor is not modified.
309         *
310         * @return a pair of start-point normalized Fourier descriptors
311         */
312        public FourierDescriptor[] makeStartPointInvariant() {
313                double phi = getStartPointPhase();
314
315                Complex[] Ga = shiftStartPointPhase(phi);
316                Complex[] Gb = shiftStartPointPhase(phi + Math.PI);
317                
318                return new FourierDescriptor[] {
319                                new FourierDescriptor(Ga, this.scale), 
320                                new FourierDescriptor(Gb, this.scale)};
321        }
322        
323        private Complex[] shiftStartPointPhase(double phi) {
324                Complex[] Gnew = G.clone();
325                for (int m = -mp; m <= mp; m++) {
326                        if (m != 0) {
327                                int k = getCoefficientIndex(m);
328                                Gnew[k] = Gnew[k].rotate(m * phi);
329                        }
330                }
331                return Gnew;
332        }
333
334        /**
335         * Calculates the 'canonical' start point. This version uses (a) a coarse search for a global maximum of fp() and
336         * subsequently (b) a numerical optimization using Brent's method (implemented with Apache Commons Math).
337         *
338         * @return the "canonical" start point phase
339         */
340        private double getStartPointPhase() {
341                
342                // the 1-dimensional target function to be optimized:
343                UnivariateFunction fp = new UnivariateFunction() {      // function to be maximized
344                        /** 
345                         * The value returned is the sum of the cross products of the FD pairs,
346                         * with all coefficients rotated to the given start point phase phi.
347                         * TODO: This could be made a lot more efficient!
348                         */
349                        @Override
350                        public double value(double phi) {
351                                double sum = 0;
352                                for (int m = 1; m <= mp; m++) {
353                                        Complex Gm = getCoefficient(-m).rotate(-m * phi);
354                                        Complex Gp = getCoefficient( m).rotate( m * phi);
355                                        sum = sum + Gp.re * Gm.im - Gp.im * Gm.re;      // "cross product" Gp x Gm
356                                }
357                                return sum;
358                        }
359                };
360                
361                // search for the global maximum in coarse steps
362                double cmax = Double.NEGATIVE_INFINITY;
363                int kmax = -1;
364                int K = 25;                                                                     // do K search steps over 0,...,PI
365                for (int k = 0; k < K; k++) {                           // find opt. phi by maximizing fp(phi)
366                        final double phi = Math.PI * k / K;     // phase to evaluate
367                        final double c = fp.value(phi);
368                        if (c > cmax) {
369                                cmax = c;
370                                kmax = k;
371                        }
372                }
373                
374                // final optimize using a BrentOptimizer and previous/next point as the bracket:
375                double minPhi = Math.PI * (kmax - 1) / K;
376                double maxPhi = Math.PI * (kmax + 1) / K;       
377                UnivariateOptimizer optimizer = new BrentOptimizer(1E-4, 1E-6);
378                int maxIter = 20;
379                UnivariatePointValuePair result = optimizer.optimize(
380                                new MaxEval(maxIter),
381                                new UnivariateObjectiveFunction(fp),
382                                GoalType.MAXIMIZE,
383                                new SearchInterval(minPhi, maxPhi)
384                                );
385                double phi0 = result.getPoint();
386                return phi0;    // the canonical start point phase
387        }
388
389        /**
390         * Returns a new rotation invariant Fourier descriptor by applying complex rotation to all coefficients (except
391         * coefficient 0). The original Fourier descriptor is not modified.
392         *
393         * @return a new rotation-normalized Fourier descriptor
394         */
395        public FourierDescriptor makeRotationInvariant() {
396                Complex z = new Complex(0,0);
397                for (int m = 1; m <= mp; m++) {
398                        final double w = 1.0 / m;       // weighting factor emphasizing low-frequency components
399                        z = z.add(getCoefficient(-m).multiply(w));
400                        z = z.add(getCoefficient(+m).multiply(w));
401                }
402                double beta = z.arg();
403                Complex[] G2 = G.clone();
404                for (int i = 1; i < G2.length; i++) {
405                        G2[i] = G[i].rotate(-beta);
406                }
407                
408                return new FourierDescriptor(G2, this.scale);
409        }
410
411        // -----------------------------------------------------------------
412
413        /**
414         * Returns a new translation invariant Fourier descriptor by setting coefficient 0 to zero. The original Fourier
415         * descriptor is not modified. This method is not used for shape invariance calculation, since position is not
416         * shape-relevant.
417         *
418         * @return a new translation-normalized Fourier descriptor
419         */
420        public FourierDescriptor makeTranslationInvariant() {
421                Complex[] G2 = G.clone();
422                G2[0] = Complex.ZERO;
423                return new FourierDescriptor(G2, this.scale);
424        }
425        
426        // ------------------------------------------------------------------
427        // Reconstruction-related methods:
428        // ------------------------------------------------------------------
429
430        /**
431         * Reconstructs the associated 2D shape (closed contour) with N sample points, using all of the Fourier descriptor's
432         * coefficient pairs. The result is returned as an array of {@link Complex} values.
433         *
434         * @param n number of shape points
435         * @return reconstructed shape points
436         */
437        public Complex[] getShapeFull(int n) {
438                return getShapePartial(n, mp);
439        }
440
441        /**
442         * Reconstructs the associated 2D shape (closed contour) with N sample points, using only a subset of the Fourier
443         * descriptor's coefficient pairs. The result is returned as an array of {@link Complex} values.
444         *
445         * @param n number of shape points
446         * @param p the number of (low frequency) coefficient pairs to be used
447         * @return the reconstructed shape points
448         */
449        public Complex[] getShapePartial(int n, int p) {
450                p = Math.min(p, this.mp);
451                Complex[] S = new Complex[n];
452                for (int i = 0; i < n; i++) {
453                        double t = (double) i / n;
454                        S[i] = getShapePointPartial(p, t);
455                }
456                return S;
457        }
458
459        /**
460         * Reconstructs a single space point of the associated shape (closed contour) at the fractional path position t
461         * &isin; [0,1], using all of the Fourier descriptor's coefficient pairs.
462         *
463         * @param t path position
464         * @return the reconstructed shape point
465         */
466        public Complex getShapePointFull(double t) {
467                return getShapePointPartial(mp, t);
468        }
469
470        /**
471         * Reconstructs a single space point of the associated shape (closed contour) at the fractional path position t
472         * &isin; [0,1], using only a subset of the Fourier descriptor's coefficient pairs.
473         *
474         * @param p the number of (low frequency) coefficient pairs to be used
475         * @param t path position &isin; [0,1]
476         * @return the reconstructed shape point
477         */
478        public Complex getShapePointPartial(int p, double t) {
479                p = Math.min(p, mp);
480                double x = G[0].re;
481                double y = G[0].im;
482                for (int m = -p; m <= +p; m++) {
483                        if (m != 0) {
484                                Complex Gm = getCoefficient(m);
485                                double A = scale * Gm.re;
486                                double B = scale * Gm.im;
487                                double phi = 2 * Math.PI * m * t;
488                                double sinPhi = Math.sin(phi);
489                                double cosPhi = Math.cos(phi);
490                                x = x + A * cosPhi - B * sinPhi;
491                                y = y + A * sinPhi + B * cosPhi;
492                        }
493                }
494                return new Complex(x, y);
495        }
496
497        /**
498         * Reconstructs the associated 2D shape (closed contour) with N sample points, using only a single coefficient
499         * pairs. The result is returned as an array of {@link Complex} values.
500         *
501         * @param n number of shape points
502         * @param m the frequency index of the coefficient pair
503         * @return the reconstructed shape points
504         */
505        public Complex[] getShapePair(int n, int m) {
506                m = Math.min(m, this.mp);
507                Complex[] S = new Complex[n];
508                for (int i = 0; i < n; i++) {
509                        double t = (double) i / n;
510                        S[i] = getShapePointPair(m, t);
511                }
512                return S;
513        }
514
515        /**
516         * Returns the spatial point reconstructed from a single DFT coefficient pair G[-m], G[+m] at position t &isin;
517         * [0,1]. Varying t creates points on an ellipse.
518         *
519         * @param m frequency index (coefficient pair number)
520         * @param t contour position &isin; [0,1]
521         * @return reconstructed shape point
522         */
523        public Complex getShapePointPair(int m, double t) {
524                Complex p1 = getShapePointSingle(-m, t);
525                Complex p2 = getShapePointSingle(+m, t);
526                return p1.add(p2);
527        }
528
529
530        /**
531         * Returns the spatial point reconstructed from a single DFT coefficient G[m] with frequency index m at position t
532         * &isin; [0,1]. Varying t creates points on a circle.
533         *
534         * @param m frequency index (pos/neg, 0 &le; m &le; this.mp)
535         * @param t contour position &isin; [0,1]
536         * @return reconstructed shape point
537         */
538        private Complex getShapePointSingle(int m, double t) {
539                Complex Gm = getCoefficient(m);
540                double wm = 2 * Math.PI * m;
541                double Am = Gm.re;
542                double Bm = Gm.im;
543                double cost = Math.cos(wm * t);
544                double sint = Math.sin(wm * t);
545                double xt = Am * cost - Bm * sint;
546                double yt = Bm * cost + Am * sint;
547                return new Complex(xt, yt);
548        }
549
550        /**
551         * <p>
552         * Returns the parameters of the geometric ellipse associated with a single Fourier coefficient pair (G[-m], G[+m]).
553         * See Sec. 26.3.5 of [1] for details. The result is in the form (ra, rb, xc, yc, theta).
554         * </p>
555         * <p>
556         * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction Using Java</em>, 2nd
557         * ed, Springer (2016).
558         * </p>
559         *
560         * @param m the frequency index of the Fourier coefficient pair
561         * @return the ellipse parameters (ra, rb, xc, yc, theta)
562         */
563        public double[] getEllipseParameters(int m) {
564                Complex Gmm = getCoefficient(-m);
565                Complex Gpm = getCoefficient(+m);
566                // see [1], Eqns. (26.50 - 52):
567                double ra = Gmm.abs() + Gpm.abs();                              // = a_m
568                double rb = Math.abs(Gmm.abs() - Gpm.abs());    // = b_m
569                double theta = 0.5 * (Gmm.arg() + Gpm.arg());   // = alpha_m
570                
571                return new double[] {ra, rb, 0, 0, theta};
572        }
573
574        /**
575         * <p>
576         * Returns the ellipse associated with a single Fourier coefficient pair (G[-m], G[+m]) as a (AWT) {@link Shape}
577         * object. The resulting ellipse is centered at (0,0). {@link AffineTransform} may be used to shift the ellipse to
578         * any other position.
579         * </p>
580         *
581         * @param m the frequency index of the Fourier coefficient pair
582         * @return the ellipse ({@link Shape})
583         * @see #getEllipseParameters(int)
584         */
585        public Shape getEllipse(int m) {
586                double[] ep = this.getEllipseParameters(m);     // = (ra, rb, 0, 0, theta)
587                double ra = ep[0]; 
588                double rb = ep[1]; 
589                double theta = ep[4];
590                AffineTransform rot = AffineTransform.getRotateInstance(theta);
591                return rot.createTransformedShape(new Ellipse2D.Double(-ra, -rb, 2 * ra, 2 * rb));
592        }
593
594        
595        // ------------------------------------------------------------------
596        // Distance methods for matching Fourier descriptors:
597        // ------------------------------------------------------------------   
598
599        /**
600         * Returns a L2-type distance between this and another {@link FourierDescriptor} instance comparing the real and
601         * imaginary parts of all coefficient pairs. The zero-frequency coefficients are ignored.
602         *
603         * @param fd2 another Fourier descriptor
604         * @return the resulting distance
605         */
606        public double distanceComplex(FourierDescriptor fd2) {
607                return distanceComplex(fd2, mp);
608        }
609
610        /**
611         * Returns a L2-type distance between this and another {@link FourierDescriptor} instance comparing the real and
612         * imaginary parts of a limited range of (low-frequency) coefficient pairs. The zero-frequency coefficients are
613         * ignored.
614         *
615         * @param fd2 another Fourier descriptor
616         * @param p the number of (low-frequency) coefficient pairs to evaluate
617         * @return the resulting distance
618         */
619        public double distanceComplex(FourierDescriptor fd2, final int p) {
620                if (this.mp < p || fd2.mp < p) {
621                        throw new IllegalArgumentException("insufficient number of Fourier coefficients");
622                }
623                FourierDescriptor fd1 = this;
624                double sum = 0;
625                for (int m = -p; m <= p; m++) {
626                        if (m != 0) {
627                                Complex G1m = fd1.getCoefficient(m);
628                                Complex G2m = fd2.getCoefficient(m);
629                                sum = sum + sqr(G1m.re - G2m.re) + sqr(G1m.im - G2m.im);
630                        }
631                }
632                return Math.sqrt(sum);
633        }
634
635        /**
636         * Returns a L2-type distance between this and another {@link FourierDescriptor} instance comparing the magnitudes
637         * of all coefficient pairs. The zero-frequency coefficients are ignored.
638         *
639         * @param fd2 another Fourier descriptor
640         * @return the resulting distance
641         */
642        public double distanceMagnitude(FourierDescriptor fd2) {
643                return distanceMagnitude(fd2, mp);
644        }
645
646        /**
647         * Returns a L2-type distance between this and another {@link FourierDescriptor} instance comparing the magnitudes
648         * of a limited range of (low-frequency) coefficient pairs. The zero-frequency coefficients are ignored.
649         *
650         * @param fd2 another Fourier descriptor
651         * @param p the number of (low-frequency) coefficient pairs to evaluate
652         * @return the resulting distance
653         */
654        public double distanceMagnitude(FourierDescriptor fd2, final int p) {
655                if (this.mp < p || fd2.mp < p) {
656                        throw new IllegalArgumentException("insufficient number of Fourier coefficients");
657                }
658                FourierDescriptor fd1 = this;
659                double sum = 0;
660                for (int m = -p; m <= p; m++) {
661                        if (m != 0) {
662                                double mag1 = fd1.getCoefficient(m).abs();
663                                double mag2 = fd2.getCoefficient(m).abs();
664                                sum = sum + sqr(mag2 - mag1);
665                        }
666                }
667                return Math.sqrt(sum);
668        }
669        
670        // ---------------------------------------------------------------------
671        
672        @Override
673        public String toString() {
674                return String.format(Locale.US, "%s: mp=%d scale=%.3f", this.getClass().getSimpleName(), mp, scale);
675        }
676        
677        // moved from Utils ------------------------------------
678
679        public static Path2D toPath(Complex[] C) {
680                Path2D path = new Path2D.Float();
681                path.moveTo(C[0].re, C[0].im);
682                for (int i = 1; i < C.length; i++) {
683                        path.lineTo(C[i].re, C[i].im);
684                }
685                path.closePath();
686                return path;
687        }
688
689        public static Complex[] toComplexArray(Pnt2d[] points) {
690                int N = points.length;
691                Complex[] samples = new Complex[N];
692                for (int i = 0; i < N; i++) {
693                        samples[i] = new Complex(points[i].getX(), points[i].getY());
694                }
695                return samples;
696        }
697        
698        // -----------------------------------------------------------------
699        // -----------------------------------------------------------------
700
701//      /**
702//       * For testing: apply shape rotation to this FourierDescriptor (phi in radians)
703//       * 
704//       * @param phi rotation angle
705//       */
706//      public void rotate(double phi) {
707//              rotate(G, phi);
708//      }
709
710//      /**
711//       * For testing: apply shape rotation to this FourierDescriptor (phi in radians)
712//       * 
713//       * @param shape complex point
714//       * @param phi angle
715//       */
716//      public static Complex[] rotateShape(Complex[] shapeOrig, double phi) {
717//              Complex[] shape = shapeOrig.clone();
718//              Complex rot = new Complex(phi);
719//              for (int m = 0; m < shape.length; m++) {
720//                      shape[m] = shape[m].multiply(rot);
721//              }
722//              return shape;
723//      }
724        
725//      public static Complex[] scaleShape(Complex[] shapeOrig, double scale) {
726//              Complex[] shape = shapeOrig.clone();
727//              for (int m = 1; m < shape.length; m++) {
728//                      shape[m] = shape[m].multiply(rot);
729//              }
730//              return shape;
731//      }
732
733//      /**
734//       * Apply a particular start-point phase shift
735//       * 
736//       * @param phi start point phase
737//       * @param mp most positive/negative frequency index
738//       */
739//      private void shiftStartPointPhase(double phi) {
740//              for (int m = -mp; m <= mp; m++) {
741//                      if (m != 0) {
742//                              setCoefficient(m, getCoefficient(m).rotate(m * phi));
743//                      }
744//              }
745//      }
746
747}