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