+++ /dev/null
-///////////////////////////////////////////////////////////////////////////
-//
-// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
-// Digital Ltd. LLC
-//
-// All rights reserved.
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are
-// met:
-// * Redistributions of source code must retain the above copyright
-// notice, this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above
-// copyright notice, this list of conditions and the following disclaimer
-// in the documentation and/or other materials provided with the
-// distribution.
-// * Neither the name of Industrial Light & Magic nor the names of
-// its contributors may be used to endorse or promote products derived
-// from this software without specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
-// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
-// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
-// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-
-#ifndef INCLUDED_IMATHROOTS_H
-#define INCLUDED_IMATHROOTS_H
-
-//---------------------------------------------------------------------
-//
-// Functions to solve linear, quadratic or cubic equations
-//
-//---------------------------------------------------------------------
-
-#include <complex>
-
-namespace Imath {
-
-//--------------------------------------------------------------------------
-// Find the real solutions of a linear, quadratic or cubic equation:
-//
-// function equation solved
-//
-// solveLinear (a, b, x) a * x + b == 0
-// solveQuadratic (a, b, c, x) a * x*x + b * x + c == 0
-// solveNormalizedCubic (r, s, t, x) x*x*x + r * x*x + s * x + t == 0
-// solveCubic (a, b, c, d, x) a * x*x*x + b * x*x + c * x + d == 0
-//
-// Return value:
-//
-// 3 three real solutions, stored in x[0], x[1] and x[2]
-// 2 two real solutions, stored in x[0] and x[1]
-// 1 one real solution, stored in x[1]
-// 0 no real solutions
-// -1 all real numbers are solutions
-//
-// Notes:
-//
-// * It is possible that an equation has real solutions, but that the
-// solutions (or some intermediate result) are not representable.
-// In this case, either some of the solutions returned are invalid
-// (nan or infinity), or, if floating-point exceptions have been
-// enabled with Iex::mathExcOn(), an Iex::MathExc exception is
-// thrown.
-//
-// * Cubic equations are solved using Cardano's Formula; even though
-// only real solutions are produced, some intermediate results are
-// complex (std::complex<T>).
-//
-//--------------------------------------------------------------------------
-
-template <class T> int solveLinear (T a, T b, T &x);
-template <class T> int solveQuadratic (T a, T b, T c, T x[2]);
-template <class T> int solveNormalizedCubic (T r, T s, T t, T x[3]);
-template <class T> int solveCubic (T a, T b, T c, T d, T x[3]);
-
-
-//---------------
-// Implementation
-//---------------
-
-template <class T>
-int
-solveLinear (T a, T b, T &x)
-{
- if (a != 0)
- {
- x = -b / a;
- return 1;
- }
- else if (b != 0)
- {
- return 0;
- }
- else
- {
- return -1;
- }
-}
-
-
-template <class T>
-int
-solveQuadratic (T a, T b, T c, T x[2])
-{
- if (a == 0)
- {
- return solveLinear (b, c, x[0]);
- }
- else
- {
- T D = b * b - 4 * a * c;
-
- if (D > 0)
- {
- T s = sqrt (D);
-
- x[0] = (-b + s) / (2 * a);
- x[1] = (-b - s) / (2 * a);
- return 2;
- }
- if (D == 0)
- {
- x[0] = -b / (2 * a);
- return 1;
- }
- else
- {
- return 0;
- }
- }
-}
-
-
-template <class T>
-int
-solveNormalizedCubic (T r, T s, T t, T x[3])
-{
- T p = (3 * s - r * r) / 3;
- T q = 2 * r * r * r / 27 - r * s / 3 + t;
- T p3 = p / 3;
- T q2 = q / 2;
- T D = p3 * p3 * p3 + q2 * q2;
-
- if (D == 0 && p3 == 0)
- {
- x[0] = -r / 3;
- x[1] = -r / 3;
- x[2] = -r / 3;
- return 1;
- }
-
- std::complex<T> u = std::pow (-q / 2 + std::sqrt (std::complex<T> (D)),
- T (1) / T (3));
-
- std::complex<T> v = -p / (T (3) * u);
-
- const T sqrt3 = T (1.73205080756887729352744634150587); // enough digits
- // for long double
- std::complex<T> y0 (u + v);
-
- std::complex<T> y1 (-(u + v) / T (2) +
- (u - v) / T (2) * std::complex<T> (0, sqrt3));
-
- std::complex<T> y2 (-(u + v) / T (2) -
- (u - v) / T (2) * std::complex<T> (0, sqrt3));
-
- if (D > 0)
- {
- x[0] = y0.real() - r / 3;
- return 1;
- }
- else if (D == 0)
- {
- x[0] = y0.real() - r / 3;
- x[1] = y1.real() - r / 3;
- return 2;
- }
- else
- {
- x[0] = y0.real() - r / 3;
- x[1] = y1.real() - r / 3;
- x[2] = y2.real() - r / 3;
- return 3;
- }
-}
-
-
-template <class T>
-int
-solveCubic (T a, T b, T c, T d, T x[3])
-{
- if (a == 0)
- {
- return solveQuadratic (b, c, d, x);
- }
- else
- {
- return solveNormalizedCubic (b / a, c / a, d / a, x);
- }
-}
-
-
-} // namespace Imath
-
-#endif
projects / opencv / blobdiff