1 /*
2 * Copyright (C) 2016 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.ar.calibration.estimators;
17
18 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.CholeskyDecomposer;
20 import com.irurueta.algebra.Complex;
21 import com.irurueta.algebra.SingularValueDecomposer;
22 import com.irurueta.ar.calibration.DualImageOfAbsoluteConic;
23 import com.irurueta.ar.epipolar.FundamentalMatrix;
24 import com.irurueta.geometry.PinholeCameraIntrinsicParameters;
25 import com.irurueta.geometry.estimators.LockedException;
26 import com.irurueta.geometry.estimators.NotReadyException;
27 import com.irurueta.numerical.NumericalException;
28 import com.irurueta.numerical.polynomials.Polynomial;
29
30 /**
31 * Estimates the DIAC (Dual Image of Absolute Conic) by solving Kruppa's
32 * equations and assuming known principal point and zero skewness.
33 * This estimator allows enforcing a known aspect ratio as well.
34 * The DIAC can be used to obtain the intrinsic parameters of a pair of
35 * cameras related by a fundamental matrix.
36 * Hence, this class can be used for auto-calibration purposes.
37 * Notice that the {@link DualAbsoluteQuadricEstimator} is a more robust method of
38 * auto-calibration.
39 * This class is based on:
40 * S.D. Hippisley-Cox & J.Porrill. Auto-calibration - Kruppa's equations and the
41 * intrinsic parameters of a camera. AI Vision Research Unit. University of
42 * Sheffield.
43 */
44 @SuppressWarnings("DuplicatedCode")
45 public class KruppaDualImageOfAbsoluteConicEstimator {
46
47 /**
48 * Degree of polynomial to solve Kruppa's equation when aspect ratio is
49 * known.
50 */
51 private static final int POLY_DEGREE_UNKNOWN_ASPECT_RATIO = 4;
52
53 /**
54 * Default value for horizontal principal point coordinate.
55 */
56 public static final double DEFAULT_PRINCIPAL_POINT_X = 0.0;
57
58 /**
59 * Default value for vertical principal point coordinate.
60 */
61 public static final double DEFAULT_PRINCIPAL_POINT_Y = 0.0;
62
63 /**
64 * Constant defining whether aspect ratio of focal distance (i.e. vertical
65 * focal distance divided by horizontal focal distance) is known or not.
66 * Notice that focal distance aspect ratio is not related to image size
67 * aspect ratio. Typically, LCD sensor cells are square and hence aspect
68 * ratio of focal distances is known and equal to 1.
69 */
70 public static final boolean DEFAULT_FOCAL_DISTANCE_ASPECT_RATIO_KNOWN = true;
71
72 /**
73 * Constant defining default aspect ratio of focal distances. This constant
74 * takes into account that typically LCD sensor cells are square and hence
75 * aspect ratio of focal distances is known and equal to 1.
76 */
77 public static final double DEFAULT_FOCAL_DISTANCE_ASPECT_RATIO = 1.0;
78
79 /**
80 * Minimum absolute value allowed for aspect ratio of focal distances.
81 */
82 public static final double MIN_ABS_FOCAL_DISTANCE_ASPECT_RATIO = 1e-6;
83
84 /**
85 * Known horizontal principal point coordinate.
86 */
87 private double principalPointX;
88
89 /**
90 * Known vertical principal point coordinate.
91 */
92 private double principalPointY;
93
94 /**
95 * Indicates whether aspect ratio of focal distances (i.e. vertical focal
96 * distance divided by horizontal focal distance) is known or not.
97 * Notice that focal distance aspect ratio is not related to image aspect
98 * ratio. Typically, LCD sensor cells are square and hence aspect ratio of
99 * focal distances is known and equal to 1.
100 */
101 private boolean focalDistanceAspectRatioKnown;
102
103 /**
104 * Contains aspect ratio of focal distances (i.e. vertical focal distance
105 * divided by horizontal focal distance).
106 * By default, this is 1.0, since it is taken into account that typically
107 * LCD sensor cells are square and hence aspect ratio focal distance is
108 * known and equal to 1.
109 * Notice that focal distance aspect ratio is not related to image size
110 * aspect ratio.
111 */
112 private double focalDistanceAspectRatio;
113
114 /**
115 * True when estimator is estimating the DIAC.
116 */
117 private boolean locked;
118
119 /**
120 * Listener to be notified of events such as when estimation starts, ends or
121 * estimation progress changes.
122 */
123 private KruppaDualImageOfAbsoluteConicEstimatorListener listener;
124
125 /**
126 * Fundamental matrix to estimate DIAC from.
127 */
128 private FundamentalMatrix fundamentalMatrix;
129
130 /**
131 * Constructor.
132 */
133 public KruppaDualImageOfAbsoluteConicEstimator() {
134 principalPointX = DEFAULT_PRINCIPAL_POINT_X;
135 principalPointY = DEFAULT_PRINCIPAL_POINT_Y;
136 focalDistanceAspectRatioKnown = DEFAULT_FOCAL_DISTANCE_ASPECT_RATIO_KNOWN;
137 focalDistanceAspectRatio = DEFAULT_FOCAL_DISTANCE_ASPECT_RATIO;
138
139 locked = false;
140 listener = null;
141 fundamentalMatrix = null;
142 }
143
144 /**
145 * Constructor.
146 *
147 * @param listener listener to be notified of events such as when estimation
148 * starts, ends or estimation progress changes.
149 */
150 public KruppaDualImageOfAbsoluteConicEstimator(final KruppaDualImageOfAbsoluteConicEstimatorListener listener) {
151 this();
152 this.listener = listener;
153 }
154
155 /**
156 * Constructor.
157 *
158 * @param fundamentalMatrix fundamental matrix to estimate DIAC from.
159 */
160 public KruppaDualImageOfAbsoluteConicEstimator(final FundamentalMatrix fundamentalMatrix) {
161 this();
162 this.fundamentalMatrix = fundamentalMatrix;
163 }
164
165 /**
166 * Constructor.
167 *
168 * @param fundamentalMatrix fundamental matrix to estimate DIAC from.
169 * @param listener listener to be notified of events such as when estimation
170 * starts, ends or estimation progress changes.
171 */
172 public KruppaDualImageOfAbsoluteConicEstimator(
173 final FundamentalMatrix fundamentalMatrix, final KruppaDualImageOfAbsoluteConicEstimatorListener listener) {
174 this(fundamentalMatrix);
175 this.listener = listener;
176 }
177
178 /**
179 * Gets known horizontal principal point coordinate.
180 *
181 * @return known horizontal principal point coordinate.
182 */
183 public double getPrincipalPointX() {
184 return principalPointX;
185 }
186
187 /**
188 * Sets known horizontal principal point coordinate.
189 *
190 * @param principalPointX known horizontal principal point coordinate to be
191 * set.
192 * @throws LockedException if estimator is locked.
193 */
194 public void setPrincipalPointX(final double principalPointX) throws LockedException {
195 if (isLocked()) {
196 throw new LockedException();
197 }
198 this.principalPointX = principalPointX;
199 }
200
201 /**
202 * Gets known vertical principal point coordinate.
203 *
204 * @return known vertical principal point coordinate.
205 */
206 public double getPrincipalPointY() {
207 return principalPointY;
208 }
209
210 /**
211 * Sets known vertical principal point coordinate.
212 *
213 * @param principalPointY known vertical principal point coordinate to be
214 * set.
215 * @throws LockedException if estimator is locked.
216 */
217 public void setPrincipalPointY(final double principalPointY) throws LockedException {
218 if (isLocked()) {
219 throw new LockedException();
220 }
221 this.principalPointY = principalPointY;
222 }
223
224 /**
225 * Returns boolean indicating whether aspect ratio of focal distances (i.e.
226 * vertical focal distance divided by horizontal focal distance) is known or
227 * not.
228 * Notice that focal distance aspect ratio is not related to image size
229 * aspect ratio. Typically, LCD sensor cells are square and hence aspect
230 * ratio of focal distances is known and equal to 1.
231 * This value is only taken into account if skewness is assumed to be zero,
232 * otherwise it is ignored.
233 *
234 * @return true if focal distance aspect ratio is known, false otherwise.
235 */
236 public boolean isFocalDistanceAspectRatioKnown() {
237 return focalDistanceAspectRatioKnown;
238 }
239
240 /**
241 * Sets value indicating whether aspect ratio of focal distances (i.e.
242 * vertical focal distance divided by horizontal focal distance) is known or
243 * not.
244 * Notice that focal distance aspect ratio is not related to image size
245 * aspect ratio. Typically, LCD sensor cells are square and hence aspect
246 * ratio of focal distances is known and equal to 1.
247 * This value is only taken into account if skewness is assumed to be
248 * otherwise it is ignored.
249 *
250 * @param focalDistanceAspectRatioKnown true if focal distance aspect ratio
251 * is known, false otherwise.
252 * @throws LockedException if estimator is locked.
253 */
254 public void setFocalDistanceAspectRatioKnown(final boolean focalDistanceAspectRatioKnown) throws LockedException {
255 if (isLocked()) {
256 throw new LockedException();
257 }
258
259 this.focalDistanceAspectRatioKnown = focalDistanceAspectRatioKnown;
260 }
261
262 /**
263 * Returns aspect ratio of focal distances (i.e. vertical focal distance
264 * divided by horizontal focal distance).
265 * By default, this is 1.0, since it is taken into account that typically
266 * LCD sensor cells are square and hence aspect ratio focal distances is
267 * known and equal to 1.
268 * Notice that focal distance aspect ratio is not related to image size
269 * aspect ratio
270 * Notice that a negative aspect ratio indicates that vertical axis is
271 * reversed. This can be useful in some situations where image vertical
272 * coordinates are reversed respect to the physical world (i.e. in computer
273 * graphics typically image vertical coordinates go downwards, while in
274 * physical world they go upwards).
275 *
276 * @return aspect ratio of focal distances.
277 */
278 public double getFocalDistanceAspectRatio() {
279 return focalDistanceAspectRatio;
280 }
281
282 /**
283 * Sets aspect ratio of focal distances (i.e. vertical focal distance
284 * divided by horizontal focal distance).
285 * This value is only taken into account if aspect ratio is marked as known,
286 * otherwise it is ignored.
287 * By default, this is 1.0, since it is taken into account that typically
288 * LCD sensor cells are square and hence aspect ratio focal distances is
289 * known and equal to 1.
290 * Notice that focal distance aspect ratio is not related to image size
291 * aspect ratio.
292 * Notice that a negative aspect ratio indicates that vertical axis is
293 * reversed. This can be useful in some situations where image vertical
294 * coordinates are reversed respect to the physical world (i.e. in computer
295 * graphics typically image vertical coordinates go downwards, while in
296 * physical world they go upwards).
297 *
298 * @param focalDistanceAspectRatio aspect ratio of focal distances to be
299 * set.
300 * @throws LockedException if estimator is locked.
301 * @throws IllegalArgumentException if focal distance aspect ratio is too
302 * close to zero, as it might produce numerical instabilities.
303 */
304 public void setFocalDistanceAspectRatio(final double focalDistanceAspectRatio) throws LockedException {
305 if (isLocked()) {
306 throw new LockedException();
307 }
308 if (Math.abs(focalDistanceAspectRatio) < MIN_ABS_FOCAL_DISTANCE_ASPECT_RATIO) {
309 throw new IllegalArgumentException();
310 }
311
312 this.focalDistanceAspectRatio = focalDistanceAspectRatio;
313 }
314
315 /**
316 * Indicates whether this instance is locked.
317 *
318 * @return true if this estimator is busy estimating a DIAC, false
319 * otherwise.
320 */
321 public boolean isLocked() {
322 return locked;
323 }
324
325 /**
326 * Returns listener to be notified of events such as when estimation starts,
327 * ends or estimation progress changes.
328 *
329 * @return listener to be notified of events.
330 */
331 public KruppaDualImageOfAbsoluteConicEstimatorListener getListener() {
332 return listener;
333 }
334
335 /**
336 * Sets listener to be notified of events such as when estimation starts,
337 * ends or estimation progress changes.
338 *
339 * @param listener listener to be notified of events.
340 * @throws LockedException if estimator is locked.
341 */
342 public void setListener(
343 final KruppaDualImageOfAbsoluteConicEstimatorListener listener) throws LockedException {
344 if (isLocked()) {
345 throw new LockedException();
346 }
347 this.listener = listener;
348 }
349
350 /**
351 * Gets fundamental matrix to estimate DIAC from.
352 *
353 * @return fundamental matrix to estimate DIAC from.
354 */
355 public FundamentalMatrix getFundamentalMatrix() {
356 return fundamentalMatrix;
357 }
358
359 /**
360 * Sets fundamental matrix to estimate DIAC from.
361 *
362 * @param fundamentalMatrix fundamental matrix to estimate DIAC from.
363 * @throws LockedException if estimator is locked.
364 */
365 public void setFundamentalMatrix(final FundamentalMatrix fundamentalMatrix) throws LockedException {
366 if (isLocked()) {
367 throw new LockedException();
368 }
369 this.fundamentalMatrix = fundamentalMatrix;
370 }
371
372 /**
373 * Returns value indicating whether required data has been provided so that
374 * DIAC estimation can start.
375 * If true, estimator is ready to compute the DIAC, otherwise more data
376 * needs to be provided.
377 *
378 * @return true if estimator is ready, false otherwise.
379 */
380 public boolean isReady() {
381 return fundamentalMatrix != null;
382 }
383
384 /**
385 * Estimates Dual Image of Absolute Conic (DIAC).
386 *
387 * @return estimated DIAC.
388 * @throws LockedException if estimator is locked.
389 * @throws NotReadyException if input has not yet been provided.
390 * @throws KruppaDualImageOfAbsoluteConicEstimatorException if an error
391 * occurs during estimation, usually because
392 * fundamental matrix corresponds to
393 * degenerate camera movements, or because of
394 * numerical instabilities.
395 */
396 public DualImageOfAbsoluteConic estimate() throws LockedException, NotReadyException,
397 KruppaDualImageOfAbsoluteConicEstimatorException {
398 final var result = new DualImageOfAbsoluteConic(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
399 estimate(result);
400 return result;
401 }
402
403 /**
404 * Estimates Dual Image of Absolute Conic (DIAC).
405 *
406 * @param result instance where estimated DIAC will be stored.
407 * @throws LockedException if estimator is locked.
408 * @throws NotReadyException if input has not yet been provided.
409 * @throws KruppaDualImageOfAbsoluteConicEstimatorException if an error
410 * occurs during estimation, usually because
411 * fundamental matrix corresponds to
412 * degenerate camera movements, or because of
413 * numerical instabilities.
414 */
415 public void estimate(final DualImageOfAbsoluteConic result) throws LockedException, NotReadyException,
416 KruppaDualImageOfAbsoluteConicEstimatorException {
417 if (isLocked()) {
418 throw new LockedException();
419 }
420 if (!isReady()) {
421 throw new NotReadyException();
422 }
423
424 try {
425 locked = true;
426
427 if (listener != null) {
428 listener.onEstimateStart(this);
429 }
430
431 if (focalDistanceAspectRatioKnown) {
432 estimateKnownAspectRatio(result);
433 } else {
434 estimateUnknownAspectRatio(result);
435 }
436
437 if (listener != null) {
438 listener.onEstimateEnd(this);
439 }
440 } finally {
441 locked = false;
442 }
443 }
444
445 /**
446 * Builds a DIAC from estimated focal length components and current
447 * principal point assuming zero skewness.
448 *
449 * @param horizontalFocalLength estimated horizontal focal length component.
450 * @param verticalFocalLength estimated vertical focal length component.
451 * @param result instance where estimated DIAC will be stored.
452 * @return true if estimated DIAC is valid, false otherwise.
453 */
454 private boolean buildDiac(final double horizontalFocalLength, final double verticalFocalLength,
455 final DualImageOfAbsoluteConic result) {
456
457 try {
458 final var intrinsic = new PinholeCameraIntrinsicParameters(horizontalFocalLength, verticalFocalLength,
459 principalPointX, principalPointY, 0.0);
460 result.setFromPinholeCameraIntrinsicParameters(intrinsic);
461
462 final var m = result.asMatrix();
463 final var decomposer = new CholeskyDecomposer(m);
464 decomposer.decompose();
465 return decomposer.isSPD();
466 } catch (final AlgebraException e) {
467 // there are numerical instabilities
468 return false;
469 }
470 }
471
472 /**
473 * Estimates the DIAC assuming unknown aspect ratio.
474 *
475 * @param result instance where estimated DIAC will be stored.
476 * @throws KruppaDualImageOfAbsoluteConicEstimatorException if an error
477 * occurs during estimation, usually because
478 * fundamental matrix corresponds to
479 * degenerate camera movements, or because of
480 * numerical instabilities.
481 */
482 private void estimateUnknownAspectRatio(final DualImageOfAbsoluteConic result)
483 throws KruppaDualImageOfAbsoluteConicEstimatorException {
484 try {
485 final var x0 = principalPointX;
486 final var y0 = principalPointY;
487
488 // SVD decompose fundamental matrix
489 fundamentalMatrix.normalize();
490 final var decomposer = new SingularValueDecomposer(fundamentalMatrix.getInternalMatrix());
491 decomposer.decompose();
492
493 final var sigmas = decomposer.getSingularValues();
494 final var u = decomposer.getU();
495 final var v = decomposer.getV();
496
497 final var sigma1 = sigmas[0];
498 final var sigma2 = sigmas[1];
499
500 // Column u1
501 final var u11 = u.getElementAt(0, 0);
502 final var u21 = u.getElementAt(1, 0);
503 final var u31 = u.getElementAt(2, 0);
504
505 // Column u2
506 final var u12 = u.getElementAt(0, 1);
507 final var u22 = u.getElementAt(1, 1);
508 final var u32 = u.getElementAt(2, 1);
509
510 // Column v1
511 final var v11 = v.getElementAt(0, 0);
512 final var v21 = v.getElementAt(1, 0);
513 final var v31 = v.getElementAt(2, 0);
514
515 // Column v2
516 final var v12 = v.getElementAt(0, 1);
517 final var v22 = v.getElementAt(1, 1);
518 final var v32 = v.getElementAt(2, 1);
519
520 // build Kruppa equations
521 final var polyA = u12 * u11;
522 final var polyB = u22 * u21;
523 final var polyC = Math.pow(x0, 2.0) * u12 * u11 + x0 * y0 * u22 * u11 + x0 * u32 * u11
524 + x0 * y0 * u12 * u21 + Math.pow(y0, 2.0) * u22 * u21 + y0 * u32 * u21
525 + x0 * u12 * u31 + y0 * u22 * u31 + u32 * u31;
526 final var polyD = Math.pow(sigma2, 2.0) * v12 * v12;
527 final var polyE = Math.pow(sigma2, 2.0) * v22 * v22;
528 final var polyF = Math.pow(sigma2 * x0, 2.0) * v12 * v12
529 + Math.pow(sigma2, 2.0) * x0 * y0 * v22 * v12
530 + Math.pow(sigma2, 2.0) * x0 * v32 * v12
531 + Math.pow(sigma2, 2.0) * x0 * y0 * v12 * v22
532 + Math.pow(sigma2 * y0, 2.0) * v22 * v22
533 + Math.pow(sigma2, 2.0) * y0 * v32 * v22
534 + Math.pow(sigma2, 2.0) * x0 * v12 * v32
535 + Math.pow(sigma2, 2.0) * y0 * v22 * v32
536 + Math.pow(sigma2, 2.0) * v32 * v32;
537 final var polyG = u11 * u11;
538 final var polyH = u21 * u21;
539 final var polyI = Math.pow(x0, 2.0) * u11 * u11 + x0 * y0 * u21 * u11 + x0 * u31 * u11
540 + x0 * y0 * u11 * u21 + Math.pow(y0, 2.0) * u21 * u21 + y0 * u31 * u21
541 + x0 * u11 * u31 + y0 * u21 * u31 + u31 * u31;
542 final var polyJ = sigma1 * sigma2 * v12 * v11;
543 final var polyK = sigma1 * sigma2 * v22 * v21;
544 final var polyL = sigma1 * sigma2 * Math.pow(x0, 2.0) * v12 * v11
545 + sigma1 * sigma2 * x0 * y0 * v22 * v11 + sigma1 * sigma2 * x0 * v32 * v11
546 + sigma1 * sigma2 * x0 * y0 * v12 * v21
547 + sigma1 * sigma2 * Math.pow(y0, 2.0) * v22 * v21
548 + sigma1 * sigma2 * y0 * v32 * v21 + sigma1 * sigma2 * x0 * v12 * v31
549 + sigma1 * sigma2 * y0 * v22 * v31 + sigma1 * sigma2 * v32 * v31;
550 final var polyM = Math.pow(sigma1, 2.0) * v11 * v11;
551 final var polyN = Math.pow(sigma1, 2.0) * v21 * v21;
552 final var polyO = Math.pow(sigma1 * x0, 2.0) * v11 * v11
553 + Math.pow(sigma1, 2.0) * x0 * y0 * v21 * v11
554 + Math.pow(sigma1, 2.0) * x0 * v31 * v11
555 + Math.pow(sigma1, 2.0) * x0 * y0 * v11 * v21
556 + Math.pow(sigma1 * y0, 2.0) * v21 * v21
557 + Math.pow(sigma1, 2.0) * y0 * v31 * v21
558 + Math.pow(sigma1, 2.0) * x0 * v11 * v31
559 + Math.pow(sigma1, 2.0) * y0 * v21 * v31
560 + Math.pow(sigma1, 2.0) * v31 * v31;
561 final var polyP = u12 * u12;
562 final var polyQ = u22 * u22;
563 final var polyR = Math.pow(x0, 2.0) * u12 * u12 + x0 * y0 * u22 * u12 + x0 * u32 * u12
564 + x0 * y0 * u12 * u22 + Math.pow(y0, 2.0) * u22 * u22 + y0 * u32 * u22
565 + x0 * u12 * u32 + y0 * u22 * u32 + u32 * u32;
566
567
568 final var tmp = (polyP * polyJ + polyA * polyM) / (polyG * polyM - polyP * polyD);
569 final var polyS = (tmp * (polyH * polyN - polyQ * polyE) - (polyQ * polyK + polyB * polyN));
570 final var polyT = (tmp * (polyG * polyN + polyH * polyM - polyP * polyE - polyQ * polyD)
571 - (polyP * polyK + polyQ * polyJ + polyA * polyN + polyB * polyM));
572 final var polyU = (tmp * (polyG * polyO + polyM * polyI - polyP * polyF - polyD * polyR)
573 - (polyP * polyL + polyJ * polyR + polyA * polyO + polyM * polyC));
574 final var polyV = (tmp * (polyH * polyO + polyN * polyI - polyQ * polyF - polyE * polyR)
575 - (polyQ * polyL + polyK * polyR + polyB * polyO + polyN * polyC));
576 final var polyW = (tmp * (polyO * polyI - polyF * polyR) - (polyL * polyR + polyO * polyC));
577
578 // assuming that x = ax^2, y = ay^2 which are the horizontal and
579 // vertical focal lengths, we obtain the following equations
580
581 // x = (-y^2*S - y*V - W) / (y*T + U)
582 // (-y^2*S - y*V - W)^2 *(-A*D - G*J) + y^2*(y*T + U)^2*(-B*E - H*K) +
583 // (-y^2*S - y*V - W)*(y*T + U)*y*(-A*E - B*D - G*K - H*J) +
584 // (-y^2*S - y*V - W)*(y*T + U)*(-A*F - D*C - G*L - J*I) +
585 // y*(y*T + U)^2*(-B*F - E*C - H*L - K*I) +
586 // (y*T + U)^2*(- F*C - L*I) = 0
587 // (-y^2*S - y*V - W)^2*(G*M - P*D) + y^2*(y*T + U)^2*(H*N - Q*E) +
588 // (-y^2*S - y*V - W)*(y*T + U)*y*(G*N + H*M - P*E - Q*D) +
589 // (-y^2*S - y*V - W)*(y*T + U)*(G*O + M*I - P*F - D*R) +
590 // y*(y*T + U)^2*(H*O + N*I - Q*F - E*R) + (y*T + U)^2*(O*I - F*R) = 0
591
592 // where we can solve y using any of the two latter equations, and
593 // then use obtained y to solve x
594 final var roots = unknownAspectRatioRoots(polyA, polyB, polyC, polyD, polyE, polyF, polyG, polyH, polyI,
595 polyJ, polyK, polyL, polyM, polyN, polyO, polyP, polyQ, polyR, polyS, polyT, polyU, polyV, polyW);
596
597 // roots contain possible y values. We use only their real part
598 // and find x = (-y^2*S - y*V - W) / (y*T + U)
599
600 // pick the best x, y values that produce a positive definite DIAC
601 // matrix
602 final var diac = new DualImageOfAbsoluteConic(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
603 var valid = false;
604 if (roots != null) {
605 for (final var root : roots) {
606 final var y = root.getReal();
607 final var x = getXFromY(y, polyS, polyT, polyU, polyV, polyW);
608
609 // build DIAC matrix and check if it is positive definite
610 if (x >= 0.0 && y >= 0.0) {
611 final var horizontalFocalLength = Math.sqrt(x);
612 final var verticalFocalLength = Math.sqrt(y);
613 valid = buildDiac(horizontalFocalLength, verticalFocalLength, diac);
614 }
615
616 if (valid) {
617 break;
618 }
619 }
620 } else {
621 throw new KruppaDualImageOfAbsoluteConicEstimatorException();
622 }
623
624 if (valid) {
625 // copy to result
626 result.setParameters(diac.getA(), diac.getB(), diac.getC(), diac.getD(), diac.getE(), diac.getF());
627 } else {
628 // no valid DIAC could be found
629 throw new KruppaDualImageOfAbsoluteConicEstimatorException();
630 }
631
632 } catch (final KruppaDualImageOfAbsoluteConicEstimatorException ex) {
633 throw ex;
634 } catch (final Exception ex) {
635 throw new KruppaDualImageOfAbsoluteConicEstimatorException(ex);
636 }
637 }
638
639 /**
640 * Gets x value from current y value.
641 * X and y values are the squared values of estimated focal length
642 * components.
643 * This method is used internally when aspect ratio is not known.
644 *
645 * @param y y value to obtain x value from.
646 * @param s internal value from Kruppa's equations.
647 * @param t internal value from Kruppa's equations.
648 * @param u internal value from Kruppa's equations.
649 * @param v internal value from Kruppa's equations.
650 * @param w internal value from Kruppa's equations.
651 * @return x value.
652 */
653 private double getXFromY(final double y, final double s, final double t, final double u, final double v,
654 final double w) {
655 return (-Math.pow(y, 2.0) * s - y * v - w) / (y * t + u);
656 }
657
658 /**
659 * One of Kruppa's equations expressed as a polynomial of degree 4 to solve
660 * y value, which is the squared value of vertical focal length.
661 * This method is only used when aspect ratio is unknown.
662 *
663 * @param a internal value from Kruppa's equations.
664 * @param b internal value from Kruppa's equations.
665 * @param c internal value from Kruppa's equations.
666 * @param d internal value from Kruppa's equations.
667 * @param e internal value from Kruppa's equations.
668 * @param f internal value from Kruppa's equations.
669 * @param g internal value from Kruppa's equations.
670 * @param h internal value from Kruppa's equations.
671 * @param i internal value from Kruppa's equations.
672 * @param j internal value from Kruppa's equations.
673 * @param k internal value from Kruppa's equations.
674 * @param l internal value from Kruppa's equations.
675 * @param s internal value from Kruppa's equations.
676 * @param t internal value from Kruppa's equations.
677 * @param u internal value from Kruppa's equations.
678 * @param v internal value from Kruppa's equations.
679 * @param w internal value from Kruppa's equations.
680 * @return a polynomial.
681 */
682 private Polynomial buildPolynomial1(
683 final double a, final double b, final double c, final double d,
684 final double e, final double f, final double g, final double h,
685 final double i, final double j, final double k, final double l,
686 final double s, final double t, final double u, final double v,
687 final double w) {
688 // (-y^2*S - y*V - W)^2 *(-A*D - G*J) + y^2*(y*T + U)^2*(-B*E - H*K) +
689 // (-y^2*S - y*V - W)*(y*T + U)*y*(-A*E - B*D - G*K - H*J) +
690 // (-y^2*S - y*V - W)*(y*T + U)*(-A*F - D*C - G*L - J*I) +
691 // y*(y*T + U)^2*(-B*F - E*C - H*L - K*I) + (y*T + U)^2*(- F*C - L*I) = 0
692
693 final var result = new Polynomial(POLY_DEGREE_UNKNOWN_ASPECT_RATIO + 1);
694
695 // (-y^2*S - y*V - W)^2 *(-A*D - G*J)
696 final var tmp = new Polynomial(-w, -v, -s);
697 final var tmp2 = new Polynomial(-w, -v, -s);
698 tmp.multiply(tmp2);
699 tmp.multiplyByScalar(-a * d - g * j);
700 result.add(tmp);
701
702 // y^2*(y*T + U)^2*(-B*E - H*K)
703 tmp.setPolyParams(0.0, 0.0, 1.0);
704 tmp2.setPolyParams(u, t);
705 tmp.multiply(tmp2);
706 tmp.multiply(tmp2);
707 tmp.multiplyByScalar(-b * e - h * k);
708 result.add(tmp);
709
710 // (-y^2*S - y*V - W)*(y*T + U)*y*(-A*E - B*D - G*K - H*J)
711 tmp.setPolyParams(-w, -v, -s);
712 tmp2.setPolyParams(u, t);
713 tmp.multiply(tmp2);
714 tmp2.setPolyParams(0.0, 1);
715 tmp.multiply(tmp2);
716 tmp.multiplyByScalar(-a * e - b * d - g * k - h * j);
717 result.add(tmp);
718
719 // (-y^2*S - y*V - W)*(y*T + U)*(-A*F - D*C - G*L - J*I)
720 tmp.setPolyParams(-w, -v, -s);
721 tmp2.setPolyParams(u, t);
722 tmp.multiply(tmp2);
723 tmp.multiplyByScalar(-a * f - d * c - g * l - j * i);
724 result.add(tmp);
725
726 // y*(y*T + U)^2*(-B*F - E*C - H*L - K*I)
727 tmp.setPolyParams(0.0, 1.0);
728 tmp2.setPolyParams(u, t);
729 tmp.multiply(tmp2);
730 tmp.multiply(tmp2);
731 tmp.multiplyByScalar(-b * f - e * c - h * l - k * i);
732 result.add(tmp);
733
734 // (y*T + U)^2*(- F*C - L*I)
735 tmp.setPolyParams(u, t);
736 tmp2.setPolyParams(u, t);
737 tmp.multiply(tmp2);
738 tmp.multiplyByScalar(-f * c - l * i);
739 result.add(tmp);
740
741 return result;
742 }
743
744 /**
745 * Another of Kruppa's equations expressed as a polynomial of degree 4 to
746 * solve y value, which is the squared value of vertical focal length.
747 * This method is only used when aspect ratio is unknown.
748 *
749 * @param d internal value from Kruppa's equations.
750 * @param e internal value from Kruppa's equations.
751 * @param f internal value from Kruppa's equations.
752 * @param g internal value from Kruppa's equations.
753 * @param h internal value from Kruppa's equations.
754 * @param i internal value from Kruppa's equations.
755 * @param m internal value from Kruppa's equations.
756 * @param n internal value from Kruppa's equations.
757 * @param o internal value from Kruppa's equations.
758 * @param p internal value from Kruppa's equations.
759 * @param q internal value from Kruppa's equations.
760 * @param r internal value from Kruppa's equations.
761 * @param s internal value from Kruppa's equations.
762 * @param t internal value from Kruppa's equations.
763 * @param u internal value from Kruppa's equations.
764 * @param v internal value from Kruppa's equations.
765 * @param w internal value from Kruppa's equations.
766 * @return a polynomial.
767 */
768 private Polynomial buildPolynomial2(
769 final double d, final double e, final double f, final double g,
770 final double h, final double i, final double m, final double n,
771 final double o, final double p, final double q, final double r,
772 final double s, final double t, final double u, final double v,
773 final double w) {
774 // (-y^2*S - y*V - W)^2*(G*M - P*D) + y^2*(y*T + U)^2*(H*N - Q*E) +
775 // (-y^2*S - y*V - W)*(y*T + U)*y*(G*N + H*M - P*E - Q*D) +
776 // (-y^2*S - y*V - W)*(y*T + U)*(G*O + M*I - P*F - D*R) +
777 // y*(y*T + U)^2*(H*O + N*I - Q*F - E*R) + (y*T + U)^2*(O*I - F*R) = 0
778 final var result = new Polynomial(POLY_DEGREE_UNKNOWN_ASPECT_RATIO + 1);
779
780 // (-y^2*S - y*V - W)^2*(G*M - P*D)
781 final var tmp = new Polynomial(-w, -v, -s);
782 final var tmp2 = new Polynomial(-w, -v, -s);
783 tmp.multiply(tmp2);
784 tmp.multiplyByScalar(g * m - p * d);
785 result.add(tmp);
786
787 // y^2*(y*T + U)^2*(H*N - Q*E)
788 tmp.setPolyParams(0.0, 0.0, 1.0);
789 tmp2.setPolyParams(u, t);
790 tmp.multiply(tmp2);
791 tmp.multiply(tmp2);
792 tmp.multiplyByScalar(h * n - q * e);
793 result.add(tmp);
794
795 // (-y^2*S - y*V - W)*(y*T + U)*y*(G*N + H*M - P*E - Q*D)
796 tmp.setPolyParams(-w, -v, -s);
797 tmp2.setPolyParams(u, t);
798 tmp.multiply(tmp2);
799 tmp2.setPolyParams(0.0, 1.0);
800 tmp.multiply(tmp2);
801 tmp.multiplyByScalar(g * n + h * m - p * e - q * d);
802 result.add(tmp);
803
804 // (-y^2*S - y*V - W)*(y*T + U)*(G*O + M*I - P*F - D*R)
805 tmp.setPolyParams(-w, -v, -s);
806 tmp2.setPolyParams(u, t);
807 tmp.multiply(tmp2);
808 tmp.multiplyByScalar(g * o + m * i - p * f - d * r);
809 result.add(tmp);
810
811 // y*(y*T + U)^2*(H*O + N*I - Q*F - E*R)
812 tmp.setPolyParams(0.0, 1.0);
813 tmp2.setPolyParams(u, t);
814 tmp.multiply(tmp2);
815 tmp.multiplyByScalar(h * o + n * i - q * f - e * r);
816 result.add(tmp);
817
818 // (y*T + U)^2*(O*I - F*R)
819 tmp.setPolyParams(u, t);
820 tmp2.setPolyParams(u, t);
821 tmp.multiply(tmp2);
822 tmp.multiplyByScalar(o * i - f * r);
823 result.add(tmp);
824
825 return result;
826 }
827
828 /**
829 * Estimates the DIAC assuming known aspect ratio.
830 *
831 * @param result instance where estimated DIAC will be stored.
832 * @throws KruppaDualImageOfAbsoluteConicEstimatorException if an error
833 * occurs during estimation, usually because
834 * fundamental matrix corresponds to
835 * degenerate camera movements, or because of
836 * numerical instabilities.
837 */
838 private void estimateKnownAspectRatio(final DualImageOfAbsoluteConic result)
839 throws KruppaDualImageOfAbsoluteConicEstimatorException {
840 try {
841 final var x0 = principalPointX;
842 final var y0 = principalPointY;
843
844 // SVD decompose fundamental matrix
845 fundamentalMatrix.normalize();
846 final var decomposer = new SingularValueDecomposer(fundamentalMatrix.getInternalMatrix());
847 decomposer.decompose();
848
849 final var sigmas = decomposer.getSingularValues();
850 final var u = decomposer.getU();
851 final var v = decomposer.getV();
852
853 final var sigma1 = sigmas[0];
854 final var sigma2 = sigmas[1];
855
856 // Column u1
857 final var u11 = u.getElementAt(0, 0);
858 final var u21 = u.getElementAt(1, 0);
859 final var u31 = u.getElementAt(2, 0);
860
861 // Column u2
862 final var u12 = u.getElementAt(0, 1);
863 final var u22 = u.getElementAt(1, 1);
864 final var u32 = u.getElementAt(2, 1);
865
866 // Column v1
867 final var v11 = v.getElementAt(0, 0);
868 final var v21 = v.getElementAt(1, 0);
869 final var v31 = v.getElementAt(2, 0);
870
871 // Column v2
872 final var v12 = v.getElementAt(0, 1);
873 final var v22 = v.getElementAt(1, 1);
874 final var v32 = v.getElementAt(2, 1);
875
876 // build Kruppa equations
877 final var polyA = u12 * u11;
878 final var polyB = u22 * u21;
879 final var polyC = Math.pow(x0, 2.0) * u12 * u11 + x0 * y0 * u22 * u11 + x0 * u32 * u11
880 + x0 * y0 * u12 * u21 + Math.pow(y0, 2.0) * u22 * u21 + y0 * u32 * u21
881 + x0 * u12 * u31 + y0 * u22 * u31 + u32 * u31;
882 final var polyD = Math.pow(sigma2, 2.0) * v12 * v12;
883 final var polyE = Math.pow(sigma2, 2.0) * v22 * v22;
884 final var polyF = Math.pow(sigma2 * x0, 2.0) * v12 * v12
885 + Math.pow(sigma2, 2.0) * x0 * y0 * v22 * v12
886 + Math.pow(sigma2, 2.0) * x0 * v32 * v12
887 + Math.pow(sigma2, 2.0) * x0 * y0 * v12 * v22
888 + Math.pow(sigma2 * y0, 2.0) * v22 * v22
889 + Math.pow(sigma2, 2.0) * y0 * v32 * v22
890 + Math.pow(sigma2, 2.0) * x0 * v12 * v32
891 + Math.pow(sigma2, 2.0) * y0 * v22 * v32
892 + Math.pow(sigma2, 2.0) * v32 * v32;
893 final var polyG = u11 * u11;
894 final var polyH = u21 * u21;
895 final var polyI = Math.pow(x0, 2.0) * u11 * u11 + x0 * y0 * u21 * u11 + x0 * u31 * u11
896 + x0 * y0 * u11 * u21 + Math.pow(y0, 2.0) * u21 * u21 + y0 * u31 * u21
897 + x0 * u11 * u31 + y0 * u21 * u31 + u31 * u31;
898 final var polyJ = sigma1 * sigma2 * v12 * v11;
899 final var polyK = sigma1 * sigma2 * v22 * v21;
900 final var polyL = sigma1 * sigma2 * Math.pow(x0, 2.0) * v12 * v11
901 + sigma1 * sigma2 * x0 * y0 * v22 * v11 + sigma1 * sigma2 * x0 * v32 * v11
902 + sigma1 * sigma2 * x0 * y0 * v12 * v21
903 + sigma1 * sigma2 * Math.pow(y0, 2.0) * v22 * v21
904 + sigma1 * sigma2 * y0 * v32 * v21 + sigma1 * sigma2 * x0 * v12 * v31
905 + sigma1 * sigma2 * y0 * v22 * v31 + sigma1 * sigma2 * v32 * v31;
906 final var polyM = Math.pow(sigma1, 2.0) * v11 * v11;
907 final var polyN = Math.pow(sigma1, 2.0) * v21 * v21;
908 final var polyO = Math.pow(sigma1 * x0, 2.0) * v11 * v11
909 + Math.pow(sigma1, 2.0) * x0 * y0 * v21 * v11
910 + Math.pow(sigma1, 2.0) * x0 * v31 * v11
911 + Math.pow(sigma1, 2.0) * x0 * y0 * v11 * v21
912 + Math.pow(sigma1 * y0, 2.0) * v21 * v21
913 + Math.pow(sigma1, 2.0) * y0 * v31 * v21
914 + Math.pow(sigma1, 2.0) * x0 * v11 * v31
915 + Math.pow(sigma1, 2.0) * y0 * v21 * v31
916 + Math.pow(sigma1, 2.0) * v31 * v31;
917 final var polyP = u12 * u12;
918 final var polyQ = u22 * u22;
919 final var polyR = Math.pow(x0, 2.0) * u12 * u12 + x0 * y0 * u22 * u12 + x0 * u32 * u12
920 + x0 * y0 * u12 * u22 + Math.pow(y0, 2.0) * u22 * u22 + y0 * u32 * u22
921 + x0 * u12 * u32 + y0 * u22 * u32 + u32 * u32;
922
923 // try to solve any of Kruppa's equations
924 final var roots = knownAspectRatioRoots(polyA, polyB, polyC, polyD, polyE, polyF, polyG, polyH, polyI,
925 polyJ, polyK, polyL, polyM, polyN, polyO, polyP, polyQ, polyR);
926
927 // roots contain possible x values. We use only their real part
928
929 // pick the best x, y values that produce a positive definite DIAC
930 // matrix
931 final var diac = new DualImageOfAbsoluteConic(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
932 var valid = false;
933 if (roots != null) {
934 final var r = focalDistanceAspectRatio;
935 final var r2 = r * r;
936 for (final var root : roots) {
937 final var x = root.getReal();
938 final var y = r2 * x;
939
940 // build DIAC matrix and check if it is positive definite
941 if (x >= 0.0 && y >= 0.0) {
942 final var horizontalFocalLength = Math.sqrt(x);
943 final var verticalFocalLength = Math.sqrt(y);
944 valid = buildDiac(horizontalFocalLength, verticalFocalLength, diac);
945 }
946
947 if (valid) {
948 break;
949 }
950 }
951 } else {
952 throw new KruppaDualImageOfAbsoluteConicEstimatorException();
953 }
954
955 if (valid) {
956 // copy to result
957 result.setParameters(diac.getA(), diac.getB(), diac.getC(), diac.getD(), diac.getE(), diac.getF());
958 } else {
959 // no valid DIAC could be found
960 throw new KruppaDualImageOfAbsoluteConicEstimatorException();
961 }
962
963 } catch (final KruppaDualImageOfAbsoluteConicEstimatorException ex) {
964 throw ex;
965 } catch (final Exception ex) {
966 throw new KruppaDualImageOfAbsoluteConicEstimatorException(ex);
967 }
968 }
969
970 /**
971 * Another of Kruppa's equations expressed as a polynomial of degree 2 to
972 * solve x value, which is the squared value of horizontal focal length.
973 * This method is only used when aspect ratio is known.
974 *
975 * @param a internal value from Kruppa's equations.
976 * @param b internal value from Kruppa's equations.
977 * @param c internal value from Kruppa's equations.
978 * @param d internal value from Kruppa's equations.
979 * @param e internal value from Kruppa's equations.
980 * @param f internal value from Kruppa's equations.
981 * @param g internal value from Kruppa's equations.
982 * @param h internal value from Kruppa's equations.
983 * @param i internal value from Kruppa's equations.
984 * @param j internal value from Kruppa's equations.
985 * @param k internal value from Kruppa's equations.
986 * @param l internal value from Kruppa's equations.
987 * @return a polynomial.
988 */
989 private Polynomial buildPolynomial3(
990 final double a, final double b, final double c, final double d,
991 final double e, final double f, final double g, final double h,
992 final double i, final double j, final double k, final double l) {
993
994 // x^2*((-A*D - G*J) + r^4*(-B*E - H*K) + r^2*(-A*E - B*D - G*K - H*J)) +
995 // x*((-A*F - D*C - G*L - J*I) + r^2*(-B*F - E*C - H*L - K*I)) +
996 // (- F*C - L*I) = 0
997 final var r = focalDistanceAspectRatio;
998 final var r2 = r * r;
999 final var r4 = r2 * r2;
1000 return new Polynomial(-f * c - l * i,
1001 (-a * f - d * c - g * l - j * i) + r2 * (-b * f - e * c - h * l - k * i),
1002 ((-a * d - g * j) + r4 * (-b * e - h * k) + r2 * (-a * e - b * d - g * k - h * j)));
1003 }
1004
1005 /**
1006 * Another of Kruppa's equations expressed as a polynomial of degree 2 to
1007 * solve x value, which is the squared value of horizontal focal length.
1008 * This method is only used when aspect ratio is known.
1009 *
1010 * @param d internal value from Kruppa's equations.
1011 * @param e internal value from Kruppa's equations.
1012 * @param f internal value from Kruppa's equations.
1013 * @param g internal value from Kruppa's equations.
1014 * @param h internal value from Kruppa's equations.
1015 * @param i internal value from Kruppa's equations.
1016 * @param m internal value from Kruppa's equations.
1017 * @param n internal value from Kruppa's equations.
1018 * @param o internal value from Kruppa's equations.
1019 * @param p internal value from Kruppa's equations.
1020 * @param q internal value from Kruppa's equations.
1021 * @param r internal value from Kruppa's equations.
1022 * @return a polynomial.
1023 */
1024 private Polynomial buildPolynomial4(
1025 final double d, final double e, final double f, final double g,
1026 final double h, final double i, final double m, final double n,
1027 final double o, final double p, final double q, final double r) {
1028
1029 // x^2*((G*M - P*D) + r^4*(H*N - Q*E) + r^2*(G*N + H*M - P*E - Q*D)) +
1030 // x*((G*O + M*I - P*F - D*R) + r^2*(H*O + N*I - Q*F - E*R)) +
1031 // (O*I - F*R) = 0
1032 final var r1 = focalDistanceAspectRatio;
1033 final var r2 = r1 * r1;
1034 final var r4 = r2 * r2;
1035 return new Polynomial(o * i - f * r,
1036 (g * o + m * i - p * f - d * r) + r2 * (h * o + n * i - q * f - e * r),
1037 (g * m - p * d) + r4 * (h * n - q * e) + r2 * (g * n + h * m - p * e - q * d));
1038 }
1039
1040 /**
1041 * Another of Kruppa's equations expressed as a polynomial of degree 2 to
1042 * solve x value, which is the squared value of horizontal focal length.
1043 * This method is only used when aspect ratio is known.
1044 *
1045 * @param a internal value from Kruppa's equations.
1046 * @param b internal value from Kruppa's equations.
1047 * @param c internal value from Kruppa's equations.
1048 * @param j internal value from Kruppa's equations.
1049 * @param k internal value from Kruppa's equations.
1050 * @param l internal value from Kruppa's equations.
1051 * @param m internal value from Kruppa's equations.
1052 * @param n internal value from Kruppa's equations.
1053 * @param o internal value from Kruppa's equations.
1054 * @param p internal value from Kruppa's equations.
1055 * @param q internal value from Kruppa's equations.
1056 * @param r internal value from Kruppa's equations.
1057 * @return a polynomial.
1058 */
1059 private Polynomial buildPolynomial5(
1060 final double a, final double b, final double c, final double j,
1061 final double k, final double l, final double m, final double n,
1062 final double o, final double p, final double q, final double r) {
1063
1064 // x^2*((P*J + A*M) + r^4*(Q*K + B*N) + r^2*(P*K + Q*J + A*N + B*M)) +
1065 // x*((P*L + J*R + A*O + M*C) + r^2*(Q*L + K*R + B*O + N*C)) +
1066 // (L*R + O*C) = 0
1067 final var r1 = focalDistanceAspectRatio;
1068 final var r2 = r1 * r1;
1069 final var r4 = r2 * r2;
1070 return new Polynomial(l * r + o * c,
1071 (p * l + j * r + a * o + m * c) + r2 * (q * l + k * r + b * o + n * c),
1072 (p * j + a * m) + r4 * (q * k + b * n) + r2 * (p * k + q * j + a * n + b * m));
1073 }
1074
1075 /**
1076 * Solves Kruppa's equations when aspect ratio is unknown
1077 *
1078 * @param polyA A parameter of Kruppa's polynomial equation.
1079 * @param polyB B parameter of Kruppa's polynomial equation.
1080 * @param polyC C parameter of Kruppa's polynomial equation.
1081 * @param polyD D parameter of Kruppa's polynomial equation.
1082 * @param polyE E parameter of Kruppa's polynomial equation.
1083 * @param polyF F parameter of Kruppa's polynomial equation.
1084 * @param polyG G parameter of Kruppa's polynomial equation.
1085 * @param polyH H parameter of Kruppa's polynomial equation.
1086 * @param polyI I parameter of Kruppa's polynomial equation.
1087 * @param polyJ J parameter of Kruppa's polynomial equation.
1088 * @param polyK K parameter of Kruppa's polynomial equation.
1089 * @param polyL L parameter of Kruppa's polynomial equation.
1090 * @param polyM M parameter of Kruppa's polynomial equation.
1091 * @param polyN N parameter of Kruppa's polynomial equation.
1092 * @param polyO O parameter of Kruppa's polynomial equation.
1093 * @param polyP P parameter of Kruppa's polynomial equation.
1094 * @param polyQ Q parameter of Kruppa's polynomial equation.
1095 * @param polyR R parameter of Kruppa's polynomial equation.
1096 * @param polyS S parameter of Kruppa's polynomial equation.
1097 * @param polyT T parameter of Kruppa's polynomial equation.
1098 * @param polyU U parameter of Kruppa's polynomial equation.
1099 * @param polyV V parameter of Kruppa's polynomial equation.
1100 * @param polyW W parameter of Kruppa's polynomial equation.
1101 * @return roots solving Kruppa's equations.
1102 * @throws NumericalException if there are numerical instabilities.
1103 */
1104 private Complex[] unknownAspectRatioRoots(
1105 final double polyA, final double polyB, final double polyC, final double polyD,
1106 final double polyE, final double polyF, final double polyG, final double polyH,
1107 final double polyI, final double polyJ, final double polyK, final double polyL,
1108 final double polyM, final double polyN, final double polyO, final double polyP,
1109 final double polyQ, final double polyR, final double polyS, final double polyT,
1110 final double polyU, final double polyV, final double polyW) throws NumericalException {
1111 Complex[] roots;
1112 try {
1113 final var poly1 = buildPolynomial1(polyA, polyB, polyC, polyD, polyE, polyF, polyG,
1114 polyH, polyI, polyJ, polyK, polyL, polyS, polyT, polyU, polyV, polyW);
1115 roots = poly1.getRoots();
1116 } catch (final NumericalException ex1) {
1117 // if solution for poly1 fails, try with second polynomial
1118 final var poly2 = buildPolynomial2(polyD, polyE, polyF, polyG, polyH, polyI,
1119 polyM, polyN, polyO, polyP, polyQ, polyR, polyS, polyT, polyU, polyV, polyW);
1120 roots = poly2.getRoots();
1121 }
1122
1123 return roots;
1124 }
1125
1126 /**
1127 * Solves Kruppa's equations when aspect ratio is known
1128 *
1129 * @param polyA A parameter of Kruppa's polynomial equation.
1130 * @param polyB B parameter of Kruppa's polynomial equation.
1131 * @param polyC C parameter of Kruppa's polynomial equation.
1132 * @param polyD D parameter of Kruppa's polynomial equation.
1133 * @param polyE E parameter of Kruppa's polynomial equation.
1134 * @param polyF F parameter of Kruppa's polynomial equation.
1135 * @param polyG G parameter of Kruppa's polynomial equation.
1136 * @param polyH H parameter of Kruppa's polynomial equation.
1137 * @param polyI I parameter of Kruppa's polynomial equation.
1138 * @param polyJ J parameter of Kruppa's polynomial equation.
1139 * @param polyK K parameter of Kruppa's polynomial equation.
1140 * @param polyL L parameter of Kruppa's polynomial equation.
1141 * @param polyM M parameter of Kruppa's polynomial equation.
1142 * @param polyN N parameter of Kruppa's polynomial equation.
1143 * @param polyO O parameter of Kruppa's polynomial equation.
1144 * @param polyP P parameter of Kruppa's polynomial equation.
1145 * @param polyQ Q parameter of Kruppa's polynomial equation.
1146 * @param polyR R parameter of Kruppa's polynomial equation.
1147 * @return roots solving Kruppa's equations.
1148 * @throws NumericalException if there are numerical instabilities.
1149 */
1150 private Complex[] knownAspectRatioRoots(
1151 final double polyA, final double polyB, final double polyC, final double polyD,
1152 final double polyE, final double polyF, final double polyG, final double polyH,
1153 final double polyI, final double polyJ, final double polyK, final double polyL,
1154 final double polyM, final double polyN, final double polyO, final double polyP,
1155 final double polyQ, final double polyR) throws NumericalException {
1156 Complex[] roots;
1157 try {
1158 final var poly3 = buildPolynomial3(polyA, polyB, polyC, polyD, polyE, polyF, polyG,
1159 polyH, polyI, polyJ, polyK, polyL);
1160 roots = poly3.getRoots();
1161 } catch (final NumericalException e3) {
1162 try {
1163 // if solution for poly3 fails, try with 4th polynomial
1164 final var poly4 = buildPolynomial4(polyD, polyE, polyF, polyG, polyH,
1165 polyI, polyM, polyN, polyO, polyP, polyQ, polyR);
1166 roots = poly4.getRoots();
1167 } catch (final NumericalException e4) {
1168 final var poly5 = buildPolynomial5(polyA, polyB, polyC,
1169 polyJ, polyK, polyL, polyM, polyN, polyO, polyP, polyQ, polyR);
1170 roots = poly5.getRoots();
1171 }
1172 }
1173
1174 return roots;
1175 }
1176 }