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-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009package imagingbook.pub.fd;
010
011import imagingbook.lib.math.Arithmetic;
012import imagingbook.lib.math.Complex;
013
014import java.awt.geom.Point2D;
015
016
017/**
018 * Subclass of FourierDescriptor whose constructors assume
019 * that input polygons are uniformly sampled.
020 * 
021 * @author W. Burger
022 * @version 2015/08/13
023 */
024public class FourierDescriptorUniform extends FourierDescriptor {
025        
026        /**
027         * Creates a new Fourier descriptor from a uniformly sampled polygon V
028         * with the maximum number of Fourier coefficient pairs.
029         * 
030         * @param V polygon
031         */
032        public FourierDescriptorUniform(Point2D[] V) {
033                g = makeComplex(V);
034                G = DFT(g);
035        }
036        
037        
038        /**
039         * Creates a new Fourier descriptor from a uniformly sampled polygon V
040         * with Mp coefficient pairs.
041         * 
042         * @param V polygon
043         * @param Mp number of coefficient pairs
044         */
045        public FourierDescriptorUniform(Point2D[] V, int Mp) {
046                g = makeComplex(V);
047                G = DFT(g, 2 * Mp + 1);
048        }
049        
050        // -------------------------------------------------------------------
051        
052
053        /**
054         * DFT with the resulting spectrum (signal, if inverse) of the same length
055         * as the input vector g. Not using sin/cos tables.
056         * 
057         * @param g signal vector
058         * @return DFT spectrum
059         */
060        private Complex[] DFT(Complex[] g) {
061                int M = g.length;
062//              double[] cosTable = makeCosTable(M);    // cosTable[m] == cos(2*pi*m/M)
063//              double[] sinTable = makeSinTable(M);
064                Complex[] G = new Complex[M];
065                double s = 1.0 / M; //common scale factor (fwd/inverse differ!)
066                for (int m = 0; m < M; m++) {
067                        double Am = 0;
068                        double Bm = 0;
069                        for (int k = 0; k < M; k++) {
070                                double x = g[k].re();
071                                double y = g[k].im();
072                                double phi = 2 * Math.PI * m * k / M;
073                                double cosPhi = Math.cos(phi);
074                                double sinPhi = Math.sin(phi);
075                                Am = Am + x * cosPhi + y * sinPhi;
076                                Bm = Bm - x * sinPhi + y * cosPhi;
077                        }
078                        G[m] = new Complex(s * Am, s * Bm);
079                }
080                return G;
081        }
082        
083
084        /**
085         * As above, but the length P of the resulting spectrum (signal, if inverse) 
086         * is explicitly specified.
087         * @param g signal vector
088         * @param P length of the resulting  DFT spectrum
089         * @return DFT spectrum
090         */
091        private Complex[] DFT(Complex[] g, int P) {
092                int M = g.length;
093//              double[] cosTable = makeCosTable(M);    // cosTable[m] == cos(2*pi*m/M)
094//              double[] sinTable = makeSinTable(M);
095                Complex[] G = new Complex[P];
096                double s = 1.0/M; //common scale factor (fwd/inverse differ!)
097                for (int m = P/2-P+1; m <= P/2; m++) {
098                        double Am = 0;
099                        double Bm = 0;
100                        for (int k = 0; k < M; k++) {
101                                double x = g[k].re();
102                                double y = g[k].im();
103                                //int mk = (m * k) % M; double phi = 2 * Math.PI * mk / M;
104                                double phi = 2 * Math.PI * m * k / M;   
105                                double cosPhi = Math.cos(phi);
106                                double sinPhi = Math.sin(phi);
107                                Am = Am + x * cosPhi + y * sinPhi;
108                                Bm = Bm - x * sinPhi + y * cosPhi;
109                        }
110                        G[Arithmetic.mod(m, P)] = new Complex(s * Am, s * Bm);
111                }
112                return G;
113        }
114
115//      private double[] makeCosTable(int M) {
116//              double[] cosTab = new double[M];
117//              for (int m = 0; m < M; m++) {
118//                      cosTab[m] = Math.cos(2 * Math.PI * m / M);
119//              }
120//              return cosTab;
121//      }
122
123//      private double[] makeSinTable(int M) {
124//              double[] sinTab = new double[M];
125//              for (int m = 0; m < M; m++) {
126//                      sinTab[m] = Math.sin(2 * Math.PI * m / M);
127//              }
128//              return sinTab;
129//      }
130
131}