View Javadoc
1   /*
2    * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *         http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  package com.irurueta.geometry;
17  
18  import com.irurueta.algebra.AlgebraException;
19  import com.irurueta.algebra.ArrayUtils;
20  import com.irurueta.algebra.LUDecomposer;
21  import com.irurueta.algebra.Matrix;
22  import com.irurueta.algebra.RQDecomposer;
23  import com.irurueta.algebra.SingularValueDecomposer;
24  import com.irurueta.algebra.Utils;
25  import com.irurueta.algebra.WrongSizeException;
26  
27  import java.io.Serializable;
28  import java.util.Arrays;
29  
30  /**
31   * This class performs projective transformations on 2D space.
32   * Projective transformations include any possible transformation that can be
33   * applied to 2D points.
34   */
35  @SuppressWarnings("DuplicatedCode")
36  public class ProjectiveTransformation2D extends Transformation2D implements Serializable {
37  
38      /**
39       * Constant indicating number of coordinates required in translation arrays.
40       */
41      public static final int NUM_TRANSLATION_COORDS = 2;
42  
43      /**
44       * Constant indicating the number of projective parameters that can be set
45       * in projective parameters array.
46       */
47      public static final int NUM_PROJECTIVE_PARAMS = 3;
48  
49      /**
50       * Constant defining number of inhomogeneous coordinates in 2D space.
51       */
52      public static final int INHOM_COORDS = 2;
53  
54      /**
55       * Constant defining number of homogeneous coordinates in 2D space.
56       */
57      public static final int HOM_COORDS = 3;
58  
59      /**
60       * Machine precision.
61       */
62      public static final double EPS = 1e-12;
63  
64      /**
65       * Constant defining a large threshold to consider a matrix valid as
66       * rotation.
67       */
68      private static final double LARGE_ROTATION_MATRIX_THRESHOLD = 1.0;
69  
70      /**
71       * Internal 3x3 matrix containing transformation.
72       */
73      private Matrix t;
74  
75      /**
76       * Indicates whether internal matrix is normalized.
77       */
78      private boolean normalized;
79  
80      /**
81       * Empty constructor.
82       * Creates transformation that has no effect.
83       */
84      public ProjectiveTransformation2D() {
85          super();
86          try {
87              t = Matrix.identity(HOM_COORDS, HOM_COORDS);
88          } catch (final WrongSizeException ignore) {
89              // never happens
90          }
91          normalize();
92      }
93  
94      /**
95       * Creates transformation with provided internal matrix.
96       * Notice that provided matrix should usually be invertible, otherwise the
97       * transformation will be degenerate and its inverse will not be available.
98       *
99       * @param t Internal 3x3 matrix.
100      * @throws NullPointerException     raised if provided matrix is null.
101      * @throws IllegalArgumentException raised if provided matrix is not 3x3
102      */
103     public ProjectiveTransformation2D(final Matrix t) {
104         setT(t);
105         normalize();
106     }
107 
108     /**
109      * Creates transformation with provided scale value.
110      *
111      * @param scale scale value. Values between 0.0 and 1.0 reduce objects,
112      *              values greater than 1.0 enlarge objects and negative values reverse
113      *              objects.
114      */
115     public ProjectiveTransformation2D(final double scale) {
116         final var diag = new double[HOM_COORDS];
117         Arrays.fill(diag, scale);
118         // set last element to 1.0
119         diag[HOM_COORDS - 1] = 1.0;
120         t = Matrix.diagonal(diag);
121         normalize();
122     }
123 
124     /**
125      * Creates transformation with provided rotation.
126      *
127      * @param rotation a 2D rotation.
128      * @throws NullPointerException raised if provided rotation is null.
129      */
130     public ProjectiveTransformation2D(final Rotation2D rotation) {
131         t = rotation.asHomogeneousMatrix();
132         normalize();
133     }
134 
135     /**
136      * Creates transformation with provided scale and rotation.
137      *
138      * @param scale    Scale value. Values between 0.0 and 1.0 reduce objects,
139      *                 values greater than 1.0 enlarge objects and negative values reverse
140      *                 objects.
141      * @param rotation a 2D rotation.
142      * @throws NullPointerException raised if provided rotation is null.
143      */
144     public ProjectiveTransformation2D(final double scale, final Rotation2D rotation) {
145         try {
146             final var diag = new double[INHOM_COORDS];
147             Arrays.fill(diag, scale);
148             final var a = Matrix.diagonal(diag);
149             a.multiply(rotation.asInhomogeneousMatrix());
150             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
151             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
152                     INHOM_COORDS - 1, a);
153         } catch (final WrongSizeException ignore) {
154             // never happens
155         }
156         normalize();
157     }
158 
159     /**
160      * Creates transformation with provided affine parameters and rotation.
161      *
162      * @param params   affine parameters including horizontal scaling, vertical
163      *                 scaling and skewness.
164      * @param rotation a 2D rotation.
165      * @throws NullPointerException raised if provided parameters are null or
166      *                              if provided rotation is null.
167      */
168     public ProjectiveTransformation2D(final AffineParameters2D params, final Rotation2D rotation) {
169         try {
170             final var a = params.asMatrix();
171             a.multiply(rotation.asInhomogeneousMatrix());
172             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
173             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
174                     INHOM_COORDS - 1, a);
175         } catch (final WrongSizeException ignore) {
176             // never happens
177         }
178         normalize();
179     }
180 
181     /**
182      * Creates transformation with provided 2D translation.
183      *
184      * @param translation array indicating 2D translation using inhomogeneous
185      *                    coordinates.
186      * @throws NullPointerException     raised if provided array is null.
187      * @throws IllegalArgumentException raised if length of array is not equal
188      *                                  to NUM_TRANSLATION_COORDS.
189      */
190     public ProjectiveTransformation2D(final double[] translation) {
191         if (translation.length != NUM_TRANSLATION_COORDS) {
192             throw new IllegalArgumentException();
193         }
194 
195         try {
196             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
197             t.setSubmatrix(0, 2, 1, 2, translation);
198         } catch (final WrongSizeException ignore) {
199             // never happens
200         }
201         normalize();
202     }
203 
204     /**
205      * Creates transformation with provided affine linear mapping and
206      * translation.
207      *
208      * @param a           affine linear mapping.
209      * @param translation array indicating 2D translation using inhomogeneous
210      *                    coordinates.
211      * @throws NullPointerException     raised if provided array is null or if
212      *                                  affine linear mapping is null.
213      * @throws IllegalArgumentException raised if length of array is not equal
214      *                                  to NUM_TRANSLATION_COORDS.
215      */
216     public ProjectiveTransformation2D(final Matrix a, final double[] translation) {
217         if (translation.length != NUM_TRANSLATION_COORDS) {
218             throw new IllegalArgumentException();
219         }
220 
221         try {
222             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
223             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
224                     INHOM_COORDS - 1, a);
225             t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
226                     HOM_COORDS - 1, translation);
227         } catch (final WrongSizeException ignore) {
228             // never happens
229         }
230         normalize();
231     }
232 
233     /**
234      * Creates transformation with provided scale and translation.
235      *
236      * @param scale       scale value. Values between 0.0 and 1.0 reduce objects,
237      *                    values greater than 1.0 enlarge objects and negative values reverse
238      *                    objects.
239      * @param translation array indicating 2D translation using inhomogeneous
240      *                    coordinates.
241      * @throws NullPointerException     raised if provided translation is null
242      * @throws IllegalArgumentException raised if provided translation does not
243      *                                  have length 2.
244      */
245     public ProjectiveTransformation2D(final double scale, final double[] translation) {
246         if (translation.length != NUM_TRANSLATION_COORDS) {
247             throw new IllegalArgumentException();
248         }
249 
250         final var diag = new double[HOM_COORDS];
251         Arrays.fill(diag, scale);
252         // set last element to 1.0
253         diag[HOM_COORDS - 1] = 1.0;
254         t = Matrix.diagonal(diag);
255 
256         // set translation
257         t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
258                 HOM_COORDS - 1, translation);
259         normalize();
260     }
261 
262     /**
263      * Creates transformation with provided rotation and translation.
264      *
265      * @param rotation    a 2D rotation.
266      * @param translation array indicating 2D translation using inhomogeneous
267      *                    coordinates.
268      * @throws NullPointerException     raised if provided rotation or translation
269      *                                  is null.
270      * @throws IllegalArgumentException raised if provided translation does not
271      *                                  have length 2.
272      */
273     public ProjectiveTransformation2D(final Rotation2D rotation, final double[] translation) {
274         if (translation.length != NUM_TRANSLATION_COORDS) {
275             throw new IllegalArgumentException();
276         }
277 
278         t = rotation.asHomogeneousMatrix();
279 
280         // set translation
281         t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
282                 HOM_COORDS - 1, translation);
283         normalize();
284     }
285 
286     /**
287      * Creates transformation with provided scale, rotation and translation.
288      *
289      * @param scale       Scale value. Values between 0.0 and 1.0 reduce objects,
290      *                    values greater than 1.0 enlarge objects and negative values reverse
291      *                    objects.
292      * @param rotation    a 2D rotation.
293      * @param translation array indicating 2D translation using inhomogeneous
294      *                    coordinates.
295      * @throws NullPointerException     raised if provided rotation or translation
296      *                                  is null.
297      * @throws IllegalArgumentException raised if provided translation does not
298      *                                  have length 2.
299      */
300     public ProjectiveTransformation2D(final double scale, final Rotation2D rotation, final double[] translation) {
301         if (translation.length != NUM_TRANSLATION_COORDS) {
302             throw new IllegalArgumentException();
303         }
304 
305         try {
306             final var diag = new double[INHOM_COORDS];
307             Arrays.fill(diag, scale);
308             final var a = Matrix.diagonal(diag);
309             a.multiply(rotation.asInhomogeneousMatrix());
310 
311             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
312             // set A
313             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
314                     INHOM_COORDS - 1, a);
315             // set translation
316             t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
317                     HOM_COORDS - 1, translation);
318         } catch (final WrongSizeException ignore) {
319             // never happens
320         }
321         normalize();
322     }
323 
324     /**
325      * Creates transformation with provided scale, rotation and translation.
326      *
327      * @param scale                scale value. Values between 0.0 and 1.0 reduce objects,
328      *                             values greater than 1.0 enlarge objects and negative values reverse
329      *                             objects.
330      * @param rotation             a 2D rotation.
331      * @param translation          array indicating 2D translation using inhomogeneous
332      *                             coordinates.
333      * @param projectiveParameters array of length 3 containing projective
334      *                             parameters.
335      * @throws NullPointerException     raised if provided rotation or translation
336      *                                  is null.
337      * @throws IllegalArgumentException raised if provided translation does not
338      *                                  have length 2 or if projective parameters array doesn't have length 3.
339      */
340     public ProjectiveTransformation2D(final double scale, final Rotation2D rotation, final double[] translation,
341                                       final double[] projectiveParameters) {
342         if (translation.length != NUM_TRANSLATION_COORDS) {
343             throw new IllegalArgumentException();
344         }
345         if (projectiveParameters.length != HOM_COORDS) {
346             throw new IllegalArgumentException();
347         }
348 
349         try {
350             final var value = projectiveParameters[HOM_COORDS - 1];
351             final var diag = new double[INHOM_COORDS];
352             Arrays.fill(diag, scale);
353             final var a = Matrix.diagonal(diag);
354             a.multiply(rotation.asInhomogeneousMatrix());
355 
356             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
357             // set A
358             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
359                     INHOM_COORDS - 1, a);
360             // set translation
361             t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
362                     HOM_COORDS - 1, translation);
363             t.multiplyByScalar(value);
364 
365             t.setSubmatrix(HOM_COORDS - 1, 0, HOM_COORDS - 1,
366                     HOM_COORDS - 1, projectiveParameters);
367         } catch (final WrongSizeException ignore) {
368             // never happens
369         }
370         normalize();
371     }
372 
373     /**
374      * Creates transformation with provided parameters, rotation and
375      * translation.
376      *
377      * @param params      Affine parameters including horizontal scaling, vertical
378      *                    scaling and skewness.
379      * @param rotation    a 2D rotation.
380      * @param translation array indicating 2D translation using inhomogeneous
381      *                    coordinates.
382      * @throws NullPointerException     raised if provided parameters, rotation or
383      *                                  translation is null.
384      * @throws IllegalArgumentException raised if provided translation does not
385      *                                  have length 2.
386      */
387     public ProjectiveTransformation2D(final AffineParameters2D params, final Rotation2D rotation,
388                                       final double[] translation) {
389         if (translation.length != NUM_TRANSLATION_COORDS) {
390             throw new IllegalArgumentException();
391         }
392 
393         try {
394             final var a = params.asMatrix();
395             a.multiply(rotation.asInhomogeneousMatrix());
396             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
397             // set A
398             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
399                     INHOM_COORDS - 1, a);
400             // set translation
401             t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
402                     HOM_COORDS - 1, translation);
403         } catch (final WrongSizeException ignore) {
404             // never happens
405         }
406         normalize();
407     }
408 
409     /**
410      * Creates transformation with provided parameters, rotation and
411      * translation.
412      *
413      * @param params               affine parameters including horizontal scaling, vertical
414      *                             scaling and skewness.
415      * @param rotation             a 2D rotation.
416      * @param translation          array indicating 2D translation using inhomogeneous
417      *                             coordinates.
418      * @param projectiveParameters array of length 3 containing projective
419      *                             parameters.
420      * @throws NullPointerException     raised if provided parameters, rotation or
421      *                                  translation is null.
422      * @throws IllegalArgumentException raised if provided translation does not
423      *                                  have length 2 or if projective parameters array doesn't have length 3.
424      */
425     public ProjectiveTransformation2D(final AffineParameters2D params, final Rotation2D rotation,
426                                       final double[] translation, final double[] projectiveParameters) {
427         if (translation.length != NUM_TRANSLATION_COORDS) {
428             throw new IllegalArgumentException();
429         }
430         if (projectiveParameters.length != HOM_COORDS) {
431             throw new IllegalArgumentException();
432         }
433 
434         try {
435             final var a = params.asMatrix();
436             a.multiply(rotation.asInhomogeneousMatrix());
437             t = Matrix.identity(HOM_COORDS, HOM_COORDS);
438             // set A
439             t.setSubmatrix(0, 0, INHOM_COORDS - 1,
440                     INHOM_COORDS - 1, a);
441             // set translation
442             t.setSubmatrix(0, HOM_COORDS - 1, translation.length - 1,
443                     HOM_COORDS - 1, translation);
444             final var value = projectiveParameters[HOM_COORDS - 1];
445             t.multiplyByScalar(value);
446 
447             t.setSubmatrix(HOM_COORDS - 1, 0, HOM_COORDS - 1,
448                     HOM_COORDS - 1, projectiveParameters);
449         } catch (final WrongSizeException ignore) {
450             // never happens
451         }
452         normalize();
453     }
454 
455     /**
456      * Creates transformation by estimating its internal matrix by providing 4
457      * corresponding original and transformed points.
458      *
459      * @param inputPoint1  1st input point.
460      * @param inputPoint2  2nd input point.
461      * @param inputPoint3  3rd input point.
462      * @param inputPoint4  4th input point.
463      * @param outputPoint1 1st transformed point corresponding to 1st input
464      *                     point.
465      * @param outputPoint2 2nd transformed point corresponding to 2nd input
466      *                     point.
467      * @param outputPoint3 3rd transformed point corresponding to 3rd input
468      *                     point.
469      * @param outputPoint4 4th transformed point corresponding to 4th input
470      *                     point.
471      * @throws CoincidentPointsException raised if transformation cannot be
472      *                                   estimated for some reason (point configuration degeneracy, duplicate
473      *                                   points or numerical instabilities).
474      */
475     public ProjectiveTransformation2D(
476             final Point2D inputPoint1, final Point2D inputPoint2, final Point2D inputPoint3, final Point2D inputPoint4,
477             final Point2D outputPoint1, final Point2D outputPoint2, final Point2D outputPoint3,
478             final Point2D outputPoint4) throws CoincidentPointsException {
479         try {
480             t = new Matrix(HOM_COORDS, HOM_COORDS);
481         } catch (final WrongSizeException ignore) {
482             // never happens
483         }
484         setTransformationFromPoints(inputPoint1, inputPoint2, inputPoint3, inputPoint4, outputPoint1, outputPoint2,
485                 outputPoint3, outputPoint4);
486     }
487 
488     /**
489      * Creates transformation by estimating its internal matrix by providing 4
490      * corresponding original and transformed lines.
491      *
492      * @param inputLine1  1st input line.
493      * @param inputLine2  2nd input line.
494      * @param inputLine3  3rd input line.
495      * @param inputLine4  4th input line.
496      * @param outputLine1 1st transformed line corresponding to 1st input line.
497      * @param outputLine2 2nd transformed line corresponding to 2nd input line.
498      * @param outputLine3 3rd transformed line corresponding to 3rd input line.
499      * @param outputLine4 4th transformed line corresponding to 4th input line.
500      * @throws CoincidentLinesException raised if transformation cannot be
501      *                                  estimated for some reason (line configuration degeneracy, duplicate lines
502      *                                  or numerical instabilities).
503      */
504     public ProjectiveTransformation2D(
505             final Line2D inputLine1, final Line2D inputLine2, final Line2D inputLine3, final Line2D inputLine4,
506             final Line2D outputLine1, final Line2D outputLine2, final Line2D outputLine3, final Line2D outputLine4)
507             throws CoincidentLinesException {
508         setTransformationFromLines(inputLine1, inputLine2, inputLine3, inputLine4, outputLine1, outputLine2,
509                 outputLine3, outputLine4);
510     }
511 
512     /**
513      * Returns internal matrix containing this transformation data.
514      * Point transformation is computed as t * x, where x is a 2D point
515      * expressed using homogeneous coordinates.
516      * Usually the internal transformation matrix will be invertible.
517      * When this is not the case, the transformation is considered degenerate
518      * and its inverse will not be available.
519      *
520      * @return internal transformation matrix.
521      */
522     public Matrix getT() {
523         return t;
524     }
525 
526     /**
527      * Sets internal matrix containing this transformation data.
528      * Point transformation is computed as t * x, where x is a 2D point
529      * expressed using homogeneous coordinates.
530      * Usually provided matrix will be invertible, when this is not the case
531      * this transformation will become degenerate and its inverse will not be
532      * available.
533      * This method does not check whether provided matrix is invertible or not.
534      *
535      * @param t transformation matrix.
536      * @throws NullPointerException     raised if provided matrix is null.
537      * @throws IllegalArgumentException raised if provided matrix is not 3x3
538      */
539     public final void setT(final Matrix t) {
540         if (t.getRows() != HOM_COORDS || t.getColumns() != HOM_COORDS) {
541             throw new IllegalArgumentException();
542         }
543 
544         this.t = t;
545         normalized = false;
546     }
547 
548     /**
549      * Returns boolean indicating whether provided matrix will produce a
550      * degenerate projective transformation or not.
551      *
552      * @param t a 3x3 matrix to be used as the internal matrix of a projective
553      *          transformation.
554      * @return true if matrix will produce a degenerate transformation, false
555      * otherwise.
556      * @throws IllegalArgumentException raised if provided matrix is not 3x3.
557      */
558     public static boolean isDegenerate(final Matrix t) {
559         if (t.getRows() != HOM_COORDS || t.getColumns() != HOM_COORDS) {
560             throw new IllegalArgumentException();
561         }
562 
563         try {
564             final var decomposer = new LUDecomposer(t);
565             decomposer.decompose();
566             return decomposer.isSingular();
567         } catch (final AlgebraException e) {
568             // if decomposition fails, assume that matrix is degenerate because
569             // of numerical instabilities
570             return true;
571         }
572     }
573 
574     /**
575      * Indicates whether this transformation is degenerate.
576      * When a transformation is degenerate, its inverse cannot be computed.
577      *
578      * @return true if transformation is degenerate, false otherwise.
579      */
580     public boolean isDegenerate() {
581         return isDegenerate(t);
582     }
583 
584     /**
585      * Returns affine linear mapping matrix.
586      *
587      * @return linear mapping matrix.
588      * @see AffineTransformation2D
589      */
590     public Matrix getA() {
591         final var a = t.getSubmatrix(0, 0,
592                 INHOM_COORDS - 1, INHOM_COORDS - 1);
593         a.multiplyByScalar(1.0 / t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1));
594         return a;
595     }
596 
597     /**
598      * Sets affine linear mapping matrix.
599      *
600      * @param a linear mapping matrix.
601      * @throws NullPointerException     raised if provided matrix is null.
602      * @throws IllegalArgumentException raised if provided matrix does not have
603      *                                  size 2x2.
604      * @see AffineTransformation2D
605      */
606     public final void setA(final Matrix a) {
607         if (a == null) {
608             throw new NullPointerException();
609         }
610         if (a.getRows() != INHOM_COORDS || a.getColumns() != INHOM_COORDS) {
611             throw new IllegalArgumentException();
612         }
613 
614         t.setSubmatrix(0, 0, INHOM_COORDS - 1, INHOM_COORDS - 1,
615                 a.multiplyByScalarAndReturnNew(t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1)));
616         normalized = false;
617     }
618 
619     /**
620      * Normalizes current matrix instance.
621      */
622     public final void normalize() {
623         if (!normalized) {
624             final var norm = Utils.normF(t);
625             if (norm > EPS) {
626                 t.multiplyByScalar(1.0 / norm);
627             }
628             normalized = true;
629         }
630     }
631 
632     /**
633      * Returns the 2D rotation component associated to this transformation.
634      * Note: if this rotation instance is modified, its changes won't be
635      * reflected on this transformation until rotation is set again.
636      *
637      * @return 2D rotation.
638      * @throws AlgebraException if for some reason rotation cannot be estimated
639      *                          (usually because of numerical instability).
640      */
641     public Rotation2D getRotation() throws AlgebraException {
642         // Use QR decomposition to retrieve rotation component of this
643         // transformation
644         normalize();
645         final var decomposer = new RQDecomposer(t.getSubmatrix(0, 0,
646                 INHOM_COORDS - 1, INHOM_COORDS - 1));
647         try {
648             decomposer.decompose();
649             return new Rotation2D(decomposer.getQ(), LARGE_ROTATION_MATRIX_THRESHOLD); //a large threshold is
650             // used because Q matrix is always assumed to be orthonormal
651         } catch (final InvalidRotationMatrixException ignore) {
652             return null;
653         }
654     }
655 
656     /**
657      * Sets 2D rotation for this transformation.
658      *
659      * @param rotation a 2D rotation.
660      * @throws NullPointerException raised if provided rotation is null.
661      * @throws AlgebraException     raised if for numerical reasons rotation cannot
662      *                              be set (usually because of numerical instability in parameters of this
663      *                              transformation).
664      */
665     public void setRotation(final Rotation2D rotation) throws AlgebraException {
666         final var rotMatrix = rotation.asInhomogeneousMatrix();
667 
668         // Use QR decomposition to retrieve parameters matrix
669         final var decomposer = new RQDecomposer(t.getSubmatrix(0, 0,
670                 INHOM_COORDS - 1, INHOM_COORDS - 1));
671         decomposer.decompose();
672         // retrieves params matrix
673         final var localA = decomposer.getR();
674         localA.multiply(rotMatrix);
675         t.setSubmatrix(0, 0, INHOM_COORDS - 1, INHOM_COORDS - 1,
676                 localA);
677         normalized = false;
678     }
679 
680     /**
681      * Adds provided rotation to current rotation assigned to this
682      * transformation.
683      *
684      * @param rotation 2D rotation to be added.
685      * @throws AlgebraException raised if for numerical reasons rotation cannot
686      *                          be set (usually because of numerical instability in parameters of this
687      *                          transformation).
688      */
689     public void addRotation(final Rotation2D rotation) throws AlgebraException {
690         final var localRotation = getRotation();
691         localRotation.combine(rotation);
692         setRotation(localRotation);
693     }
694 
695     /**
696      * Sets scale of this transformation.
697      *
698      * @param scale scale value to be set. A value between 0.0 and 1.0 indicates
699      *              that objects will be reduced, a value greater than 1.0 indicates that
700      *              objects will be enlarged, and a negative value indicates that objects
701      *              will be reversed.
702      * @throws AlgebraException raised if for numerical reasons scale cannot
703      *                          be set (usually because of numerical instability in parameters of this
704      *                          transformation).
705      */
706     public void setScale(final double scale) throws AlgebraException {
707         normalize();
708         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
709         final var decomposer = new RQDecomposer(t.getSubmatrix(0, 0,
710                 INHOM_COORDS - 1, INHOM_COORDS - 1));
711         decomposer.decompose();
712         final var localA = decomposer.getR(); //params
713         localA.setElementAt(0, 0, scale * value);
714         localA.setElementAt(1, 1, scale * value);
715         localA.multiply(decomposer.getQ());
716         t.setSubmatrix(0, 0, INHOM_COORDS - 1, INHOM_COORDS - 1,
717                 localA);
718         normalized = false;
719     }
720 
721     /**
722      * Gets affine parameters of associated to this instance.
723      * Affine parameters contain horizontal scale, vertical scale and skewness
724      * of axes.
725      *
726      * @return affine parameters.
727      * @throws AlgebraException raised if for numerical reasons affine
728      *                          parameters cannot be retrieved (usually because of numerical instability
729      *                          of the internal matrix of this instance).
730      */
731     public AffineParameters2D getAffineParameters() throws AlgebraException {
732         final var parameters = new AffineParameters2D();
733         getAffineParameters(parameters);
734         return parameters;
735     }
736 
737     /**
738      * Computes affine parameters associated to this instance and stores the
739      * result in provided instance.
740      * Affine parameters contain horizontal scale, vertical scale and skewness
741      * of axes.
742      *
743      * @param result instance where affine parameters will be stored.
744      * @throws AlgebraException raised if for numerical reasons affine
745      *                          parameters cannot be retrieved (usually because of numerical instability
746      *                          of the internal matrix of this instance).
747      */
748     public void getAffineParameters(final AffineParameters2D result) throws AlgebraException {
749         normalize();
750         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
751         final var decomposer = new RQDecomposer(t.getSubmatrix(0, 0,
752                 INHOM_COORDS - 1, INHOM_COORDS - 1));
753         decomposer.decompose();
754         final var r = decomposer.getR();
755         r.multiplyByScalar(1.0 / value);
756         result.fromMatrix(r);
757     }
758 
759     /**
760      * Sets affine parameters associated to this instance.
761      * Affine parameters contain horizontal scale, vertical scale and skewness
762      * of axes.
763      *
764      * @param parameters affine parameters to be set.
765      * @throws AlgebraException raised if for numerical reasons affine
766      *                          parameters cannot be set (usually because of numerical instability of
767      *                          the internal matrix of this instance).
768      */
769     public void setAffineParameters(final AffineParameters2D parameters) throws AlgebraException {
770         normalize();
771         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
772         final var decomposer = new RQDecomposer(t.getSubmatrix(0, 0,
773                 INHOM_COORDS - 1, INHOM_COORDS - 1));
774         decomposer.decompose();
775         final var params = parameters.asMatrix();
776         final var rotation = decomposer.getQ();
777 
778         // params is equivalent to A because it
779         // has been multiplied by rotation
780         params.multiply(rotation);
781         // normalize
782         params.multiplyByScalar(value);
783         t.setSubmatrix(0, 0, INHOM_COORDS - 1, INHOM_COORDS - 1,
784                 params);
785         normalized = false;
786     }
787 
788     /**
789      * Returns the projective parameters associated to this instance.
790      * These parameters are the located in the last row of the internal
791      * transformation matrix.
792      * For affine, metric or Euclidean transformations this last row is always
793      * [0, 0, 1] (taking into account that transformation matrix is defined
794      * up to scale).
795      *
796      * @return Projective parameters returned as the array containing the values
797      * of the last row of the internal transformation matrix.
798      */
799     public double[] getProjectiveParameters() {
800         // return last row of matrix t
801         return t.getSubmatrixAsArray(HOM_COORDS - 1, 0, HOM_COORDS - 1,
802                 HOM_COORDS - 1, true);
803     }
804 
805     /**
806      * Sets the projective parameters associated to this instance.
807      * These parameters will be set in the last row of the internal
808      * transformation matrix.
809      * For affine, matrix or Euclidean transformations parameters are always
810      * [0, 0, 1] (taking into account that transformation matrix is defined up
811      * to scale).
812      *
813      * @param params projective parameters to be set. It must be an array of
814      *               length 3.
815      * @throws IllegalArgumentException raised if provided array does not have
816      *                                  length 3.
817      */
818     public final void setProjectiveParameters(final double[] params) {
819         if (params.length != HOM_COORDS) {
820             throw new IllegalArgumentException();
821         }
822 
823         t.setSubmatrix(HOM_COORDS - 1, 0, HOM_COORDS - 1,
824                 HOM_COORDS - 1, params);
825         normalized = false;
826     }
827 
828     /**
829      * Returns 2D translation assigned to this transformation as an array
830      * expressed in inhomogeneous coordinates.
831      * Note: Updating the values of the returned array will not update the
832      * translation of this instance. To do so, translation needs to be set
833      * again.
834      *
835      * @return 2D translation array.
836      */
837     public double[] getTranslation() {
838         normalize();
839         final var translation = t.getSubmatrixAsArray(0, HOM_COORDS - 1,
840                 INHOM_COORDS - 1, HOM_COORDS - 1);
841         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
842         ArrayUtils.multiplyByScalar(translation, 1.0 / value, translation);
843         return translation;
844     }
845 
846     /**
847      * Obtains 2D translation assigned to this transformation and stores result
848      * into provided array.
849      * Note: updating the values of the returned array will not update the
850      * translation of this instance. To do so, translation needs to be set.
851      *
852      * @param out array where translation values will be stored.
853      * @throws WrongSizeException if provided array does not have length 2.
854      */
855     public void getTranslation(final double[] out) throws WrongSizeException {
856         t.getSubmatrixAsArray(0, HOM_COORDS - 1,
857                 INHOM_COORDS - 1, HOM_COORDS - 1, out);
858         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
859         ArrayUtils.multiplyByScalar(out, 1.0 / value, out);
860     }
861 
862     /**
863      * Sets 2D translation assigned to this transformation as an array expressed
864      * in inhomogeneous coordinates.
865      *
866      * @param translation 2D translation array.
867      * @throws IllegalArgumentException raised if provided array does not have
868      *                                  length equal to NUM_TRANSLATION_COORDS.
869      */
870     public void setTranslation(final double[] translation) {
871         if (translation.length != NUM_TRANSLATION_COORDS) {
872             throw new IllegalArgumentException();
873         }
874 
875         final var value = t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
876         final var translation2 = ArrayUtils.multiplyByScalarAndReturnNew(translation, value);
877         t.setSubmatrix(0, HOM_COORDS - 1, translation2.length - 1,
878                 HOM_COORDS - 1, translation2);
879         normalized = false;
880     }
881 
882     /**
883      * Adds provided translation to current translation on this transformation.
884      * Provided translation must be expressed as an array of inhomogeneous
885      * coordinates.
886      *
887      * @param translation 2D translation array.
888      * @throws IllegalArgumentException raised if provided array does not have
889      *                                  length equal to NUM_TRANSLATION_COORDS.
890      */
891     public void addTranslation(final double[] translation) {
892         final var currentTranslation = getTranslation();
893         ArrayUtils.sum(currentTranslation, translation, currentTranslation);
894         setTranslation(currentTranslation);
895     }
896 
897     /**
898      * Returns current x coordinate translation assigned to this transformation.
899      *
900      * @return X coordinate translation.
901      */
902     public double getTranslationX() {
903         normalize();
904         return t.getElementAt(0, HOM_COORDS - 1)
905                 / t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
906     }
907 
908     /**
909      * Sets x coordinate translation to be made by this transformation.
910      *
911      * @param translationX X coordinate translation to be set.
912      */
913     public void setTranslationX(final double translationX) {
914         t.setElementAt(0, HOM_COORDS - 1,
915                 translationX * t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1));
916         normalized = false;
917     }
918 
919     /**
920      * Returns current y coordinate translation assigned to this transformation.
921      *
922      * @return Y coordinate translation.
923      */
924     public double getTranslationY() {
925         normalize();
926         return t.getElementAt(1, HOM_COORDS - 1)
927                 / t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1);
928     }
929 
930     /**
931      * Sets y coordinate translation to be made by this transformation.
932      *
933      * @param translationY Y coordinate translation to be set.
934      */
935     public void setTranslationY(final double translationY) {
936         t.setElementAt(1, HOM_COORDS - 1,
937                 translationY * t.getElementAt(HOM_COORDS - 1, HOM_COORDS - 1));
938         normalized = false;
939     }
940 
941     /**
942      * Sets x, y coordinates of translation to be made by this transformation.
943      *
944      * @param translationX translation x coordinate to be set.
945      * @param translationY translation y coordinate to be set.
946      */
947     public void setTranslation(final double translationX, final double translationY) {
948         setTranslationX(translationX);
949         setTranslationY(translationY);
950     }
951 
952     /**
953      * Sets x, y coordinates of translation to be made by this transformation.
954      *
955      * @param translation translation to be set.
956      */
957     public void setTranslation(final Point2D translation) {
958         setTranslation(translation.getInhomX(), translation.getInhomY());
959     }
960 
961     /**
962      * Gets x, y coordinates of translation to be made by this transformation
963      * as a new point.
964      *
965      * @return a new point containing translation coordinates.
966      */
967     public Point2D getTranslationPoint() {
968         final var out = Point2D.create();
969         getTranslationPoint(out);
970         return out;
971     }
972 
973     /**
974      * Gets x, y coordinates of translation to be made by this transformation
975      * and stores them into provided point.
976      *
977      * @param out point where translation coordinates will be stored.
978      */
979     public void getTranslationPoint(final Point2D out) {
980         out.setInhomogeneousCoordinates(getTranslationX(), getTranslationY());
981     }
982 
983     /**
984      * Adds provided x coordinate to current translation assigned to this
985      * transformation.
986      *
987      * @param translationX X coordinate to be added to current translation.
988      */
989     public void addTranslationX(final double translationX) {
990         setTranslationX(getTranslationX() + translationX);
991     }
992 
993     /**
994      * Adds provided y coordinate to current translation assigned to this
995      * transformation.
996      *
997      * @param translationY Y coordinate to be added to current translation.
998      */
999     public void addTranslationY(final double translationY) {
1000         setTranslationY(getTranslationY() + translationY);
1001     }
1002 
1003     /**
1004      * Adds provided coordinates to current translation assigned ot this
1005      * transformation.
1006      *
1007      * @param translationX x coordinate to be added to current translation.
1008      * @param translationY y coordinate to be added to current translation.
1009      */
1010     public void addTranslation(final double translationX, final double translationY) {
1011         addTranslationX(translationX);
1012         addTranslationY(translationY);
1013     }
1014 
1015     /**
1016      * Adds provided coordinates to current translation assigned to this
1017      * transformation.
1018      *
1019      * @param translation x, y coordinates to be added to current translation.
1020      */
1021     public void addTranslation(final Point2D translation) {
1022         addTranslation(translation.getInhomX(), translation.getInhomY());
1023     }
1024 
1025     /**
1026      * Represents this transformation as a 3x3 matrix.
1027      * A point can be transformed as t * p, where t is the transformation matrix
1028      * and p is a point expressed as an homogeneous vector.
1029      *
1030      * @return This transformation in matrix form.
1031      */
1032     @Override
1033     public Matrix asMatrix() {
1034         return new Matrix(t);
1035     }
1036 
1037     /**
1038      * Represents this transformation as a 3x3 matrix and stores the result in
1039      * provided instance.
1040      *
1041      * @param m Instance where transformation matrix will be stored.
1042      * @throws IllegalArgumentException Raised if provided instance is not a 3x3
1043      *                                  matrix.
1044      */
1045     @Override
1046     public void asMatrix(final Matrix m) {
1047         if (m.getRows() != HOM_COORDS || m.getColumns() != HOM_COORDS) {
1048             throw new IllegalArgumentException();
1049         }
1050 
1051         m.copyFrom(t);
1052     }
1053 
1054     /**
1055      * Transforms input point using this transformation and stores the result in
1056      * provided output points.
1057      *
1058      * @param inputPoint  point to be transformed.
1059      * @param outputPoint instance where transformed point data will be stored.
1060      */
1061     @Override
1062     public void transform(final Point2D inputPoint, final Point2D outputPoint) {
1063 
1064         inputPoint.normalize();
1065         normalize();
1066         try {
1067             final var point = new Matrix(Point2D.POINT2D_HOMOGENEOUS_COORDINATES_LENGTH, 1);
1068             point.setElementAtIndex(0, inputPoint.getHomX());
1069             point.setElementAtIndex(1, inputPoint.getHomY());
1070             point.setElementAtIndex(2, inputPoint.getHomW());
1071 
1072             final var transformedPoint = t.multiplyAndReturnNew(point);
1073 
1074             outputPoint.setHomogeneousCoordinates(
1075                     transformedPoint.getElementAtIndex(0),
1076                     transformedPoint.getElementAtIndex(1),
1077                     transformedPoint.getElementAtIndex(2));
1078         } catch (final WrongSizeException ignore) {
1079             // never happens
1080         }
1081     }
1082 
1083     /**
1084      * Transforms a conic using this transformation and stores the result into
1085      * provided output conic.
1086      *
1087      * @param inputConic  conic to be transformed.
1088      * @param outputConic instance where data of transformed conic will be
1089      *                    stored.
1090      * @throws NonSymmetricMatrixException raised if due to numerical precision
1091      *                                     the resulting output conic matrix is not considered to be symmetric.
1092      * @throws AlgebraException            raised if transform cannot be computed because of
1093      *                                     numerical instabilities.
1094      */
1095     @Override
1096     public void transform(final Conic inputConic, final Conic outputConic) throws NonSymmetricMatrixException,
1097             AlgebraException {
1098         // point' * conic * point = 0
1099         // point' * t' * transformedConic * t * point = 0
1100         // where:
1101         // - transformedPoint = t * point
1102 
1103         // Hence:
1104         // transformedConic = t^-' * conic * t^-1
1105 
1106         inputConic.normalize();
1107 
1108         final var c = inputConic.asMatrix();
1109         normalize();
1110 
1111         final var invT = inverseAndReturnNew().asMatrix();
1112         // normalize transformation matrix invT to increase accuracy
1113         var norm = Utils.normF(invT);
1114         invT.multiplyByScalar(1.0 / norm);
1115 
1116         final var m = invT.transposeAndReturnNew();
1117         try {
1118             m.multiply(c);
1119             m.multiply(invT);
1120         } catch (final WrongSizeException ignore) {
1121             // never happens
1122         }
1123 
1124         // normalize resulting m matrix to increase accuracy so that it can be
1125         // considered symmetric
1126         norm = Utils.normF(m);
1127         m.multiplyByScalar(1.0 / norm);
1128 
1129         outputConic.setParameters(m);
1130     }
1131 
1132     /**
1133      * Transforms a dual conic using this transformation and stores the result
1134      * into provided output dual conic.
1135      *
1136      * @param inputDualConic  Dual conic to be transformed.
1137      * @param outputDualConic Instance where data of transformed dual conic will
1138      *                        be stored.
1139      * @throws NonSymmetricMatrixException raised if due to numerical precision
1140      *                                     the resulting output dual conic matrix is not considered to be symmetric.
1141      * @throws AlgebraException            Raised if transform cannot be computed because
1142      *                                     of numerical instabilities.
1143      */
1144     @Override
1145     public void transform(final DualConic inputDualConic, final DualConic outputDualConic)
1146             throws NonSymmetricMatrixException, AlgebraException {
1147         // line' * dualConic * line = 0
1148         // line' * t^-1 * t * dualConic * t' * t^-1' * line
1149 
1150         // Hence:
1151         // transformed line : t^-1'*line
1152         // transformed dual conic: t * dualConic * t'
1153 
1154         inputDualConic.normalize();
1155         normalize();
1156 
1157         final var dualC = inputDualConic.asMatrix();
1158         final var transT = t.transposeAndReturnNew();
1159 
1160         final var m = t.multiplyAndReturnNew(dualC);
1161         m.multiply(transT);
1162 
1163         // normalize resulting m matrix to increase accuracy so that it can be
1164         // considered symmetric
1165         final var norm = Utils.normF(m);
1166         m.multiplyByScalar(1.0 / norm);
1167 
1168         outputDualConic.setParameters(m);
1169     }
1170 
1171     /**
1172      * Transforms provided input line using this transformation and stores the
1173      * result into provided output line instance.
1174      *
1175      * @param inputLine  line to be transformed.
1176      * @param outputLine instance where data of transformed line will be stored
1177      * @throws AlgebraException Raised if transform cannot be computed because
1178      *                          of numerical instabilities.
1179      */
1180     @Override
1181     public void transform(final Line2D inputLine, final Line2D outputLine) throws AlgebraException {
1182         // line' * point = 0 --> line' * t^-1 * t * point
1183         // (line' * t^-1)*(t*point) = (t^-1'*line)'*(t*point)
1184         // where:
1185         // - transformedLine = t^-1'*line
1186         // - transformedPoint = t*point
1187 
1188         inputLine.normalize();
1189         normalize();
1190 
1191         final var invT = inverseAndReturnNew().asMatrix();
1192         final var l = Matrix.newFromArray(inputLine.asArray());
1193 
1194         invT.transpose();
1195         invT.multiply(l);
1196 
1197         outputLine.setParameters(invT.toArray());
1198     }
1199 
1200     /**
1201      * Inverses this transformation.
1202      *
1203      * @throws AlgebraException if inverse transform cannot be computed because
1204      *                          of numerical instabilities.
1205      */
1206     public void inverse() throws AlgebraException {
1207         inverse(this);
1208     }
1209 
1210     /**
1211      * Computes the inverse of this transformation and returns the result as a
1212      * new transformation instance.
1213      *
1214      * @return inverse transformation.
1215      * @throws AlgebraException if inverse transform cannot be computed because
1216      *                          of numerical instabilities.
1217      */
1218     public Transformation2D inverseAndReturnNew() throws AlgebraException {
1219         final var result = new ProjectiveTransformation2D();
1220         inverse(result);
1221         return result;
1222     }
1223 
1224     /**
1225      * Computes the inverse of this transformation and stores the result in
1226      * provided instance.
1227      *
1228      * @param result Instance where inverse transformation will be stored.
1229      * @throws AlgebraException if inverse transform cannot be computed because
1230      *                          of numerical instabilities.
1231      */
1232     protected void inverse(final ProjectiveTransformation2D result) throws AlgebraException {
1233         result.t = Utils.inverse(t);
1234         result.normalized = false;
1235     }
1236 
1237     /**
1238      * Combines this transformation with provided transformation.
1239      * The combination is equivalent to multiplying the matrix of this
1240      * transformation with the matrix of provided transformation.
1241      *
1242      * @param transformation transformation to be combined with.
1243      */
1244     public void combine(final ProjectiveTransformation2D transformation) {
1245         combine(transformation, this);
1246     }
1247 
1248     /**
1249      * Combines this transformation with provided transformation and returns
1250      * the result as a new transformation instance.
1251      * The combination is equivalent to multiplying the matrix of this
1252      * transformation with the matrix of provided transformation.
1253      *
1254      * @param transformation transformation to be combined with.
1255      * @return a new transformation resulting of the combination with this
1256      * transformation and provided transformation.
1257      */
1258     public ProjectiveTransformation2D combineAndReturnNew(final ProjectiveTransformation2D transformation) {
1259         final var result = new ProjectiveTransformation2D();
1260         combine(transformation, result);
1261         return result;
1262     }
1263 
1264     /**
1265      * Combines this transformation with provided input transformation and
1266      * stores the result into provided output transformation.
1267      * The combination is equivalent to multiplying the matrix of this
1268      * transformation with the matrix of provided input transformation.
1269      *
1270      * @param inputTransformation  transformation to be combined with.
1271      * @param outputTransformation transformation where result will be stored.
1272      */
1273     private void combine(final ProjectiveTransformation2D inputTransformation,
1274                          final ProjectiveTransformation2D outputTransformation) {
1275         // combination in matrix representation is: T1 * T2
1276         normalize();
1277         inputTransformation.normalize();
1278 
1279         try {
1280             outputTransformation.t = this.t.multiplyAndReturnNew(inputTransformation.t);
1281             outputTransformation.normalized = false;
1282 
1283         } catch (final WrongSizeException ignore) {
1284             // never happens
1285         }
1286     }
1287 
1288     /**
1289      * Estimates this transformation internal matrix by providing 4
1290      * corresponding original and transformed points.
1291      *
1292      * @param inputPoint1  1st input point.
1293      * @param inputPoint2  2nd input point.
1294      * @param inputPoint3  3rd input point.
1295      * @param inputPoint4  4th input point.
1296      * @param outputPoint1 1st transformed point corresponding to 1st input
1297      *                     point.
1298      * @param outputPoint2 2nd transformed point corresponding to 2nd input
1299      *                     point.
1300      * @param outputPoint3 3rd transformed point corresponding to 3rd input
1301      *                     point.
1302      * @param outputPoint4 4th transformed point corresponding to 4th input
1303      *                     point.
1304      * @throws CoincidentPointsException raised if transformation cannot be
1305      *                                   estimated for some reason (point configuration degeneracy, duplicate
1306      *                                   points or numerical instabilities).
1307      */
1308     public final void setTransformationFromPoints(
1309             final Point2D inputPoint1, final Point2D inputPoint2, final Point2D inputPoint3, final Point2D inputPoint4,
1310             final Point2D outputPoint1, final Point2D outputPoint2, final Point2D outputPoint3,
1311             final Point2D outputPoint4) throws CoincidentPointsException {
1312 
1313         // normalize points to increase accuracy
1314         inputPoint1.normalize();
1315         inputPoint2.normalize();
1316         inputPoint3.normalize();
1317         inputPoint4.normalize();
1318 
1319         outputPoint1.normalize();
1320         outputPoint2.normalize();
1321         outputPoint3.normalize();
1322         outputPoint4.normalize();
1323 
1324         // matrix of homogeneous linear system of equations.
1325         // There are 9 unknowns and 8 equations (2 for each pair of corresponding
1326         // points)
1327         Matrix m = null;
1328         try {
1329             // build matrix initialized to zero
1330             m = new Matrix(8, 9);
1331 
1332             // 1st pair of points
1333             var iX = inputPoint1.getHomX();
1334             var iY = inputPoint1.getHomY();
1335             var iW = inputPoint1.getHomW();
1336 
1337             var oX = outputPoint1.getHomX();
1338             var oY = outputPoint1.getHomY();
1339             var oW = outputPoint1.getHomW();
1340 
1341             var oWiX = oW * iX;
1342             var oWiY = oW * iY;
1343             var oWiW = oW * iW;
1344 
1345             var oXiX = oX * iX;
1346             var oXiY = oX * iY;
1347             var oXiW = oX * iW;
1348 
1349             var oYiX = oY * iX;
1350             var oYiY = oY * iY;
1351             var oYiW = oY * iW;
1352 
1353             var tmp = oWiX * oWiX + oWiY * oWiY + oWiW * oWiW;
1354             var norm = Math.sqrt(tmp + oXiX * oXiX + oXiY * oXiY + oXiW * oXiW);
1355 
1356             m.setElementAt(0, 0, oWiX / norm);
1357             m.setElementAt(0, 1, oWiY / norm);
1358             m.setElementAt(0, 2, oWiW / norm);
1359             m.setElementAt(0, 6, -oXiX / norm);
1360             m.setElementAt(0, 7, -oXiY / norm);
1361             m.setElementAt(0, 8, -oXiW / norm);
1362 
1363             norm = Math.sqrt(tmp + oYiX * oYiX + oYiY * oYiY + oYiW * oYiW);
1364 
1365             m.setElementAt(1, 3, oWiX / norm);
1366             m.setElementAt(1, 4, oWiY / norm);
1367             m.setElementAt(1, 5, oWiW / norm);
1368             m.setElementAt(1, 6, -oYiX / norm);
1369             m.setElementAt(1, 7, -oYiY / norm);
1370             m.setElementAt(1, 8, -oYiW / norm);
1371 
1372             // 2nd pair of points
1373             iX = inputPoint2.getHomX();
1374             iY = inputPoint2.getHomY();
1375             iW = inputPoint2.getHomW();
1376 
1377             oX = outputPoint2.getHomX();
1378             oY = outputPoint2.getHomY();
1379             oW = outputPoint2.getHomW();
1380 
1381             oWiX = oW * iX;
1382             oWiY = oW * iY;
1383             oWiW = oW * iW;
1384 
1385             oXiX = oX * iX;
1386             oXiY = oX * iY;
1387             oXiW = oX * iW;
1388 
1389             oYiX = oY * iX;
1390             oYiY = oY * iY;
1391             oYiW = oY * iW;
1392 
1393             tmp = oWiX * oWiX + oWiY * oWiY + oWiW * oWiW;
1394             norm = Math.sqrt(tmp + oXiX * oXiX + oXiY * oXiY + oXiW * oXiW);
1395 
1396             m.setElementAt(2, 0, oWiX / norm);
1397             m.setElementAt(2, 1, oWiY / norm);
1398             m.setElementAt(2, 2, oWiW / norm);
1399             m.setElementAt(2, 6, -oXiX / norm);
1400             m.setElementAt(2, 7, -oXiY / norm);
1401             m.setElementAt(2, 8, -oXiW / norm);
1402 
1403             norm = Math.sqrt(tmp + oYiX * oYiX + oYiY * oYiY + oYiW * oYiW);
1404 
1405             m.setElementAt(3, 3, oWiX / norm);
1406             m.setElementAt(3, 4, oWiY / norm);
1407             m.setElementAt(3, 5, oWiW / norm);
1408             m.setElementAt(3, 6, -oYiX / norm);
1409             m.setElementAt(3, 7, -oYiY / norm);
1410             m.setElementAt(3, 8, -oYiW / norm);
1411 
1412             // 3rd pair of points
1413             iX = inputPoint3.getHomX();
1414             iY = inputPoint3.getHomY();
1415             iW = inputPoint3.getHomW();
1416 
1417             oX = outputPoint3.getHomX();
1418             oY = outputPoint3.getHomY();
1419             oW = outputPoint3.getHomW();
1420 
1421             oWiX = oW * iX;
1422             oWiY = oW * iY;
1423             oWiW = oW * iW;
1424 
1425             oXiX = oX * iX;
1426             oXiY = oX * iY;
1427             oXiW = oX * iW;
1428 
1429             oYiX = oY * iX;
1430             oYiY = oY * iY;
1431             oYiW = oY * iW;
1432 
1433             tmp = oWiX * oWiX + oWiY * oWiY + oWiW * oWiW;
1434             norm = Math.sqrt(tmp + oXiX * oXiX + oXiY * oXiY + oXiW * oXiW);
1435 
1436             m.setElementAt(4, 0, oWiX / norm);
1437             m.setElementAt(4, 1, oWiY / norm);
1438             m.setElementAt(4, 2, oWiW / norm);
1439             m.setElementAt(4, 6, -oXiX / norm);
1440             m.setElementAt(4, 7, -oXiY / norm);
1441             m.setElementAt(4, 8, -oXiW / norm);
1442 
1443             norm = Math.sqrt(tmp + oYiX * oYiX + oYiY * oYiY + oYiW * oYiW);
1444 
1445             m.setElementAt(5, 3, oWiX / norm);
1446             m.setElementAt(5, 4, oWiY / norm);
1447             m.setElementAt(5, 5, oWiW / norm);
1448             m.setElementAt(5, 6, -oYiX / norm);
1449             m.setElementAt(5, 7, -oYiY / norm);
1450             m.setElementAt(5, 8, -oYiW / norm);
1451 
1452             // 4th pair of points
1453             iX = inputPoint4.getHomX();
1454             iY = inputPoint4.getHomY();
1455             iW = inputPoint4.getHomW();
1456 
1457             oX = outputPoint4.getHomX();
1458             oY = outputPoint4.getHomY();
1459             oW = outputPoint4.getHomW();
1460 
1461             oWiX = oW * iX;
1462             oWiY = oW * iY;
1463             oWiW = oW * iW;
1464 
1465             oXiX = oX * iX;
1466             oXiY = oX * iY;
1467             oXiW = oX * iW;
1468 
1469             oYiX = oY * iX;
1470             oYiY = oY * iY;
1471             oYiW = oY * iW;
1472 
1473             tmp = oWiX * oWiX + oWiY * oWiY + oWiW * oWiW;
1474             norm = Math.sqrt(tmp + oXiX * oXiX + oXiY * oXiY + oXiW * oXiW);
1475 
1476             m.setElementAt(6, 0, oWiX / norm);
1477             m.setElementAt(6, 1, oWiY / norm);
1478             m.setElementAt(6, 2, oWiW / norm);
1479             m.setElementAt(6, 6, -oXiX / norm);
1480             m.setElementAt(6, 7, -oXiY / norm);
1481             m.setElementAt(6, 8, -oXiW / norm);
1482 
1483             norm = Math.sqrt(tmp + oYiX * oYiX + oYiY * oYiY + oYiW * oYiW);
1484 
1485             m.setElementAt(7, 3, oWiX / norm);
1486             m.setElementAt(7, 4, oWiY / norm);
1487             m.setElementAt(7, 5, oWiW / norm);
1488             m.setElementAt(7, 6, -oYiX / norm);
1489             m.setElementAt(7, 7, -oYiY / norm);
1490             m.setElementAt(7, 8, -oYiW / norm);
1491         } catch (final WrongSizeException ignore) {
1492             // never happens
1493         }
1494 
1495         // use SVD to decompose matrix m
1496         Matrix v;
1497         try {
1498             final var decomposer = new SingularValueDecomposer(m);
1499             decomposer.decompose();
1500 
1501             // ensure that matrix m has enough rank and there is a unique
1502             // solution (up to scale)
1503             if (decomposer.getRank() < 8) {
1504                 throw new CoincidentPointsException();
1505             }
1506             // V is 9x9
1507             v = decomposer.getV();
1508 
1509             // last column of V will contain parameters of transformation
1510             t.setSubmatrix(0, 0, HOM_COORDS - 1, HOM_COORDS - 1,
1511                     v.getSubmatrix(0, 8, 8, 8).toArray(), false);
1512             normalized = true; //because columns of V are normalized after SVD
1513 
1514         } catch (final AlgebraException e) {
1515             throw new CoincidentPointsException(e);
1516         }
1517     }
1518 
1519     /**
1520      * Estimates this transformation internal matrix by providing 4
1521      * corresponding original and transformed lines.
1522      *
1523      * @param inputLine1  1st input line.
1524      * @param inputLine2  2nd input line.
1525      * @param inputLine3  3rd input line.
1526      * @param inputLine4  4th input line.
1527      * @param outputLine1 1st transformed line corresponding to 1st input
1528      *                    line.
1529      * @param outputLine2 2nd transformed line corresponding to 2nd input
1530      *                    line.
1531      * @param outputLine3 3rd transformed line corresponding to 3rd input
1532      *                    line.
1533      * @param outputLine4 4th transformed line corresponding to 4th input
1534      *                    line.
1535      * @throws CoincidentLinesException raised if transformation cannot be
1536      *                                  estimated for some reason (line configuration degeneracy, duplicate
1537      *                                  line or numerical instabilities).
1538      */
1539     public final void setTransformationFromLines(
1540             final Line2D inputLine1, final Line2D inputLine2, final Line2D inputLine3, final Line2D inputLine4,
1541             final Line2D outputLine1, final Line2D outputLine2, final Line2D outputLine3, final Line2D outputLine4)
1542             throws CoincidentLinesException {
1543 
1544         // normalize lines to increase accuracy
1545         inputLine1.normalize();
1546         inputLine2.normalize();
1547         inputLine3.normalize();
1548         inputLine4.normalize();
1549 
1550         outputLine1.normalize();
1551         outputLine2.normalize();
1552         outputLine3.normalize();
1553         outputLine4.normalize();
1554 
1555         // matrix of homogeneous linear system of equations.
1556         // There are 9 unknowns and 8 equations (2 for each pair of corresponding
1557         // points)
1558         Matrix m = null;
1559         try {
1560             // build matrix initialized to zero
1561             m = new Matrix(8, 9);
1562 
1563             // 1st pair of lines
1564             var iA = inputLine1.getA();
1565             var iB = inputLine1.getB();
1566             var iC = inputLine1.getC();
1567 
1568             var oA = outputLine1.getA();
1569             var oB = outputLine1.getB();
1570             var oC = outputLine1.getC();
1571 
1572             var oCiA = oC * iA;
1573             var oCiB = oC * iB;
1574             var oCiC = oC * iC;
1575 
1576             var oAiA = oA * iA;
1577             var oAiB = oA * iB;
1578             var oAiC = oA * iC;
1579 
1580             var oBiA = oB * iA;
1581             var oBiB = oB * iB;
1582             var oBiC = oB * iC;
1583 
1584             var tmp = oCiA * oCiA + oCiB * oCiB + oCiC * oCiC;
1585             var norm = Math.sqrt(tmp + oAiA * oAiA + oAiB * oAiB + oAiC * oAiC);
1586 
1587             m.setElementAt(0, 0, oCiA / norm);
1588             m.setElementAt(0, 1, oCiB / norm);
1589             m.setElementAt(0, 2, oCiC / norm);
1590             m.setElementAt(0, 6, -oAiA / norm);
1591             m.setElementAt(0, 7, -oAiB / norm);
1592             m.setElementAt(0, 8, -oAiC / norm);
1593 
1594             norm = Math.sqrt(tmp + oBiA * oBiA + oBiB * oBiB + oBiC * oBiC);
1595 
1596             m.setElementAt(1, 3, oCiA / norm);
1597             m.setElementAt(1, 4, oCiB / norm);
1598             m.setElementAt(1, 5, oCiC / norm);
1599             m.setElementAt(1, 6, -oBiA / norm);
1600             m.setElementAt(1, 7, -oBiB / norm);
1601             m.setElementAt(1, 8, -oBiC / norm);
1602 
1603             // 2nd pair of lines
1604             iA = inputLine2.getA();
1605             iB = inputLine2.getB();
1606             iC = inputLine2.getC();
1607 
1608             oA = outputLine2.getA();
1609             oB = outputLine2.getB();
1610             oC = outputLine2.getC();
1611 
1612             oCiA = oC * iA;
1613             oCiB = oC * iB;
1614             oCiC = oC * iC;
1615 
1616             oAiA = oA * iA;
1617             oAiB = oA * iB;
1618             oAiC = oA * iC;
1619 
1620             oBiA = oB * iA;
1621             oBiB = oB * iB;
1622             oBiC = oB * iC;
1623 
1624             tmp = oCiA * oCiA + oCiB * oCiB + oCiC * oCiC;
1625             norm = Math.sqrt(tmp + oAiA * oAiA + oAiB * oAiB + oAiC * oAiC);
1626 
1627             m.setElementAt(2, 0, oCiA / norm);
1628             m.setElementAt(2, 1, oCiB / norm);
1629             m.setElementAt(2, 2, oCiC / norm);
1630             m.setElementAt(2, 6, -oAiA / norm);
1631             m.setElementAt(2, 7, -oAiB / norm);
1632             m.setElementAt(2, 8, -oAiC / norm);
1633 
1634             norm = Math.sqrt(tmp + oBiA * oBiA + oBiB * oBiB + oBiC * oBiC);
1635 
1636             m.setElementAt(3, 3, oCiA / norm);
1637             m.setElementAt(3, 4, oCiB / norm);
1638             m.setElementAt(3, 5, oCiC / norm);
1639             m.setElementAt(3, 6, -oBiA / norm);
1640             m.setElementAt(3, 7, -oBiB / norm);
1641             m.setElementAt(3, 8, -oBiC / norm);
1642 
1643             // 3rd pair of points
1644             iA = inputLine3.getA();
1645             iB = inputLine3.getB();
1646             iC = inputLine3.getC();
1647 
1648             oA = outputLine3.getA();
1649             oB = outputLine3.getB();
1650             oC = outputLine3.getC();
1651 
1652             oCiA = oC * iA;
1653             oCiB = oC * iB;
1654             oCiC = oC * iC;
1655 
1656             oAiA = oA * iA;
1657             oAiB = oA * iB;
1658             oAiC = oA * iC;
1659 
1660             oBiA = oB * iA;
1661             oBiB = oB * iB;
1662             oBiC = oB * iC;
1663 
1664             tmp = oCiA * oCiA + oCiB * oCiB + oCiC * oCiC;
1665             norm = Math.sqrt(tmp + oAiA * oAiA + oAiB * oAiB + oAiC * oAiC);
1666 
1667             m.setElementAt(4, 0, oCiA / norm);
1668             m.setElementAt(4, 1, oCiB / norm);
1669             m.setElementAt(4, 2, oCiC / norm);
1670             m.setElementAt(4, 6, -oAiA / norm);
1671             m.setElementAt(4, 7, -oAiB / norm);
1672             m.setElementAt(4, 8, -oAiC / norm);
1673 
1674             norm = Math.sqrt(tmp + oBiA * oBiA + oBiB * oBiB + oBiC * oBiC);
1675 
1676             m.setElementAt(5, 3, oCiA / norm);
1677             m.setElementAt(5, 4, oCiB / norm);
1678             m.setElementAt(5, 5, oCiC / norm);
1679             m.setElementAt(5, 6, -oBiA / norm);
1680             m.setElementAt(5, 7, -oBiB / norm);
1681             m.setElementAt(5, 8, -oBiC / norm);
1682 
1683             // 4th pair of points
1684             iA = inputLine4.getA();
1685             iB = inputLine4.getB();
1686             iC = inputLine4.getC();
1687 
1688             oA = outputLine4.getA();
1689             oB = outputLine4.getB();
1690             oC = outputLine4.getC();
1691 
1692             oCiA = oC * iA;
1693             oCiB = oC * iB;
1694             oCiC = oC * iC;
1695 
1696             oAiA = oA * iA;
1697             oAiB = oA * iB;
1698             oAiC = oA * iC;
1699 
1700             oBiA = oB * iA;
1701             oBiB = oB * iB;
1702             oBiC = oB * iC;
1703 
1704             tmp = oCiA * oCiA + oCiB * oCiB + oCiC * oCiC;
1705             norm = Math.sqrt(tmp + oAiA * oAiA + oAiB * oAiB + oAiC * oAiC);
1706 
1707             m.setElementAt(6, 0, oCiA / norm);
1708             m.setElementAt(6, 1, oCiB / norm);
1709             m.setElementAt(6, 2, oCiC / norm);
1710             m.setElementAt(6, 6, -oAiA / norm);
1711             m.setElementAt(6, 7, -oAiB / norm);
1712             m.setElementAt(6, 8, -oAiC / norm);
1713 
1714             norm = Math.sqrt(tmp + oBiA * oBiA + oBiB * oBiB + oBiC * oBiC);
1715 
1716             m.setElementAt(7, 3, oCiA / norm);
1717             m.setElementAt(7, 4, oCiB / norm);
1718             m.setElementAt(7, 5, oCiC / norm);
1719             m.setElementAt(7, 6, -oBiA / norm);
1720             m.setElementAt(7, 7, -oBiB / norm);
1721             m.setElementAt(7, 8, -oBiC / norm);
1722         } catch (final WrongSizeException ignore) {
1723             // never happens
1724         }
1725 
1726         // use SVD to decompose matrix m
1727         Matrix v;
1728         try {
1729             final var decomposer = new SingularValueDecomposer(m);
1730             decomposer.decompose();
1731 
1732             // ensure that matrix m has enough rank and there is a unique
1733             // solution (up to scale)
1734             if (decomposer.getRank() < 8) {
1735                 throw new CoincidentLinesException();
1736             }
1737             // V is 9x9
1738             v = decomposer.getV();
1739 
1740             // last column of V will contain parameters of transformation
1741             final var transInvT = new Matrix(HOM_COORDS, HOM_COORDS);
1742             transInvT.setSubmatrix(0, 0, HOM_COORDS - 1,
1743                     HOM_COORDS - 1, v.getSubmatrix(0, 8, 8,
1744                             8).toArray(), false);
1745             // this is now invT
1746             transInvT.transpose();
1747             t = Utils.inverse(transInvT);
1748             // invT is normalized, but not t
1749             normalized = false;
1750 
1751         } catch (final AlgebraException e) {
1752             throw new CoincidentLinesException(e);
1753         }
1754     }
1755 }