1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2016 Imperial College London
5 * Copyright 2016 Andreas Schuh
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/PolynomialSolvers.h"
21
22 #include "mirtk/Math.h"
23
24
25 namespace mirtk {
26
27
28 // =============================================================================
29 // Root finding
30 // =============================================================================
31
32 // -----------------------------------------------------------------------------
SolveCubicEquation(double x[3],double a,double b,double c)33 int SolveCubicEquation(double x[3], double a, double b, double c)
34 {
35 double a2 = a*a;
36 double q = (a2 - 3.*b)/9.;
37 double r = (a*(2.*a2 - 9.*b) + 27.*c) / 54.;
38 double r2 = r*r;
39 double q3 = q*q*q;
40 double A, B;
41 if(r2<q3) {
42 double t = acos(clamp(r / sqrt(q3), -1., 1.));
43 a /= 3.;
44 q = -2. * sqrt(q);
45 x[0] = q * cos( t / 3.) - a;
46 x[1] = q * cos((t + two_pi) / 3.) - a;
47 x[2] = q * cos((t - two_pi) / 3.) - a;
48 return 3;
49 } else {
50 A = -pow(abs(r) + sqrt(r2 - q3), 1./3.);
51 if (r < 0.) A = -A;
52 B = (A == 0. ? 0 : B = q / A);
53 a /= 3.;
54 x[0] = (A + B) - a;
55 x[1] = -.5 * (A + B) -a;
56 x[2] = .5 * sqrt(3.) * (A - B);
57 if (abs(x[2]) < 1e-14) {
58 x[2] = x[1];
59 return 2;
60 }
61 return 1;
62 }
63 }
64
65 // =============================================================================
66 // Minimization
67 // =============================================================================
68
69 // -----------------------------------------------------------------------------
MinimumOf4thDegreePolynomial(double a,double b,double c)70 double MinimumOf4thDegreePolynomial(double a, double b, double c)
71 {
72 int i, j, n;
73 double x[3], x2;
74
75 // Find roots of derivative
76 n = SolveCubicEquation(x, 4., 3. * a, 2. * b, c);
77
78 // Determine global minimum
79 double value, min_value = numeric_limits<double>::max();
80 for (i = j = 0; i < n; ++i) {
81 x2 = x[i] * x[i];
82 value = x2 * x2 + a * x2 * x[i] + b * x2 + c * x[i];
83 if (value < min_value) {
84 min_value = value, j = i;
85 }
86 }
87
88 return x[j];
89 }
90
91
92 } // namespace mirtk
93