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.math.Arithmetic;
013
014/**
015 * <p>
016 * This class defines methods for statistical moment calculations on 2D point sets. See Sec. 8.5 of [1] for details.
017 * This abstract class defines static methods only.
018 * </p>
019 * <p>
020 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
021 * (2022).
022 * </p>
023 *
024 * @author WB
025 * @version 2022/11/17
026 */
027public abstract class Moments2D {       //TODO: make more efficient versions!
028        
029        private Moments2D() {}
030
031        /**
032         * Calculates and returns the ordinary moment of order (p,q) for the specified set of 2D points.
033         *
034         * @param points a set of 2D points
035         * @param p order index p
036         * @param q order index q
037         * @return the moment value
038         */
039        public static double ordinaryMoment(Iterable<Pnt2d> points, int p, int q) {
040                double mpq = 0.0;
041                if (p == 0 && q == 0) { // just count the number of points
042                        for (@SuppressWarnings("unused") Pnt2d pnt : points) {
043                                mpq += 1;
044                        }
045                }
046                else {
047                        for (Pnt2d pnt : points) {
048                                mpq += Math.pow(pnt.getX(), p) * Math.pow(pnt.getY(), q);
049                        }
050                }
051                return mpq;
052        }
053
054        /**
055         * Calculates and returns the central moment of order (p,q) for the specified set of 2D points.
056         *
057         * @param points a set of 2D points
058         * @param p order index p
059         * @param q order index q
060         * @return the moment value
061         */
062        public static double centralMoment(Iterable<Pnt2d> points, int p, int q) {
063                double m00 = ordinaryMoment(points, 0, 0); // region area
064                if (Arithmetic.isZero(m00)) {
065                        throw new RuntimeException("empty point set");
066                }
067                double xc = ordinaryMoment(points, 1, 0) / m00;
068                double yc = ordinaryMoment(points, 0, 1) / m00;
069                double mupq = 0.0;
070                for (Pnt2d pnt : points) {
071                        mupq += Math.pow(pnt.getX() - xc, p) * Math.pow(pnt.getY() - yc, q);
072                }
073                return mupq;
074        }
075
076        /**
077         * Calculates and returns the normalized central moment of order (p,q) for the specified set of 2D points.
078         *
079         * @param points a set of 2D points
080         * @param p order index p
081         * @param q order index q
082         * @return the moment value
083         */
084        public static double normalizedCentralMoment(Iterable<Pnt2d> points, int p, int q) {
085                double m00 = ordinaryMoment(points, 0, 0);
086                double scale = 1.0 / Math.pow(m00, 0.5 * (p + q) + 1);
087                return scale * centralMoment(points, p, q);
088        }
089
090}
091