~paparazzi-uav/paparazzi/v5.0-manual

« back to all changes in this revision

Viewing changes to sw/ext/opencv_bebop/opencv/3rdparty/openexr/Imath/ImathRoots.h

  • Committer: Paparazzi buildbot
  • Date: 2016-05-18 15:00:29 UTC
  • Revision ID: felix.ruess+docbot@gmail.com-20160518150029-e8lgzi5kvb4p7un9
Manual import commit 4b8bbb730080dac23cf816b98908dacfabe2a8ec from v5.0 branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
///////////////////////////////////////////////////////////////////////////
 
2
//
 
3
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
 
4
// Digital Ltd. LLC
 
5
//
 
6
// All rights reserved.
 
7
//
 
8
// Redistribution and use in source and binary forms, with or without
 
9
// modification, are permitted provided that the following conditions are
 
10
// met:
 
11
// *       Redistributions of source code must retain the above copyright
 
12
// notice, this list of conditions and the following disclaimer.
 
13
// *       Redistributions in binary form must reproduce the above
 
14
// copyright notice, this list of conditions and the following disclaimer
 
15
// in the documentation and/or other materials provided with the
 
16
// distribution.
 
17
// *       Neither the name of Industrial Light & Magic nor the names of
 
18
// its contributors may be used to endorse or promote products derived
 
19
// from this software without specific prior written permission.
 
20
//
 
21
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 
22
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 
23
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 
24
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 
25
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 
26
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 
27
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 
28
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 
29
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 
30
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 
31
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
32
//
 
33
///////////////////////////////////////////////////////////////////////////
 
34
 
 
35
 
 
36
 
 
37
#ifndef INCLUDED_IMATHROOTS_H
 
38
#define INCLUDED_IMATHROOTS_H
 
39
 
 
40
//---------------------------------------------------------------------
 
41
//
 
42
//      Functions to solve linear, quadratic or cubic equations
 
43
//
 
44
//---------------------------------------------------------------------
 
45
 
 
46
#include <ImathMath.h>
 
47
#include <complex>
 
