1 /*****************************************************************************/ 2 // Copyright 2006-2007 Adobe Systems Incorporated 3 // All Rights Reserved. 4 // 5 // NOTICE: Adobe permits you to use, modify, and distribute this file in 6 // accordance with the terms of the Adobe license agreement accompanying it. 7 /*****************************************************************************/ 8 9 /* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_spline.cpp#1 $ */ 10 /* $DateTime: 2012/05/30 13:28:51 $ */ 11 /* $Change: 832332 $ */ 12 /* $Author: tknoll $ */ 13 14 /*****************************************************************************/ 15 16 #include "dng_spline.h" 17 18 #include "dng_assertions.h" 19 #include "dng_exceptions.h" 20 21 /******************************************************************************/ 22 23 dng_spline_solver::dng_spline_solver () 24 25 : X () 26 , Y () 27 , S () 28 29 { 30 31 } 32 33 /******************************************************************************/ 34 35 dng_spline_solver::~dng_spline_solver () 36 { 37 38 } 39 40 /******************************************************************************/ 41 42 void dng_spline_solver::Reset () 43 { 44 45 X.clear (); 46 Y.clear (); 47 48 S.clear (); 49 50 } 51 52 /******************************************************************************/ 53 54 void dng_spline_solver::Add (real64 x, real64 y) 55 { 56 57 X.push_back (x); 58 Y.push_back (y); 59 60 } 61 62 /******************************************************************************/ 63 64 void dng_spline_solver::Solve () 65 { 66 67 // This code computes the unique curve such that: 68 // It is C0, C1, and C2 continuous 69 // The second derivative is zero at the end points 70 71 int32 count = (int32) X.size (); 72 73 DNG_ASSERT (count >= 2, "Too few points"); 74 75 int32 start = 0; 76 int32 end = count; 77 78 real64 A = X [start+1] - X [start]; 79 real64 B = (Y [start+1] - Y [start]) / A; 80 81 S.resize (count); 82 83 S [start] = B; 84 85 int32 j; 86 87 // Slopes here are a weighted average of the slopes 88 // to each of the adjcent control points. 89 90 for (j = start + 2; j < end; ++j) 91 { 92 93 real64 C = X [j] - X [j-1]; 94 real64 D = (Y [j] - Y [j-1]) / C; 95 96 S [j-1] = (B * C + D * A) / (A + C); 97 98 A = C; 99 B = D; 100 101 } 102 103 S [end-1] = 2.0 * B - S [end-2]; 104 S [start] = 2.0 * S [start] - S [start+1]; 105 106 if ((end - start) > 2) 107 { 108 109 dng_std_vector<real64> E; 110 dng_std_vector<real64> F; 111 dng_std_vector<real64> G; 112 113 E.resize (count); 114 F.resize (count); 115 G.resize (count); 116 117 F [start] = 0.5; 118 E [end-1] = 0.5; 119 G [start] = 0.75 * (S [start] + S [start+1]); 120 G [end-1] = 0.75 * (S [end-2] + S [end-1]); 121 122 for (j = start+1; j < end - 1; ++j) 123 { 124 125 A = (X [j+1] - X [j-1]) * 2.0; 126 127 E [j] = (X [j+1] - X [j]) / A; 128 F [j] = (X [j] - X [j-1]) / A; 129 G [j] = 1.5 * S [j]; 130 131 } 132 133 for (j = start+1; j < end; ++j) 134 { 135 136 A = 1.0 - F [j-1] * E [j]; 137 138 if (j != end-1) F [j] /= A; 139 140 G [j] = (G [j] - G [j-1] * E [j]) / A; 141 142 } 143 144 for (j = end - 2; j >= start; --j) 145 G [j] = G [j] - F [j] * G [j+1]; 146 147 for (j = start; j < end; ++j) 148 S [j] = G [j]; 149 150 } 151 152 } 153 154 /******************************************************************************/ 155 156 bool dng_spline_solver::IsIdentity () const 157 { 158 159 int32 count = (int32) X.size (); 160 161 if (count != 2) 162 return false; 163 164 if (X [0] != 0.0 || X [1] != 1.0 || 165 Y [0] != 0.0 || Y [1] != 1.0) 166 return false; 167 168 return true; 169 170 } 171 172 /******************************************************************************/ 173 174 real64 dng_spline_solver::Evaluate (real64 x) const 175 { 176 177 int32 count = (int32) X.size (); 178 179 // Check for off each end of point list. 180 181 if (x <= X [0]) 182 return Y [0]; 183 184 if (x >= X [count-1]) 185 return Y [count-1]; 186 187 // Binary search for the index. 188 189 int32 lower = 1; 190 int32 upper = count - 1; 191 192 while (upper > lower) 193 { 194 195 int32 mid = (lower + upper) >> 1; 196 197 if (x == X [mid]) 198 { 199 return Y [mid]; 200 } 201 202 if (x > X [mid]) 203 lower = mid + 1; 204 else 205 upper = mid; 206 207 } 208 209 DNG_ASSERT (upper == lower, "Binary search error in point list"); 210 211 int32 j = lower; 212 213 // X [j - 1] < x <= X [j] 214 // A is the distance between the X [j] and X [j - 1] 215 // B and C describe the fractional distance to either side. B + C = 1. 216 217 // We compute a cubic spline between the two points with slopes 218 // S[j-1] and S[j] at either end. Specifically, we compute the 1-D Bezier 219 // with control values: 220 // 221 // Y[j-1], Y[j-1] + S[j-1]*A, Y[j]-S[j]*A, Y[j] 222 223 return EvaluateSplineSegment (x, 224 X [j - 1], 225 Y [j - 1], 226 S [j - 1], 227 X [j ], 228 Y [j ], 229 S [j ]); 230 231 } 232 233 /*****************************************************************************/ 234
