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.mappings.linear;
011
012import imagingbook.common.geometry.basic.Pnt2d;
013import imagingbook.common.geometry.basic.Pnt2d.PntDouble;
014import imagingbook.common.geometry.fitting.points.AffineFit2d;
015import imagingbook.common.math.Arithmetic;
016import imagingbook.common.math.Matrix;
017
018/**
019 * <p>
020 * This class represents an affine transformation in 2D, which can be defined by three pairs of corresponding points. It
021 * can be assumed that every instance of this class is indeed an affine mapping. Instances are immutable. See Secs.
022 * 21.1.3 and 21.3.1 of [1] for details.
023 * </p>
024 * <p>
025 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
026 * (2022).
027 * </p>
028 *
029 * @author WB
030 */
031public class AffineMapping2D extends ProjectiveMapping2D {
032
033        /**
034         * Creates an affine 2D mapping from two sequences of corresponding points. If 3 point pairs are specified, the
035         * mapping is exact, otherwise a minimum least-squares fit is calculated.
036         *
037         * @param P the source points
038         * @param Q the target points
039         * @return a new {@link AffineMapping2D} instance for the two point sets
040         */
041        public static AffineMapping2D fromPoints(Pnt2d[] P, Pnt2d[] Q) {
042                if (P.length != Q.length) {
043                        throw new IllegalArgumentException("point sets P, Q must have the same size");
044                }
045                if (P.length < 3) {
046                        throw new IllegalArgumentException("at least 3 point pairs are required");
047                }
048                if (P.length == 3) {
049                        // exact fit
050                        return fromPoints(P[0], P[1], P[2], Q[0], Q[1], Q[2]);
051                }
052                else {
053                        // minimum least-squares fit
054                        AffineFit2d fit = new AffineFit2d(P, Q);
055                        return new AffineMapping2D(fit.getTransformationMatrix());
056                }
057        }
058
059        /**
060         * <p>
061         * Creates an affine mapping from 3 pairs of corresponding 2D points (p0, p1, p2) &rarr; (q0, q1, q2). The solution
062         * is found in closed form (see [1], Sec. 21.1.3, eq. 21.31).
063         * </p>
064         * <p>
065         * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed,
066         * Springer (2022).
067         * </p>
068         *
069         * @param p0 point 1 of source triangle A
070         * @param p1 point 2 of source triangle A
071         * @param p2 point 3 of source triangle A
072         * @param q0 point 1 of source triangle B
073         * @param q1 point 2 of source triangle B
074         * @param q2 point 3 of source triangle B
075         * @return a new {@link AffineMapping2D} instance
076         */
077        public static AffineMapping2D fromPoints(Pnt2d p0, Pnt2d p1, Pnt2d p2, Pnt2d q0, Pnt2d q1, Pnt2d q2) {
078                final double px0 = p0.getX(), px1 = p1.getX(), px2 = p2.getX();
079                final double py0 = p0.getY(), py1 = p1.getY(), py2 = p2.getY();
080                final double qx0 = q0.getX(), qx1 = q1.getX(), qx2 = q2.getX();
081                final double qy0 = q0.getY(), qy1 = q1.getY(), qy2 = q2.getY();
082
083                final double d = px0 * (py2 - py1) + px1 * (py0 - py2) + px2 * (py1 - py0); 
084                if (Arithmetic.isZero(d)) {
085                        throw new ArithmeticException("affine mapping is undefined (d=0)");
086                }
087                double a00 = (py0 * (qx1 - qx2) + py1 * (qx2 - qx0) + py2 * (qx0 - qx1)) / d;
088                double a01 = (px0 * (qx2 - qx1) + px1 * (qx0 - qx2) + px2 * (qx1 - qx0)) / d;
089                double a10 = (py0 * (qy1 - qy2) + py1 * (qy2 - qy0) + py2 * (qy0 - qy1)) / d;
090                double a11 = (px0 * (qy2 - qy1) + px1 * (qy0 - qy2) + px2 * (qy1 - qy0)) / d;
091                
092                double a02 = (px0*(py2*qx1-py1*qx2) + px1*(py0*qx2-py2*qx0) + px2*(py1*qx0-py0*qx1)) / d;
093                double a12 = (px0*(py2*qy1-py1*qy2) + px1*(py0*qy2-py2*qy0) + px2*(py1*qy0-py0*qy1)) / d;
094                
095                return new AffineMapping2D(a00, a01, a02, a10, a11, a12);
096        }
097        
098//      /**
099//       * Creates an affine mapping from an arbitrary 2D triangle A to another triangle B.
100//       * @param A the first triangle
101//       * @param B the second triangle
102//       * @return a new affine mapping
103//       * @deprecated
104//       */
105//      public static AffineMapping2D from3Points(Point[] A, Point[] B) {
106//              return AffineMapping2D.from3Points(A[0], A[1], A[2], B[0], B[1], B[2]);
107//      }
108        // ---------------------------------------------------------------------------
109        
110        /**
111         * Creates the identity mapping.
112         */
113        public AffineMapping2D() {
114                super();
115        }
116
117        /**
118         * Creates a linear mapping from a transformation matrix A, which must be at least of size 2 x 3. The elements of A
119         * are copied into a 3x3 identity matrix. If A is larger than 2 x 3, the remaining elements are ignored.
120         *
121         * @param A a 2x3 (or larger) matrix
122         */
123        public AffineMapping2D(double[][] A) {
124                //super(A[0][0], A[0][1], A[0][2], A[1][0], A[1][1], A[1][2], 0, 0);
125                super(extractAffineMatrix(A));
126        }
127        
128        private static double[][] extractAffineMatrix(double[][] A) {
129                double[][] M = Matrix.idMatrix(3);
130                final int m = Math.min(2, A.length);    // max. 2 rows
131                for (int i = 0; i < m; i++) {
132                        final int n = Math.min(3, A[i].length); // max. 3 columns
133                        for (int j = 0; j < n; j++) {
134                                M[i][j] = A[i][j];
135                        }
136                }
137                return M;
138        }
139
140        /**
141         * Creates an affine mapping from the specified matrix elements.
142         *
143         * @param a00 matrix element A_00
144         * @param a01 matrix element A_01
145         * @param a02 matrix element A_02
146         * @param a10 matrix element A_10
147         * @param a11 matrix element A_11
148         * @param a12 matrix element A_12
149         */
150        public AffineMapping2D (
151                        double a00, double a01, double a02, 
152                        double a10, double a11, double a12) {
153                super(a00, a01, a02, a10, a11, a12, 0, 0);
154        }
155
156        /**
157         * Creates a new affine mapping from an existing affine mapping.
158         * @param m an affine mapping
159         */
160        public AffineMapping2D(AffineMapping2D m) {
161                this(m.a00, m.a01, m.a02, m.a10, m.a11, m.a11);
162        }
163        
164        // ----------------------------------------------------------
165
166
167        /**
168         * Checks if the given linear mapping could be affine, i.e. if the bottom row of its transformation matrix is (0, 0,
169         * 1). Note that this is a necessary but not sufficient requirement. The threshold {@link Arithmetic#EPSILON_DOUBLE}
170         * is used in this check.
171         *
172         * @param map a linear mapping
173         * @return true if the mapping could be affine
174         */
175        public static boolean isAffine(LinearMapping2D map) {
176                final double tol = Arithmetic.EPSILON_DOUBLE; // max. deviation for 0/1 values
177                if (Math.abs(map.a20) > tol) return false;
178                if (Math.abs(map.a21) > tol) return false;
179                if (Math.abs(map.a22 - 1.0) > tol) return false;
180                return true;
181        }
182
183        /**
184         * Concatenates this mapping A with another affine mapping B and returns a new mapping C, such that C(x) = B(A(x)).
185         *
186         * @param B the second mapping
187         * @return the concatenated affine mapping
188         */
189        public AffineMapping2D concat(AffineMapping2D B) {      // use some super method instead?
190                double[][] C = Matrix.multiply(B.getTransformationMatrix(), this.getTransformationMatrix());
191                return new AffineMapping2D(C[0][0], C[0][1], C[0][2], C[1][0], C[1][1], C[1][2]);
192        }
193        
194        @Override       // more efficient than generic version
195        public Pnt2d applyTo(Pnt2d pnt) {
196                final double x = pnt.getX();
197                final double y = pnt.getY();
198                double x1 = (a00 * x + a01 * y + a02);
199                double y1 = (a10 * x + a11 * y + a12);
200                return PntDouble.from(x1, y1);
201        }
202
203        /**
204         * {@inheritDoc} Note that inverting an affine transformation always yields another affine transformation.
205         */
206        @Override
207        public AffineMapping2D getInverse() {
208                double det = a00 * a11 - a01 * a10;
209                double b00 = a11 / det;
210                double b01 = -a01 / det;
211                double b02 = (a01 * a12 - a02 * a11) / det;
212                double b10 = -a10 / det;
213                double b11 = a00 / det;
214                double b12 = (a02 * a10 - a00 * a12) / det;
215                return new AffineMapping2D(b00, b01, b02, b10, b11, b12);
216        }
217
218        /**
219         * {@inheritDoc}
220         * @return a new affine mapping
221         */
222        @Override
223        public AffineMapping2D duplicate() {
224                return new AffineMapping2D(this);
225        }
226
227        @Override
228        public double[][] getJacobian(Pnt2d xy) {
229                final double x = xy.getX();
230                final double y = xy.getY();
231                return new double[][]
232                        {{x, y, 0, 0, 1, 0},
233                         {0, 0, x, y, 0, 1}};
234        }
235        
236        // ----------------------------------------------------------------------
237        
238//      /**
239//       * For testing only.
240//       * @param args ignored
241//       */
242//      public static void main(String[] args) {
243//              PrintPrecision.set(6);
244//              double[][] A = 
245//                      {{-2, 4, -3}, 
246//                      {3, 7, 2}, 
247//                      {0, 0, 1}};
248//              System.out.println("a = \n" + Matrix.toString(A));
249//              System.out.println();
250//              double[][] ai = Matrix.inverse(A);
251//              System.out.println("ai = \n" + Matrix.toString(ai));
252//              
253//              LinearMapping2D Ai = new LinearMapping2D(ai);
254//              System.out.println("Ai is affine: " + isAffine(Ai));
255//              
256//              double[][] I = Matrix.multiply(A, ai);
257//              System.out.println("\ntest: should be the  identity matrix: = \n" + Matrix.toString(I));
258//              
259//              double[][] B = 
260//                      {{-2, 4, -3}, 
261//                      {3, 7, 2}};
262//              
263//              LinearMapping2D am = new AffineMapping2D(B);
264//              System.out.println("an = \n" + Matrix.toString(am.getTransformationMatrix()));
265//      }
266
267}
268
269
270
271