48
 
 
49
namespace Imath {
 
50
 
 
51
//--------------------------------------------------------------------------
 
52
// Find the real solutions of a linear, quadratic or cubic equation:
 
53
//
 
54
//      function                                   equation solved
 
55
//
 
56
//   solveLinear (a, b, x)                                    a * x + b == 0
 
57
//   solveQuadratic (a, b, c, x)                    a * x*x + b * x + c == 0
 
58
//   solveNormalizedCubic (r, s, t, x)      x*x*x + r * x*x + s * x + t == 0
 
59
//   solveCubic (a, b, c, d, x)         a * x*x*x + b * x*x + c * x + d == 0
 
60
//
 
61
// Return value:
 
62
//
 
63
//       3      three real solutions, stored in x[0], x[1] and x[2]
 
64
//       2      two real solutions, stored in x[0] and x[1]
 
65
//       1      one real solution, stored in x[1]
 
66
//       0      no real solutions
 
67
//      -1      all real numbers are solutions
 
68
//
 
69
// Notes:
 
70
//
 
71
//    * It is possible that an equation has real solutions, but that the
 
72
//      solutions (or some intermediate result) are not representable.
 
73
//      In this case, either some of the solutions returned are invalid
 
74
//      (nan or infinity), or, if floating-point exceptions have been
 
75
//      enabled with Iex::mathExcOn(), an Iex::MathExc exception is
 
76
//      thrown.
 
77
//
 
78
//    * Cubic equations are solved using Cardano's Formula; even though
 
79
//      only real solutions are produced, some intermediate results are
 
80
//      complex (std::complex<T>).
 
81
//
 
82
//--------------------------------------------------------------------------
 
83
 
 
84
template <class T> int  solveLinear (T a, T b, T &x);
 
85
template <class T> int  solveQuadratic (T a, T b, T c, T x[2]);
 
86
template <class T> int  solveNormalizedCubic (T r, T s, T t, T x[3]);
 
87
template <class T> int  solveCubic (T a, T b, T c, T d, T x[3]);
 
88
 
 
89
 
 
90
//---------------
 
91
// Implementation
 
92
//---------------
 
93
 
 
94
template <class T>
 
95
int
 
96
solveLinear (T a, T b, T &x)
 
97
{
 
98
    if (a != 0)
 
99
    {
 
100
    x = -b / a;
 
101
    return 1;
 
102
    }
 
103
    else if (b != 0)
 
104
    {
 
105
    return 0;
 
106
    }
 
107
    else
 
108
    {
 
109
    return -1;
 
110
    }
 
111
}
 
112
 
 
113
 
 
114
template <class T>
 
115
int
 
116
solveQuadratic (T a, T b, T c, T x[2])
 
117
{
 
118
    if (a == 0)
 
119
    {
 
120
    return solveLinear (b, c, x[0]);
 
121
    }
 
122
    else
 
123
    {
 
124
    T D = b * b - 4 * a * c;
 
125
 
 
126
    if (D > 0)
 
127
    {
 
128
        T s = Math<T>::sqrt (D);
 
129
        T q = -(b + (b > 0 ? 1 : -1) * s) / T(2);
 
130
 
 
131
        x[0] = q / a;
 
132
        x[1] = c / q;
 
133
        return 2;
 
134
    }
 
135
    if (D == 0)
 
136
    {
 
137
        x[0] = -b / (2 * a);
 
138
        return 1;
 
139
    }
 
140
    else
 
141
    {
 
142
        return 0;
 
143
    }
 
144
    }
 
145
}
 
146
 
 
147
 
 
148
template <class T>
 
149
int
 
150
solveNormalizedCubic (T r, T s, T t, T x[3])
 
151
{
 
152
    T p  = (3 * s - r * r) / 3;
 
153
    T q  = 2 * r * r * r / 27 - r * s / 3 + t;
 
154
    T p3 = p / 3;
 
155
    T q2 = q / 2;
 
156
    T D  = p3 * p3 * p3 + q2 * q2;
 
157
 
 
158
    if (D == 0 && p3 == 0)
 
159
    {
 
160
    x[0] = -r / 3;
 
161
    x[1] = -r / 3;
 
162
    x[2] = -r / 3;
 
163
    return 1;
 
164
    }
 
165
 
 
166
    std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
 
167
                  T (1) / T (3));
 
168
 
 
169
    std::complex<T> v = -p / (T (3) * u);
 
170
 
 
171
    const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
 
172
                                // for long double
 
173
    std::complex<T> y0 (u + v);
 
174
 
 
175
    std::complex<T> y1 (-(u + v) / T (2) +
 
176
             (u - v) / T (2) * std::complex<T> (0, sqrt3));
 
177
 
 
178
    std::complex<T> y2 (-(u + v) / T (2) -
 
179
             (u - v) / T (2) * std::complex<T> (0, sqrt3));
 
180
 
 
181
    if (D > 0)
 
182
    {
 
183
    x[0] = y0.real() - r / 3;
 
184
    return 1;
 
185
    }
 
186
    else if (D == 0)
 
187
    {
 
188
    x[0] = y0.real() - r / 3;
 
189
    x[1] = y1.real() - r / 3;
 
190
    return 2;
 
191
    }
 
192
    else
 
193
    {
 
194
    x[0] = y0.real() - r / 3;
 
195
    x[1] = y1.real() - r / 3;
 
196
    x[2] = y2.real() - r / 3;
 
197
    return 3;
 
198
    }
 
199
}
 
200
 
 
201
 
 
202
template <class T>
 
203
int
 
204
solveCubic (T a, T b, T c, T d, T x[3])
 
205
{
 
206
    if (a == 0)
 
207
    {
 
208
    return solveQuadratic (b, c, d, x);
 
209
    }
 
210
    else
 
211
    {
 
212
    return solveNormalizedCubic (b / a, c / a, d / a, x);
 
213
    }
 
214
}
 
215
 
 
216
 
 
217
} // namespace Imath
 
218
 
 
219
#endif