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.ProjectiveFit2d;
015import imagingbook.common.geometry.mappings.Jacobian;
016import imagingbook.common.math.Arithmetic;
017import imagingbook.common.math.Matrix;
018import imagingbook.common.math.PrintPrecision;
019
020/**
021 * <p>
022 * This class represents a projective transformation in 2D (also known as a "homography"). It can be specified uniquely
023 * by four pairs of corresponding points. It can be assumed that every instance of this class is indeed a projective
024 * mapping. See Secs. 21.1.4 and 21.3.1 of [1] for details.
025 * </p>
026 * <p>
027 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
028 * (2022).
029 * </p>
030 *
031 * @author WB
032 */
033public class ProjectiveMapping2D extends LinearMapping2D implements Jacobian {
034        
035        //  static methods -----------------------------------------------------
036
037        /**
038         * Creates a projective 2D mapping from two sequences of corresponding points. If 4 point pairs are specified, the
039         * mapping is exact, otherwise a minimum least-squares fit is calculated.
040         *
041         * @param P the source points
042         * @param Q the target points
043         * @return a new projective mapping for the two point sets
044         */
045        public static ProjectiveMapping2D fromPoints(Pnt2d[] P, Pnt2d[] Q) {
046                ProjectiveFit2d fit = new ProjectiveFit2d(P, Q);
047                return new ProjectiveMapping2D(fit.getTransformationMatrix());
048        }
049        
050//      /**
051//       * Creates the most specific linear mapping from two sequences of corresponding
052//       * 2D points.
053//       * @param P first point sequence
054//       * @param Q second point sequence
055//       * @return a linear mapping derived from point correspondences
056//       */
057//      public static ProjectiveMapping2D makeMapping(Point[] P, Point[] Q) {
058//              int n = Math.min(P.length, Q.length);
059//              if (n < 1) {
060//                      throw new IllegalArgumentException("cannot create a mapping from zero points");
061//              }
062//              else if (n == 1) {
063//                      return new Translation2D(P[0], Q[0]);
064//              }
065//              else if (n == 2) {      // TODO: similarity transformation?
066//                      throw new UnsupportedOperationException("makeMapping: don't know yet how to handle 2 points");
067//              }
068//              else if (n == 3) {
069//                      return AffineMapping2D.fromPoints(P, Q);
070//              }
071//              else if (n == 4) {
072//                      return ProjectiveMapping2D.from4Points(P, Q);
073//              }
074//              else {  // n > 4 (over-determined)
075//                      return ProjectiveMapping2D.fromNPoints(P, Q);
076//              }
077//      }
078        
079//      /**
080//       * Creates the projective mapping from the unit square S to
081//       * the arbitrary quadrilateral P, specified by four points.
082//       * @param p0 point 0
083//       * @param p1 point 1
084//       * @param p2 point 2
085//       * @param p3 point 3
086//       * @return a new projective mapping
087//       */
088//      public static ProjectiveMapping2D fromUnitSquareTo4Points(Point p0, Point p1, Point p2, Point p3) {
089//              double x0 = p0.getX(), x1 = p1.getX(), x2 = p2.getX(), x3 = p3.getX();
090//              double y0 = p0.getY(), y1 = p1.getY(), y2 = p2.getY(), y3 = p3.getY();
091//              double S = (x1 - x2) * (y3 - y2) - (x3 - x2) * (y1 - y2);
092//              if (Arithmetic.isZero(S)) {
093//                      throw new ArithmeticException("fromUnitSquareTo4Points(): division by zero!");
094//              }
095//              double a20 = ((x0 - x1 + x2 - x3) * (y3 - y2) - (y0 - y1 + y2 - y3) * (x3 - x2)) / S;
096//              double a21 = ((y0 - y1 + y2 - y3) * (x1 - x2) - (x0 - x1 + x2 - x3) * (y1 - y2)) / S;
097//              double a00 = x1 - x0 + a20 * x1;
098//              double a01 = x3 - x0 + a21 * x3;
099//              double a02 = x0;
100//              double a10 = y1 - y0 + a20 * y1;
101//              double a11 = y3 - y0 + a21 * y3;
102//              double a12 = y0;
103//              return new ProjectiveMapping2D(a00, a01, a02, a10, a11, a12, a20, a21);
104//      }
105        
106//      /**
107//       * Creates a new projective mapping between arbitrary two quadrilaterals P, Q.
108//       * @param p0 point 0 of source quad P
109//       * @param p1 point 1 of source quad P
110//       * @param p2 point 2 of source quad P
111//       * @param p3 point 3 of source quad P
112//       * @param q0 point 0 of target quad Q
113//       * @param q1 point 1 of target quad Q
114//       * @param q2 point 2 of target quad Q
115//       * @param q3 point 3 of target quad Q
116//       * @return a new projective mapping
117//       */
118//      public static ProjectiveMapping2D from4Points(
119//                      Point p0, Point p1, Point p2, Point p3, 
120//                      Point q0, Point q1, Point q2, Point q3) {
121//              ProjectiveMapping2D T1 = ProjectiveMapping2D.fromUnitSquareTo4Points(p0, p1, p2, p3);
122//              ProjectiveMapping2D T2 = ProjectiveMapping2D.fromUnitSquareTo4Points(q0, q1, q2, q3);
123//              ProjectiveMapping2D T1i = T1.getInverse();
124//              return T1i.concat(T2);          
125//      }
126        
127//      /**
128//       * Creates a new projective mapping between arbitrary two quadrilaterals P, Q.
129//       * @param P source quad
130//       * @param Q target quad
131//       * @return a new projective mapping
132//       */
133//      public static final ProjectiveMapping2D from4Points(Point[] P, Point[] Q) {
134//              return ProjectiveMapping2D.from4Points(P[0], P[1], P[2], P[3], Q[0], Q[1], Q[2], Q[3]);
135//      }
136        
137//      /**
138//       * Maps between n &gt; 4 point pairs, finds a least-squares solution
139//       * for the homography parameters.
140//       * NOTE: this is UNFINISHED code! check against DLT estimation of homography
141//       * @param P sequence of points (source)
142//       * @param Q sequence of points (target)
143//       * @return a new projective mapping
144//       */
145//      public static ProjectiveMapping2D fromNPoints(Point[] P, Point[] Q) {
146//              final int n = P.length;
147//              if (n < 4) {
148//                      throw new IllegalArgumentException(ProjectiveMapping2D.class.getName() + ": fromNPoints() needs at least 4 points pairs");
149//              }
150//              double[] ba = new double[2 * n];
151//              double[][] Ma = new double[2 * n][];
152//              for (int i = 0; i < n; i++) {
153//                      double x = P[i].getX();
154//                      double y = P[i].getY();
155//                      double u = Q[i].getX();
156//                      double v = Q[i].getY();
157//                      ba[2 * i + 0] = u;
158//                      ba[2 * i + 1] = v;
159//                      Ma[2 * i + 0] = new double[] { x, y, 1, 0, 0, 0, -u * x, -u * y };
160//                      Ma[2 * i + 1] = new double[] { 0, 0, 0, x, y, 1, -v * x, -v * y };
161//              }
162//              
163//              RealMatrix M = MatrixUtils.createRealMatrix(Ma);
164//              RealVector b = MatrixUtils.createRealVector(ba);
165//              DecompositionSolver solver = new SingularValueDecomposition(M).getSolver();
166//              RealVector h = solver.solve(b);
167//              double a00 = h.getEntry(0);
168//              double a01 = h.getEntry(1);
169//              double a02 = h.getEntry(2);
170//              double a10 = h.getEntry(3);
171//              double a11 = h.getEntry(4);
172//              double a12 = h.getEntry(5);
173//              double a20 = h.getEntry(6);
174//              double a21 = h.getEntry(7);
175//              return new ProjectiveMapping2D(a00, a01, a02, a10, a11, a12, a20, a21);
176//      }
177        
178        //  constructors -----------------------------------------------------
179        
180        /**
181         * Creates the identity mapping.
182         */
183        public ProjectiveMapping2D() {
184                super();
185        }
186
187        /**
188         * Creates a projective mapping from a transformation matrix A, of arbitrary size. All relevant elements of A are
189         * inserted into actual 3x3 transformation matrix, except element (2,2) which is ignored and always set to 1 (to
190         * keep the transformation projective). If A is larger than 3 x 3, the remaining elements are ignored.
191         *
192         * @param A a transformation matrix of arbitrary size
193         */
194        public ProjectiveMapping2D(double[][] A) {
195                //this(new LinearMapping2D(A));
196                super(extractProjectiveMatrix(A));
197        }
198        
199        private static double[][] extractProjectiveMatrix(double[][] A) {
200                double[][] M = Matrix.idMatrix(3);
201                final int m = Math.min(3, A.length);    // max. 3 rows
202                for (int i = 0; i < m; i++) {
203                        final int n = Math.min(3, A[i].length); // max. 3 columns
204                        for (int j = 0; j < n; j++) {
205                                M[i][j] = A[i][j];
206                        }
207                }
208                M[2][2] = 1;
209                return M;
210        }
211
212
213        /**
214         * Creates a projective mapping from the specified matrix elements.
215         *
216         * @param a00 matrix element A_00
217         * @param a01 matrix element A_01
218         * @param a02 matrix element A_02
219         * @param a10 matrix element A_10
220         * @param a11 matrix element A_11
221         * @param a12 matrix element A_12
222         * @param a20 matrix element A_20
223         * @param a21 matrix element A_21
224         */
225        public ProjectiveMapping2D(
226                        double a00, double a01, double a02, 
227                        double a10, double a11, double a12, 
228                        double a20, double a21) {
229                super(a00, a01, a02, a10, a11, a12, a20, a21, 1);
230        }
231
232        /**
233         * Creates a projective mapping from any linear mapping. The transformation matrix gets normalized to a22 = 1.
234         *
235         * @param m a linear mapping
236         */
237        public ProjectiveMapping2D(LinearMapping2D m) {
238                this(m.normalize());
239        }
240        
241        /**
242         * Creates a projective mapping from an existing instance.
243         * @param m a projective mapping
244         */
245        public ProjectiveMapping2D(ProjectiveMapping2D m) {
246                //this(m.getTransformationMatrix());
247                this(m.a00, m.a01, m.a02, m.a10, m.a11, m.a12, m.a20, m.a21);
248        }
249        
250        // ----------------------------------------------------------
251
252        /**
253         * Concatenates this mapping A with another linear mapping B and returns a new mapping C, such that C(x) = B(A(x)).
254         *
255         * @param B the second mapping
256         * @return the concatenated mapping
257         */
258        public ProjectiveMapping2D concat(ProjectiveMapping2D B) {
259                LinearMapping2D C = LinearMapping2D.concatenate(B, this);
260                return new ProjectiveMapping2D(C);
261        }
262        
263        /**
264         * {@inheritDoc}
265         * @return a new projective mapping
266         */
267        @Override
268        public ProjectiveMapping2D duplicate() {
269                return new ProjectiveMapping2D(this);
270        }
271
272        /**
273         * {@inheritDoc} Note that the inverse A' of a projective transformation matrix A is again a linear transformation
274         * but its a'2' element is generally not 1. Scaling A' to A'' = A' / a22' yields a projective transformation that
275         * reverses A. While A * A' = I, the result of A * A'' is a scaled identity matrix.
276         *
277         * @return the inverse projective transformation
278         */
279        @Override
280        public ProjectiveMapping2D getInverse() {
281                return new ProjectiveMapping2D(super.getInverse());
282        }
283        
284        @Override
285        public double[][] getJacobian(Pnt2d xy) {
286                // see Baker 2003 "20 Years" Part 1, Eq. 99 (p. 46)
287                final double x = xy.getX();
288                final double y = xy.getY();
289                final double a = a00 * x + a01 * y + a02;       // = alpha
290                final double b = a10 * x + a11 * y + a12;       // = beta
291                final double c = a20 * x + a21 * y + 1;         // = gamma
292                if (Arithmetic.isZero(c)) {
293                        throw new ArithmeticException("getJacobian(): division by zero!");
294                }
295                final double cc = c * c;
296                return new double[][]
297                        {{x/c, y/c, 0,   0,   -(x*a)/cc, -(y*a)/cc, 1/c, 0  },
298                         {0,   0,   x/c, y/c, -(x*b)/cc, -(y*b)/cc, 0,   1/c}};
299        }
300        
301        // -----------------------------------------------------------------
302        
303        /**
304         * For testing only.
305         * @param args ignored
306         */
307        public static void main(String[] args) {
308                PrintPrecision.set(6);
309
310                // book example:
311                Pnt2d[] P = {
312                                PntDouble.from(2, 5),
313                                PntDouble.from(4, 6),
314                                PntDouble.from(7, 9),
315                                PntDouble.from(5, 9),
316                                PntDouble.from(5.2, 9.1)        // 5 points, overdetermined!
317                                };
318                
319                Pnt2d[] Q = {
320                                PntDouble.from(4, 3),
321                                PntDouble.from(5, 2),
322                                PntDouble.from(9, 3),
323                                PntDouble.from(7, 5),
324                                PntDouble.from(7, 4.9)  // 5 points, overdetermined!
325                                };
326                
327                ProjectiveMapping2D pm = ProjectiveMapping2D.fromPoints(P, Q);
328                
329                System.out.println("\nprojective mapping = \n" + pm.toString());
330                
331                for (int i = 0; i < P.length; i++) {
332                        Pnt2d Bi = pm.applyTo(P[i]);
333                        System.out.println(P[i].toString() + " -> " + Bi.toString());
334                }
335                
336                System.out.println("pm is of class " + pm.getClass().getName());
337                ProjectiveMapping2D pmi = pm.getInverse();
338                pmi = pmi.normalize();
339                System.out.println("\ninverse projective mapping (normalized) = \n" + pmi.toString());
340                
341                for (int i = 0; i < Q.length; i++) {
342                        Pnt2d Ai = pmi.applyTo(Q[i]);
343                        System.out.println(Q[i].toString() + " -> " + Ai.toString());
344                }
345                
346                ProjectiveMapping2D testId = pm.concat(pmi);
347                System.out.println("\ntest: should be a scaled identity matrix: = \n" + testId.toString());
348        }
349        /*
350         projective mapping = 
351        {{-0.732424, 1.208712, 0.244379}, 
352        {-1.702388, 1.725545, -1.654658}, 
353        {-0.206047, 0.123181, 1.000000}}
354        Point[2,000, 5,000] -> Point[4,007, 2,964]
355        Point[4,000, 6,000] -> Point[4,992, 2,065]
356        Point[7,000, 9,000] -> Point[8,999, 2,939]
357        Point[5,000, 9,000] -> Point[6,918, 4,973]
358        Point[5,200, 9,100] -> Point[7,084, 4,950]
359        pm is of class imagingbook.pub.geometry.mappings.linear.ProjectiveMapping2D
360        
361        inverse projective mapping (normalized) = 
362        {{2.430345, -1.484645, -3.050507}, 
363        {2.573894, -0.859176, -2.050649}, 
364        {0.183712, -0.200073, 1.000000}}
365        Point[4,000, 3,000] -> Point[1,954, 4,995]
366        Point[5,000, 2,000] -> Point[4,038, 5,993]
367        Point[9,000, 3,000] -> Point[6,998, 9,028]
368        Point[7,000, 5,000] -> Point[5,086, 9,078]
369        Point[7,000, 4,900] -> Point[5,122, 9,005]
370        
371        test: should be a scaled identity matrix: = 
372        {{1.000000, -0.000000, -0.000000}, 
373        {0.000000, 1.000000, -0.000000}, 
374        {0.000000, 0.000000, 1.000000}}
375         */
376}