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.fd; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.math.Arithmetic; 013import imagingbook.common.math.Complex; 014 015/** 016 * <p> 017 * Defines static methods to create Fourier descriptors from uniformly sampled point sequences. See Ch. 26 of [1] for 018 * additional details. 019 * </p> 020 * <p> 021 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction Using Java</em>, 2nd ed, 022 * Springer (2016). 023 * </p> 024 * 025 * @author WB 026 * @version 2022/10/24 027 */ 028public abstract class FourierDescriptorUniform { 029 030 private FourierDescriptorUniform() {} 031 032 /** 033 * Creates a new Fourier descriptor from a uniformly sampled polygon V with the maximum number of Fourier 034 * coefficient pairs. The length of the resulting DFT spectrum equals V.length. 035 * 036 * @param V a polygon 037 * @return a new {@link FourierDescriptorUniform} instance 038 */ 039 public static FourierDescriptor from(Pnt2d[] V) { 040 return from(FourierDescriptor.toComplexArray(V)); 041 } 042 043 /** 044 * Creates a new Fourier descriptor from a uniformly sampled polygon V with Mp coefficient pairs. The length of the 045 * resulting DFT spectrum is 2 * mp + 1, i.e., it must be assured that Mp < (V.length - 1) ÷ 2. 046 * 047 * @param V a polygon 048 * @param mp number of coefficient pairs 049 * @return a new {@link FourierDescriptorUniform} instance 050 */ 051 public static FourierDescriptor from(Pnt2d[] V, int mp) { 052 return from(FourierDescriptor.toComplexArray(V), mp); 053 } 054 055 056 // ------------------------------------------------------------------- 057 058 /** 059 * Creates a new Fourier descriptor from a uniformly sampled polygon V with the maximum number of Fourier 060 * coefficient pairs. The length of the resulting DFT spectrum equals V.length. 061 * 062 * @param V a polygon 063 * @return a new {@link FourierDescriptorUniform} instance 064 */ 065 public static FourierDescriptor from(Complex[] V) { 066 return new FourierDescriptor(DFT(V)); 067 } 068 069 /** 070 * Creates a new Fourier descriptor from a uniformly sampled polygon V with Mp coefficient pairs. The length of the 071 * resulting DFT spectrum is 2 * mp + 1, i.e., it must be assured that Mp < (V.length - 1) ÷ 2. 072 * 073 * @param V a sequence of uniformly-spaced 2D sample points 074 * @param mp number of coefficient pairs 075 * @return a new {@link FourierDescriptorUniform} instance 076 */ 077 public static FourierDescriptor from(Complex[] V, int mp) { 078 if (mp > (V.length - 1) / 2) { 079 throw new IllegalArgumentException("number of Fourier pairs (mp) may not be greater than " 080 + ((V.length - 1) / 2)); 081 } 082 return new FourierDescriptor(DFT(V, 2 * mp + 1)); 083 } 084 085 // ------------------------------------------------------------------- 086 087 /** 088 * DFT with the resulting spectrum (signal, if inverse) of the same length as the input vector g. Not using sin/cos 089 * tables. 090 * 091 * @param g signal vector 092 * @return DFT spectrum 093 */ 094 private static Complex[] DFT(Complex[] g) { 095 int M = g.length; 096// double[] cosTable = makeCosTable(M); // cosTable[m] == cos(2*pi*m/M) 097// double[] sinTable = makeSinTable(M); // sinTable[m] == sin(2*pi*m/M) 098 Complex[] G = new Complex[M]; 099 double s = 1.0 / M; //common scale factor (fwd/inverse differ!) 100 for (int m = 0; m < M; m++) { 101 double Am = 0; 102 double Bm = 0; 103 for (int k = 0; k < M; k++) { 104 double x = g[k].re; 105 double y = g[k].im; 106 double phi = 2 * Math.PI * m * k / M; 107 double cosPhi = Math.cos(phi); 108 double sinPhi = Math.sin(phi); 109 Am = Am + x * cosPhi + y * sinPhi; 110 Bm = Bm - x * sinPhi + y * cosPhi; 111 } 112 G[m] = new Complex(s * Am, s * Bm); 113 } 114 return G; 115 } 116 117 118 /** 119 * As {@link #DFT(Complex[])}, with the length of the resulting spectrum (or signal, if inverse) explicitly 120 * specified. 121 * 122 * @param g signal vector 123 * @param MM length of the resulting DFT spectrum 124 * @return DFT spectrum 125 */ 126 private static Complex[] DFT(Complex[] g, int MM) { 127 int M = g.length; 128 if (MM > M) { 129 throw new IllegalArgumentException("truncated spectrum must be shorter than original MM = " + MM); 130 } 131// double[] cosTable = makeCosTable(M); // cosTable[m] == cos(2*pi*m/M) 132// double[] sinTable = makeSinTable(M); // sinTable[m] == sin(2*pi*m/M) 133 Complex[] G = new Complex[MM]; 134 double s = 1.0 / M; //common scale factor (fwd/inverse differ!) 135 136// for (int j = Mp/2-Mp+1; j <= Mp/2; j++) { 137 for (int j = -MM/2; j <= (MM-1)/2; j++) { 138 double Am = 0; 139 double Bm = 0; 140 for (int k = 0; k < M; k++) { 141 double x = g[k].re; 142 double y = g[k].im; 143 //int mk = (m * k) % M; double phi = 2 * Math.PI * mk / M; 144 double phi = 2 * Math.PI * j * k / M; 145 double cosPhi = Math.cos(phi); 146 double sinPhi = Math.sin(phi); 147 Am = Am + x * cosPhi + y * sinPhi; 148 Bm = Bm - x * sinPhi + y * cosPhi; 149 } 150 G[Arithmetic.mod(j, MM)] = new Complex(s * Am, s * Bm); 151 } 152 return G; 153 } 154 155// private double[] makeCosTable(int M) { 156// double[] cosTab = new double[M]; 157// for (int m = 0; m < M; m++) { 158// cosTab[m] = Math.cos(2 * Math.PI * m / M); 159// } 160// return cosTab; 161// } 162 163// private double[] makeSinTable(int M) { 164// double[] sinTab = new double[M]; 165// for (int m = 0; m < M; m++) { 166// sinTab[m] = Math.sin(2 * Math.PI * m / M); 167// } 168// return sinTab; 169// } 170 171}