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 | |