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 ******************************************************************************/
009package imagingbook.common.geometry.moments;
010
011import imagingbook.common.geometry.basic.Pnt2d;
012import imagingbook.common.geometry.basic.PntUtils;
013import imagingbook.common.math.Complex;
014
015/**
016 * <p>
017 * Naive implementation of Flusser's complex invariant moments [1]. See Sec. 8.6.5 (Eq. 8.51 - 8.54) of [2] for
018 * additional details.
019 * </p>
020 * <p>
021 * [1] J. Flusser, B. Zitova, and T. Suk. "Moments and Moment Invariants in Pattern Recognition". John Wiley and Sons
022 * (2009).
023 * <br>
024 * [2] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
025 * (2022).
026 * </p>
027 *
028 * @author WB
029 * @version 2022/12/28
030 */
031public class FlusserMoments {
032
033    private final Iterable<Pnt2d> points;
034    private final int n;
035    private final double xc, yc;
036
037    public FlusserMoments(Iterable<Pnt2d> points) {
038        this.points = points;
039        double sx = 0;
040        double sy = 0;
041        int i = 0;
042        for (Pnt2d p : points) {
043            sx = sx + p.getX();
044            sy = sy + p.getY();
045            i++;
046        }
047        this.n = i;
048        if (n == 0) {
049            throw new IllegalArgumentException("at least one point is required");
050        }
051        this.xc = sx / n;
052        this.yc = sy / n;
053    }
054    /**
055     * Returns the scale-normalized complex moment of order (p,q) for the 2D point set associated with
056     * this {@link FlusserMoments} instance.
057     *
058     * @param p order index p
059     * @param q order index q
060     * @return the complex moment of order (p,q)
061     */
062    public Complex getComplexMoment(int p, int q) {
063        Complex sum = new Complex(0, 0);
064        for (Pnt2d pnt : points) {
065            double x = pnt.getX() - xc;
066            double y = pnt.getY() - yc;
067            Complex zp = (x == 0.0 && y == 0 && p == 0) ?                       // beware: 0^0 is undefined!
068                    Complex.ZERO : new Complex(x, y).pow(p);
069            Complex zq = (x == 0.0 && y == 0 && q == 0) ?
070                    Complex.ZERO : new Complex(x, -y).pow(q);
071            sum = sum.add(zp.multiply(zq));
072        }
073        checkForNaN(sum);
074        return sum;
075    }
076
077    /**
078     * Returns the scale-normalized complex moment of order (p,q) for the specified set of 2D points.
079     *
080     * @param p order index p
081     * @param q order index q
082     * @return the complex moment of order (p,q)
083     */
084    public Complex getScaleNormalizedMoment(int p, int q) {
085        Complex cpq = getComplexMoment(p, q);
086        return cpq.multiply(1.0 / Math.pow(n, 0.5 * (p + q) + 1));
087    }
088
089    /**
090     * Calculates and returns a vector of 11 invariant moments for the specified set of 2D points.
091     *
092     * @return the vector of invariant moments
093     */
094    public double[] getInvariantMoments() {
095        Complex c11 = getScaleNormalizedMoment(1, 1);
096        Complex c12 = getScaleNormalizedMoment(1, 2);
097        Complex c21 = getScaleNormalizedMoment(2, 1);
098        Complex c20 = getScaleNormalizedMoment(2, 0);
099        Complex c22 = getScaleNormalizedMoment(2, 2);
100        Complex c30 = getScaleNormalizedMoment(3, 0);
101        Complex c31 = getScaleNormalizedMoment(3, 1);
102        Complex c40 = getScaleNormalizedMoment(4, 0);
103        double p1 = c11.getRe();
104        double p2 = c21.multiply(c12).getRe();
105        double p3 = c20.multiply(c12.pow(2)).getRe();
106        double p4 = c20.multiply(c12.pow(2)).getIm();
107        double p5 = c30.multiply(c12.pow(3)).getRe();
108        double p6 = c30.multiply(c12.pow(3)).getIm();
109        double p7 = c22.getRe();
110        double p8 = c31.multiply(c12.pow(2)).getRe();
111        double p9 = c31.multiply(c12.pow(2)).getIm();
112        double p10 = c40.multiply(c12.pow(4)).getRe();
113        double p11 = c40.multiply(c12.pow(4)).getIm();
114        return new double[] {p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11};
115    }
116
117    private void checkForNaN(Complex z) {
118        if (z.isNaN()) {
119            throw new RuntimeException("NaN encountered in complex quantity");
120        }
121    }
122
123}