Source code

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 *******************************************************************************/
009
010package imagingbook.pub.geometry.mappings.linear;
011
012import imagingbook.lib.settings.PrintPrecision;
013import imagingbook.pub.geometry.mappings.WarpParameters;
014
015import java.awt.geom.Point2D;
016
017import org.apache.commons.math3.linear.DecompositionSolver;
018import org.apache.commons.math3.linear.MatrixUtils;
019import org.apache.commons.math3.linear.RealMatrix;
020import org.apache.commons.math3.linear.RealVector;
021import org.apache.commons.math3.linear.SingularValueDecomposition;
022
023
024/**
025 * 2015-02-16: Added preliminary constructor for more than 4 point pairs 
026 * (overdetermined, least-squares).
027 * 
028 * @author WB
029 *
030 */
031public class ProjectiveMapping extends LinearMapping implements WarpParameters {
032        
033        
034        public static ProjectiveMapping makeMapping (Point2D[] P, Point2D[] Q) {
035                int minLen = Math.min(P.length, Q.length);
036                if (minLen < 1) {
037                        throw new IllegalArgumentException("cannot create a mapping from zero points");
038                }
039                else if (minLen <= 2) {                         // TODO: special case for minLen == 2?
040                        return new Translation(P, Q);
041                }
042                else if (minLen <= 3) {
043                        return new AffineMapping(P, Q);
044                }
045                else {
046                        return new ProjectiveMapping(P, Q);     
047                }
048        }
049        
050        // creates the identity mapping:
051        public ProjectiveMapping() {
052                super();
053        }
054        
055        public ProjectiveMapping(
056                        double a00, double a01, double a02, 
057                        double a10, double a11, double a12, 
058                        double a20, double a21, 
059                        boolean inv) {
060                super(a00, a01, a02, a10, a11, a12, a20, a21, 1, inv);
061        }
062        
063        public ProjectiveMapping(LinearMapping lm) {
064                super(lm);
065                //this.normalize();     // needed??
066        }
067        
068        // creates the projective mapping from the unit square S to
069        // the arbitrary quadrilateral P given by points p0,...,p3:
070        public ProjectiveMapping(Point2D p0, Point2D p1, Point2D p2, Point2D p3) {
071                super();
072                double x0 = p0.getX(), x1 = p1.getX(), x2 = p2.getX(), x3 = p3.getX(); 
073                double y0 = p0.getY(), y1 = p1.getY(), y2 = p2.getY(), y3 = p3.getY();
074                double S = (x1-x2)*(y3-y2) - (x3-x2)*(y1-y2);
075                // TODO: check S for zero value and throw exception
076                a20 = ((x0-x1+x2-x3)*(y3-y2)-(y0-y1+y2-y3)*(x3-x2)) / S;
077                a21 = ((y0-y1+y2-y3)*(x1-x2)-(x0-x1+x2-x3)*(y1-y2)) / S;
078                a00 = x1 - x0 + a20*x1;
079                a01 = x3 - x0 + a21*x3;
080                a02 = x0;
081                a10 = y1 - y0 + a20*y1;
082                a11 = y3 - y0 + a21*y3;
083                a12 = y0;
084        }
085        
086        
087        /**
088         * Creates a new {@link ProjectiveMapping} between arbitrary quadrilaterals P, Q.
089         * @param p0 point 1 of source quad P.
090         * @param p1 point 2 of source quad P.
091         * @param p2 point 3 of source quad P.
092         * @param p3 point 4 of source quad P.
093         * @param q0 point 1 of target quad Q.
094         * @param q1 point 2 of target quad Q.
095         * @param q2 point 3 of target quad Q.
096         * @param q3 point 4 of target quad Q.
097         */
098        public ProjectiveMapping(
099                        Point2D p0, Point2D p1, Point2D p2, Point2D p3, 
100                        Point2D q0, Point2D q1, Point2D q2, Point2D q3) {
101                super();        // initialized to identity
102                ProjectiveMapping T1 = new ProjectiveMapping(p0, p1, p2, p3);
103                ProjectiveMapping T2 = new ProjectiveMapping(q0, q1, q2, q3);
104                ProjectiveMapping T1i = T1.invert();
105                ProjectiveMapping T12 = T1i.concat(T2);         
106                this.concatDestructive(T12);    // transfer T12 -> this
107        }
108        
109        /**
110         * Creates a new {@link ProjectiveMapping} between arbitrary quadrilaterals P, Q.
111         * @param P source quad.
112         * @param Q target quad.
113         */
114        public ProjectiveMapping(Point2D[] P, Point2D[] Q) {
115                this(P[0], P[1], P[2], P[3], Q[0], Q[1], Q[2], Q[3]);
116        }
117        
118        /**
119         * Constructor for more than 4 point pairs, finds a least-squares solution
120         * for the homography parameters.
121         * NOTE: this is UNFINISHED code!
122         * @param P sequence of points (source)
123         * @param Q sequence of points (target)
124         * @param dummy unused (only to avoid duplicate signature)
125         */
126        public ProjectiveMapping(Point2D[] P, Point2D[] Q, boolean dummy) {
127                final int n = P.length;
128                double[] ba = new double[2 * n];
129                double[][] Ma = new double[2 * n][];
130                for (int i = 0; i < n; i++) {
131                        double x = P[i].getX();
132                        double y = P[i].getY();
133                        double u = Q[i].getX();
134                        double v = Q[i].getY();
135                        ba[2 * i + 0] = u;
136                        ba[2 * i + 1] = v;
137                        Ma[2 * i + 0] = new double[] { x, y, 1, 0, 0, 0, -u * x, -u * y };
138                        Ma[2 * i + 1] = new double[] { 0, 0, 0, x, y, 1, -v * x, -v * y };
139                }
140                
141                RealMatrix M = MatrixUtils.createRealMatrix(Ma);
142                RealVector b = MatrixUtils.createRealVector(ba);
143                DecompositionSolver solver = new SingularValueDecomposition(M).getSolver();
144                RealVector h = solver.solve(b);
145                a00 = h.getEntry(0);
146                a01 = h.getEntry(1);
147                a02 = h.getEntry(2);
148                a10 = h.getEntry(3);
149                a11 = h.getEntry(4);
150                a12 = h.getEntry(5);
151                a20 = h.getEntry(6);
152                a21 = h.getEntry(7);
153                a22 = 1;
154        }
155        
156        // -----------------------------------------------------------
157        
158        
159        public ProjectiveMapping concat(ProjectiveMapping B) {
160                ProjectiveMapping A = new ProjectiveMapping(this);
161                A.concatDestructive(B);
162                return A;
163        }
164        
165        public ProjectiveMapping invert() {
166                ProjectiveMapping pm = new ProjectiveMapping(this);
167                pm.invertDestructive();
168                return pm;
169        }
170        
171        void normalize() {
172                // scales the matrix such that a22 becomes 1
173                // TODO: check a22 for zero value and throw exception
174                a00 = a00/a22;          a01 = a01/a22;          a02 = a02/a22;
175                a10 = a10/a22;          a11 = a11/a22;          a12 = a12/a22;
176                a20 = a20/a22;          a21 = a21/a22;          a22 = 1;
177        }
178        
179        @Override
180        public ProjectiveMapping duplicate() {
181                return (ProjectiveMapping) this.clone();
182        }
183        
184        // Warp parameter support -------------------------------------
185        
186        public int getWarpParameterCount() {
187                return 8;
188        }
189        
190        public double[] getWarpParameters() {
191                double[] p = new double[] {
192                        a00 - 1, a01, 
193                        a10, a11 - 1, 
194                        a20, a21,
195                        a02, a12,
196                        };
197                return p;
198        }
199        
200//      p[0] = M3x3[0][0] - 1;  // = a
201//      p[1] = M3x3[0][1];              // = b
202//      p[2] = M3x3[1][0];              // = c
203//      p[3] = M3x3[1][1] - 1;  // = d
204//      p[4] = M3x3[2][0];              // = e
205//      p[5] = M3x3[2][1];              // = f
206//      p[6] = M3x3[0][2];              // = tx
207//      p[7] = M3x3[1][2];              // = ty
208
209
210        public void setWarpParameters(double[] p) {
211                a00 = p[0] + 1;   a01 = p[1];        a02 = p[6];
212                a10 = p[2];       a11 = p[3] + 1;    a12 = p[7];
213                a20 = p[4];       a21 = p[5];        a22 = 1;
214        }
215        
216
217        public double[][] getWarpJacobian(double[] xy) {
218                // see Baker 2003 "20 Years" Part 1, Eq. 99 (p. 46)
219                final double x = xy[0];
220                final double y = xy[1];
221                double a = a00 * x + a01 * y + a02;     // = alpha
222                double b = a10 * x + a11 * y + a12;     // = beta
223                double c = a20 * x + a21 * y + 1;       // = gamma
224                double cc = c * c;
225                // TODO: check c for zero-value and throw exception, make more efficient
226                return new double[][]
227                        {{x/c, y/c, 0,   0,   -(x*a)/cc, -(y*a)/cc, 1/c, 0  },
228                         {0,   0,   x/c, y/c, -(x*b)/cc, -(y*b)/cc, 0,   1/c}};
229        }
230        
231        // for testing only -----------------------------------------------------------------
232        
233        public static void main(String[] args) {
234                PrintPrecision.set(6);
235
236                // book example:
237                Point2D[] A = new Point2D[] {
238                                new Point2D.Double(2,5),
239                                new Point2D.Double(4,6),
240                                new Point2D.Double(7,9),
241                                new Point2D.Double(5,9),
242                                new Point2D.Double(5.2,9.1)     // 5 points, overdetermined!
243                                };
244                
245                Point2D[] B = new Point2D[] {
246                                new Point2D.Double(4,3),
247                                new Point2D.Double(5,2),
248                                new Point2D.Double(9,3),
249                                new Point2D.Double(7,5),
250                                new Point2D.Double(7,4.9)       // 5 points, overdetermined!
251                                };
252                
253                ProjectiveMapping pm = new ProjectiveMapping(A, B, true);
254                
255                System.out.println("\nprojective mapping = " + pm.toString());
256                
257                for (int i = 0; i < A.length; i++) {
258                        Point2D Bi = pm.applyTo(A[i]);
259                        System.out.println(A[i].toString() + " -> " + Bi.toString());
260                }
261                
262                System.out.println();
263                ProjectiveMapping pmi = pm.invert();
264                System.out.println("\ninverse projective mapping = " + pmi.toString());
265                for (int i = 0; i < B.length; i++) {
266                        Point2D Ai = pmi.applyTo(B[i]);
267                        System.out.println(B[i].toString() + " -> " + Ai.toString());
268                }
269        }
270        
271        
272}