| /*****************************************************************************/ |
| // Copyright 2006-2007 Adobe Systems Incorporated |
| // All Rights Reserved. |
| // |
| // NOTICE: Adobe permits you to use, modify, and distribute this file in |
| // accordance with the terms of the Adobe license agreement accompanying it. |
| /*****************************************************************************/ |
| |
| /* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_spline.cpp#1 $ */ |
| /* $DateTime: 2012/05/30 13:28:51 $ */ |
| /* $Change: 832332 $ */ |
| /* $Author: tknoll $ */ |
| |
| /*****************************************************************************/ |
| |
| #include "dng_spline.h" |
| |
| #include "dng_assertions.h" |
| #include "dng_exceptions.h" |
| |
| /******************************************************************************/ |
| |
| dng_spline_solver::dng_spline_solver () |
| |
| : X () |
| , Y () |
| , S () |
| |
| { |
| |
| } |
| |
| /******************************************************************************/ |
| |
| dng_spline_solver::~dng_spline_solver () |
| { |
| |
| } |
| |
| /******************************************************************************/ |
| |
| void dng_spline_solver::Reset () |
| { |
| |
| X.clear (); |
| Y.clear (); |
| |
| S.clear (); |
| |
| } |
| |
| /******************************************************************************/ |
| |
| void dng_spline_solver::Add (real64 x, real64 y) |
| { |
| |
| X.push_back (x); |
| Y.push_back (y); |
| |
| } |
| |
| /******************************************************************************/ |
| |
| void dng_spline_solver::Solve () |
| { |
| |
| // This code computes the unique curve such that: |
| // It is C0, C1, and C2 continuous |
| // The second derivative is zero at the end points |
| |
| int32 count = (int32) X.size (); |
| |
| DNG_ASSERT (count >= 2, "Too few points"); |
| |
| int32 start = 0; |
| int32 end = count; |
| |
| real64 A = X [start+1] - X [start]; |
| real64 B = (Y [start+1] - Y [start]) / A; |
| |
| S.resize (count); |
| |
| S [start] = B; |
| |
| int32 j; |
| |
| // Slopes here are a weighted average of the slopes |
| // to each of the adjcent control points. |
| |
| for (j = start + 2; j < end; ++j) |
| { |
| |
| real64 C = X [j] - X [j-1]; |
| real64 D = (Y [j] - Y [j-1]) / C; |
| |
| S [j-1] = (B * C + D * A) / (A + C); |
| |
| A = C; |
| B = D; |
| |
| } |
| |
| S [end-1] = 2.0 * B - S [end-2]; |
| S [start] = 2.0 * S [start] - S [start+1]; |
| |
| if ((end - start) > 2) |
| { |
| |
| dng_std_vector<real64> E; |
| dng_std_vector<real64> F; |
| dng_std_vector<real64> G; |
| |
| E.resize (count); |
| F.resize (count); |
| G.resize (count); |
| |
| F [start] = 0.5; |
| E [end-1] = 0.5; |
| G [start] = 0.75 * (S [start] + S [start+1]); |
| G [end-1] = 0.75 * (S [end-2] + S [end-1]); |
| |
| for (j = start+1; j < end - 1; ++j) |
| { |
| |
| A = (X [j+1] - X [j-1]) * 2.0; |
| |
| E [j] = (X [j+1] - X [j]) / A; |
| F [j] = (X [j] - X [j-1]) / A; |
| G [j] = 1.5 * S [j]; |
| |
| } |
| |
| for (j = start+1; j < end; ++j) |
| { |
| |
| A = 1.0 - F [j-1] * E [j]; |
| |
| if (j != end-1) F [j] /= A; |
| |
| G [j] = (G [j] - G [j-1] * E [j]) / A; |
| |
| } |
| |
| for (j = end - 2; j >= start; --j) |
| G [j] = G [j] - F [j] * G [j+1]; |
| |
| for (j = start; j < end; ++j) |
| S [j] = G [j]; |
| |
| } |
| |
| } |
| |
| /******************************************************************************/ |
| |
| bool dng_spline_solver::IsIdentity () const |
| { |
| |
| int32 count = (int32) X.size (); |
| |
| if (count != 2) |
| return false; |
| |
| if (X [0] != 0.0 || X [1] != 1.0 || |
| Y [0] != 0.0 || Y [1] != 1.0) |
| return false; |
| |
| return true; |
| |
| } |
| |
| /******************************************************************************/ |
| |
| real64 dng_spline_solver::Evaluate (real64 x) const |
| { |
| |
| int32 count = (int32) X.size (); |
| |
| // Check for off each end of point list. |
| |
| if (x <= X [0]) |
| return Y [0]; |
| |
| if (x >= X [count-1]) |
| return Y [count-1]; |
| |
| // Binary search for the index. |
| |
| int32 lower = 1; |
| int32 upper = count - 1; |
| |
| while (upper > lower) |
| { |
| |
| int32 mid = (lower + upper) >> 1; |
| |
| if (x == X [mid]) |
| { |
| return Y [mid]; |
| } |
| |
| if (x > X [mid]) |
| lower = mid + 1; |
| else |
| upper = mid; |
| |
| } |
| |
| DNG_ASSERT (upper == lower, "Binary search error in point list"); |
| |
| int32 j = lower; |
| |
| // X [j - 1] < x <= X [j] |
| // A is the distance between the X [j] and X [j - 1] |
| // B and C describe the fractional distance to either side. B + C = 1. |
| |
| // We compute a cubic spline between the two points with slopes |
| // S[j-1] and S[j] at either end. Specifically, we compute the 1-D Bezier |
| // with control values: |
| // |
| // Y[j-1], Y[j-1] + S[j-1]*A, Y[j]-S[j]*A, Y[j] |
| |
| return EvaluateSplineSegment (x, |
| X [j - 1], |
| Y [j - 1], |
| S [j - 1], |
| X [j ], |
| Y [j ], |
| S [j ]); |
| |
| } |
| |
| /*****************************************************************************/ |