1 /**********************************************************************
2 * File: quspline.cpp (Formerly qspline.c)
3 * Description: Code for the QSPLINE class.
4 * Author: Ray Smith
5 * Created: Tue Oct 08 17:16:12 BST 1991
6 *
7 * (C) Copyright 1991, Hewlett-Packard Ltd.
8 ** Licensed under the Apache License, Version 2.0 (the "License");
9 ** you may not use this file except in compliance with the License.
10 ** You may obtain a copy of the License at
11 ** http://www.apache.org/licenses/LICENSE-2.0
12 ** Unless required by applicable law or agreed to in writing, software
13 ** distributed under the License is distributed on an "AS IS" BASIS,
14 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15 ** See the License for the specific language governing permissions and
16 ** limitations under the License.
17 *
18 **********************************************************************/
19
20 #include "mfcpch.h"
21 #include "memry.h"
22 #include "quadlsq.h"
23 #include "quspline.h"
24
25 #define QSPLINE_PRECISION 16 //no of steps to draw
26
27 /**********************************************************************
28 * QSPLINE::QSPLINE
29 *
30 * Constructor to build a QSPLINE given the components used in the old code.
31 **********************************************************************/
32
QSPLINE(inT32 count,inT32 * xstarts,double * coeffs)33 QSPLINE::QSPLINE( //constructor
34 inT32 count, //no of segments
35 inT32 *xstarts, //start coords
36 double *coeffs //coefficients
37 ) {
38 inT32 index; //segment index
39
40 //get memory
41 xcoords = (inT32 *) alloc_mem ((count + 1) * sizeof (inT32));
42 quadratics = (QUAD_COEFFS *) alloc_mem (count * sizeof (QUAD_COEFFS));
43 segments = count;
44 for (index = 0; index < segments; index++) {
45 //copy them
46 xcoords[index] = xstarts[index];
47 quadratics[index] = QUAD_COEFFS (coeffs[index * 3],
48 coeffs[index * 3 + 1],
49 coeffs[index * 3 + 2]);
50 }
51 //right edge
52 xcoords[index] = xstarts[index];
53 }
54
55
56 /**********************************************************************
57 * QSPLINE::QSPLINE
58 *
59 * Constructor to build a QSPLINE by appproximation of points.
60 **********************************************************************/
61
QSPLINE(int xstarts[],int segcount,int xpts[],int ypts[],int pointcount,int degree)62 QSPLINE::QSPLINE ( //constructor
63 int xstarts[], //spline boundaries
64 int segcount, //no of segments
65 int xpts[], //points to fit
66 int ypts[], int pointcount, //no of pts
67 int degree //fit required
68 ) {
69 register int pointindex; /*no along text line */
70 register int segment; /*segment no */
71 inT32 *ptcounts; //no in each segment
72 QLSQ qlsq; /*accumulator */
73
74 segments = segcount;
75 xcoords = (inT32 *) alloc_mem ((segcount + 1) * sizeof (inT32));
76 ptcounts = (inT32 *) alloc_mem ((segcount + 1) * sizeof (inT32));
77 quadratics = (QUAD_COEFFS *) alloc_mem (segcount * sizeof (QUAD_COEFFS));
78 memmove (xcoords, xstarts, (segcount + 1) * sizeof (inT32));
79 ptcounts[0] = 0; /*none in any yet */
80 for (segment = 0, pointindex = 0; pointindex < pointcount; pointindex++) {
81 while (segment < segcount && xpts[pointindex] >= xstarts[segment]) {
82 segment++; /*try next segment */
83 /*cumulative counts */
84 ptcounts[segment] = ptcounts[segment - 1];
85 }
86 ptcounts[segment]++; /*no in previous partition */
87 }
88 while (segment < segcount) {
89 segment++;
90 /*zero the rest */
91 ptcounts[segment] = ptcounts[segment - 1];
92 }
93
94 for (segment = 0; segment < segcount; segment++) {
95 qlsq.clear ();
96 /*first blob */
97 pointindex = ptcounts[segment];
98 if (pointindex > 0
99 && xpts[pointindex] != xpts[pointindex - 1]
100 && xpts[pointindex] != xstarts[segment])
101 qlsq.add (xstarts[segment],
102 ypts[pointindex - 1]
103 + (ypts[pointindex] - ypts[pointindex - 1])
104 * (xstarts[segment] - xpts[pointindex - 1])
105 / (xpts[pointindex] - xpts[pointindex - 1]));
106 for (; pointindex < ptcounts[segment + 1]; pointindex++) {
107 qlsq.add (xpts[pointindex], ypts[pointindex]);
108 }
109 if (pointindex > 0 && pointindex < pointcount
110 && xpts[pointindex] != xstarts[segment + 1])
111 qlsq.add (xstarts[segment + 1],
112 ypts[pointindex - 1]
113 + (ypts[pointindex] - ypts[pointindex - 1])
114 * (xstarts[segment + 1] - xpts[pointindex - 1])
115 / (xpts[pointindex] - xpts[pointindex - 1]));
116 qlsq.fit (degree);
117 quadratics[segment].a = qlsq.get_a ();
118 quadratics[segment].b = qlsq.get_b ();
119 quadratics[segment].c = qlsq.get_c ();
120 }
121 free_mem(ptcounts);
122 }
123
124
125 /**********************************************************************
126 * QSPLINE::QSPLINE
127 *
128 * Constructor to build a QSPLINE from another.
129 **********************************************************************/
130
QSPLINE(const QSPLINE & src)131 QSPLINE::QSPLINE( //constructor
132 const QSPLINE &src) {
133 segments = 0;
134 xcoords = NULL;
135 quadratics = NULL;
136 *this = src;
137 }
138
139
140 /**********************************************************************
141 * QSPLINE::~QSPLINE
142 *
143 * Destroy a QSPLINE.
144 **********************************************************************/
145
~QSPLINE()146 QSPLINE::~QSPLINE ( //constructor
147 ) {
148 if (xcoords != NULL) {
149 free_mem(xcoords);
150 xcoords = NULL;
151 }
152 if (quadratics != NULL) {
153 free_mem(quadratics);
154 quadratics = NULL;
155 }
156 }
157
158
159 /**********************************************************************
160 * QSPLINE::operator=
161 *
162 * Copy a QSPLINE
163 **********************************************************************/
164
operator =(const QSPLINE & source)165 QSPLINE & QSPLINE::operator= ( //assignment
166 const QSPLINE & source) {
167 if (xcoords != NULL)
168 free_mem(xcoords);
169 if (quadratics != NULL)
170 free_mem(quadratics);
171
172 segments = source.segments;
173 xcoords = (inT32 *) alloc_mem ((segments + 1) * sizeof (inT32));
174 quadratics = (QUAD_COEFFS *) alloc_mem (segments * sizeof (QUAD_COEFFS));
175 memmove (xcoords, source.xcoords, (segments + 1) * sizeof (inT32));
176 memmove (quadratics, source.quadratics, segments * sizeof (QUAD_COEFFS));
177 return *this;
178 }
179
180
181 /**********************************************************************
182 * QSPLINE::step
183 *
184 * Return the total of the step functions between the given coords.
185 **********************************************************************/
186
step(double x1,double x2)187 double QSPLINE::step( //find step functions
188 double x1, //between coords
189 double x2) {
190 int index1, index2; //indices of coords
191 double total; /*total steps */
192
193 index1 = spline_index (x1);
194 index2 = spline_index (x2);
195 total = 0;
196 while (index1 < index2) {
197 total +=
198 (double) quadratics[index1 + 1].y ((float) xcoords[index1 + 1]);
199 total -= (double) quadratics[index1].y ((float) xcoords[index1 + 1]);
200 index1++; /*next segment */
201 }
202 return total; /*total steps */
203 }
204
205
206 /**********************************************************************
207 * QSPLINE::y
208 *
209 * Return the y value at the given x value.
210 **********************************************************************/
211
y(double x) const212 double QSPLINE::y( //evaluate
213 double x //coord to evaluate at
214 ) const {
215 inT32 index; //segment index
216
217 index = spline_index (x);
218 return quadratics[index].y (x);//in correct segment
219 }
220
221
222 /**********************************************************************
223 * QSPLINE::spline_index
224 *
225 * Return the index to the largest xcoord not greater than x.
226 **********************************************************************/
227
spline_index(double x) const228 inT32 QSPLINE::spline_index( //evaluate
229 double x //coord to evaluate at
230 ) const {
231 inT32 index; //segment index
232 inT32 bottom; //bottom of range
233 inT32 top; //top of range
234
235 bottom = 0;
236 top = segments;
237 while (top - bottom > 1) {
238 index = (top + bottom) / 2; //centre of range
239 if (x >= xcoords[index])
240 bottom = index; //new min
241 else
242 top = index; //new max
243 }
244 return bottom;
245 }
246
247
248 /**********************************************************************
249 * QSPLINE::move
250 *
251 * Reposition spline by vector
252 **********************************************************************/
253
move(ICOORD vec)254 void QSPLINE::move( // reposition spline
255 ICOORD vec // by vector
256 ) {
257 inT32 segment; //index of segment
258 inT16 x_shift = vec.x ();
259
260 for (segment = 0; segment < segments; segment++) {
261 xcoords[segment] += x_shift;
262 quadratics[segment].move (vec);
263 }
264 xcoords[segment] += x_shift;
265 }
266
267
268 /**********************************************************************
269 * QSPLINE::overlap
270 *
271 * Return TRUE if spline2 overlaps this by no more than fraction less
272 * than the bounds of this.
273 **********************************************************************/
274
overlap(QSPLINE * spline2,double fraction)275 BOOL8 QSPLINE::overlap( //test overlap
276 QSPLINE *spline2, //2 cannot be smaller
277 double fraction //by more than this
278 ) {
279 int leftlimit; /*common left limit */
280 int rightlimit; /*common right limit */
281
282 leftlimit = xcoords[1];
283 rightlimit = xcoords[segments - 1];
284 /*or too non-overlap */
285 if (spline2->segments < 3 || spline2->xcoords[1] > leftlimit + fraction * (rightlimit - leftlimit)
286 || spline2->xcoords[spline2->segments - 1] < rightlimit
287 - fraction * (rightlimit - leftlimit))
288 return FALSE;
289 else
290 return TRUE;
291 }
292
293
294 /**********************************************************************
295 * extrapolate_spline
296 *
297 * Extrapolates the spline linearly using the same gradient as the
298 * quadratic has at either end.
299 **********************************************************************/
300
extrapolate(double gradient,int xmin,int xmax)301 void QSPLINE::extrapolate( //linear extrapolation
302 double gradient, //gradient to use
303 int xmin, //new left edge
304 int xmax //new right edge
305 ) {
306 register int segment; /*current segment of spline */
307 int dest_segment; //dest index
308 int *xstarts; //new boundaries
309 QUAD_COEFFS *quads; //new ones
310 int increment; //in size
311
312 increment = xmin < xcoords[0] ? 1 : 0;
313 if (xmax > xcoords[segments])
314 increment++;
315 if (increment == 0)
316 return;
317 xstarts = (int *) alloc_mem ((segments + 1 + increment) * sizeof (int));
318 quads =
319 (QUAD_COEFFS *) alloc_mem ((segments + increment) * sizeof (QUAD_COEFFS));
320 if (xmin < xcoords[0]) {
321 xstarts[0] = xmin;
322 quads[0].a = 0;
323 quads[0].b = gradient;
324 quads[0].c = y (xcoords[0]) - quads[0].b * xcoords[0];
325 dest_segment = 1;
326 }
327 else
328 dest_segment = 0;
329 for (segment = 0; segment < segments; segment++) {
330 xstarts[dest_segment] = xcoords[segment];
331 quads[dest_segment] = quadratics[segment];
332 dest_segment++;
333 }
334 xstarts[dest_segment] = xcoords[segment];
335 if (xmax > xcoords[segments]) {
336 quads[dest_segment].a = 0;
337 quads[dest_segment].b = gradient;
338 quads[dest_segment].c = y (xcoords[segments])
339 - quads[dest_segment].b * xcoords[segments];
340 dest_segment++;
341 xstarts[dest_segment] = xmax + 1;
342 }
343 segments = dest_segment;
344 free_mem(xcoords);
345 free_mem(quadratics);
346 xcoords = (inT32 *) xstarts;
347 quadratics = quads;
348 }
349
350
351 /**********************************************************************
352 * QSPLINE::plot
353 *
354 * Draw the QSPLINE in the given colour.
355 **********************************************************************/
356
357 #ifndef GRAPHICS_DISABLED
plot(ScrollView * window,ScrollView::Color colour) const358 void QSPLINE::plot( //draw it
359 ScrollView* window, //window to draw in
360 ScrollView::Color colour //colour to draw in
361 ) const {
362 inT32 segment; //index of segment
363 inT16 step; //index of poly piece
364 double increment; //x increment
365 double x; //x coord
366
367 window->Pen(colour);
368 for (segment = 0; segment < segments; segment++) {
369 increment =
370 (double) (xcoords[segment + 1] -
371 xcoords[segment]) / QSPLINE_PRECISION;
372 x = xcoords[segment];
373 for (step = 0; step <= QSPLINE_PRECISION; step++) {
374 if (segment == 0 && step == 0)
375 window->SetCursor(x, quadratics[segment].y (x));
376 else
377 window->DrawTo(x, quadratics[segment].y (x));
378 x += increment;
379 }
380 }
381 }
382 #endif
383