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
23dng_spline_solver::dng_spline_solver ()
24
25 : X ()
26 , Y ()
27 , S ()
28
29 {
30
31 }
32
33/******************************************************************************/
34
35dng_spline_solver::~dng_spline_solver ()
36 {
37
38 }
39
40/******************************************************************************/
41
42void dng_spline_solver::Reset ()
43 {
44
45 X.clear ();
46 Y.clear ();
47
48 S.clear ();
49
50 }
51
52/******************************************************************************/
53
54void 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
64void 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
156bool 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
174real64 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