00001
00020 #include "mfcpch.h"
00021 #include "memry.h"
00022 #include "quadlsq.h"
00023 #include "quspline.h"
00024
00025 #define QSPLINE_PRECISION 16 //no of steps to draw
00026
00030 QSPLINE::QSPLINE(
00031 INT32 count,
00032 INT32 *xstarts,
00033 double *coeffs
00034 ) {
00035 INT32 index;
00036
00037
00038 xcoords = (INT32 *) alloc_mem ((count + 1) * sizeof (INT32));
00039 quadratics = (QUAD_COEFFS *) alloc_mem (count * sizeof (QUAD_COEFFS));
00040 segments = count;
00041 for (index = 0; index < segments; index++) {
00042
00043 xcoords[index] = xstarts[index];
00044 quadratics[index] = QUAD_COEFFS (coeffs[index * 3],
00045 coeffs[index * 3 + 1],
00046 coeffs[index * 3 + 2]);
00047 }
00048
00049 xcoords[index] = xstarts[index];
00050 }
00051
00052
00056 QSPLINE::QSPLINE (
00057 int xstarts[],
00058 int segcount,
00059 int xpts[],
00060 int ypts[], int pointcount,
00061 int degree
00062 ) {
00063 register int pointindex;
00064 register int segment;
00065 INT32 *ptcounts;
00066 QLSQ qlsq;
00067
00068 segments = segcount;
00069 xcoords = (INT32 *) alloc_mem ((segcount + 1) * sizeof (INT32));
00070 ptcounts = (INT32 *) alloc_mem ((segcount + 1) * sizeof (INT32));
00071 quadratics = (QUAD_COEFFS *) alloc_mem (segcount * sizeof (QUAD_COEFFS));
00072 memmove (xcoords, xstarts, (segcount + 1) * sizeof (INT32));
00073 ptcounts[0] = 0;
00074 for (segment = 0, pointindex = 0; pointindex < pointcount; pointindex++) {
00075 while (segment < segcount && xpts[pointindex] >= xstarts[segment]) {
00076 segment++;
00077
00078 ptcounts[segment] = ptcounts[segment - 1];
00079 }
00080 ptcounts[segment]++;
00081 }
00082 while (segment < segcount) {
00083 segment++;
00084
00085 ptcounts[segment] = ptcounts[segment - 1];
00086 }
00087
00088 for (segment = 0; segment < segcount; segment++) {
00089 qlsq.clear ();
00090
00091 pointindex = ptcounts[segment];
00092 if (pointindex > 0
00093 && xpts[pointindex] != xpts[pointindex - 1]
00094 && xpts[pointindex] != xstarts[segment])
00095 qlsq.add (xstarts[segment],
00096 ypts[pointindex - 1]
00097 + (ypts[pointindex] - ypts[pointindex - 1])
00098 * (xstarts[segment] - xpts[pointindex - 1])
00099 / (xpts[pointindex] - xpts[pointindex - 1]));
00100 for (; pointindex < ptcounts[segment + 1]; pointindex++) {
00101 qlsq.add (xpts[pointindex], ypts[pointindex]);
00102 }
00103 if (pointindex > 0 && pointindex < pointcount
00104 && xpts[pointindex] != xstarts[segment + 1])
00105 qlsq.add (xstarts[segment + 1],
00106 ypts[pointindex - 1]
00107 + (ypts[pointindex] - ypts[pointindex - 1])
00108 * (xstarts[segment + 1] - xpts[pointindex - 1])
00109 / (xpts[pointindex] - xpts[pointindex - 1]));
00110 qlsq.fit (degree);
00111 quadratics[segment].a = qlsq.get_a ();
00112 quadratics[segment].b = qlsq.get_b ();
00113 quadratics[segment].c = qlsq.get_c ();
00114 }
00115 free_mem(ptcounts);
00116 }
00117
00118
00122 QSPLINE::QSPLINE(
00123 const QSPLINE &src) {
00124 segments = 0;
00125 xcoords = NULL;
00126 quadratics = NULL;
00127 *this = src;
00128 }
00129
00130
00134 QSPLINE::~QSPLINE (
00135 ) {
00136 if (xcoords != NULL) {
00137 free_mem(xcoords);
00138 xcoords = NULL;
00139 }
00140 if (quadratics != NULL) {
00141 free_mem(quadratics);
00142 quadratics = NULL;
00143 }
00144 }
00145
00146
00150 QSPLINE & QSPLINE::operator= (
00151 const QSPLINE & source) {
00152 if (xcoords != NULL)
00153 free_mem(xcoords);
00154 if (quadratics != NULL)
00155 free_mem(quadratics);
00156
00157 segments = source.segments;
00158 xcoords = (INT32 *) alloc_mem ((segments + 1) * sizeof (INT32));
00159 quadratics = (QUAD_COEFFS *) alloc_mem (segments * sizeof (QUAD_COEFFS));
00160 memmove (xcoords, source.xcoords, (segments + 1) * sizeof (INT32));
00161 memmove (quadratics, source.quadratics, segments * sizeof (QUAD_COEFFS));
00162 return *this;
00163 }
00164
00165
00169 double QSPLINE::step(
00170 double x1,
00171 double x2) {
00172 int index1, index2;
00173 double total;
00174
00175 index1 = spline_index (x1);
00176 index2 = spline_index (x2);
00177 total = 0;
00178 while (index1 < index2) {
00179 total +=
00180 (double) quadratics[index1 + 1].y ((float) xcoords[index1 + 1]);
00181 total -= (double) quadratics[index1].y ((float) xcoords[index1 + 1]);
00182 index1++;
00183 }
00184 return total;
00185 }
00186
00187
00191 double QSPLINE::y(
00192 double x
00193 ) const {
00194 INT32 index;
00195
00196 index = spline_index (x);
00197 return quadratics[index].y (x);
00198 }
00199
00200
00204 INT32 QSPLINE::spline_index(
00205 double x
00206 ) const {
00207 INT32 index;
00208 INT32 bottom;
00209 INT32 top;
00210
00211 bottom = 0;
00212 top = segments;
00213 while (top - bottom > 1) {
00214 index = (top + bottom) / 2;
00215 if (x >= xcoords[index])
00216 bottom = index;
00217 else
00218 top = index;
00219 }
00220 return bottom;
00221 }
00222
00223
00227 void QSPLINE::move(
00228 ICOORD vec
00229 ) {
00230 INT32 segment;
00231 INT16 x_shift = vec.x ();
00232
00233 for (segment = 0; segment < segments; segment++) {
00234 xcoords[segment] += x_shift;
00235 quadratics[segment].move (vec);
00236 }
00237 xcoords[segment] += x_shift;
00238 }
00239
00240
00245 BOOL8 QSPLINE::overlap(
00246 QSPLINE *spline2,
00247 double fraction
00248 ) {
00249 int leftlimit;
00250 int rightlimit;
00251
00252 leftlimit = xcoords[1];
00253 rightlimit = xcoords[segments - 1];
00254
00255 if (spline2->segments < 3 || spline2->xcoords[1] > leftlimit + fraction * (rightlimit - leftlimit)
00256 || spline2->xcoords[spline2->segments - 1] < rightlimit
00257 - fraction * (rightlimit - leftlimit))
00258 return FALSE;
00259 else
00260 return TRUE;
00261 }
00262
00263
00268 void QSPLINE::extrapolate(
00269 double gradient,
00270 int xmin,
00271 int xmax
00272 ) {
00273 register int segment;
00274 int dest_segment;
00275 int *xstarts;
00276 QUAD_COEFFS *quads;
00277 int increment;
00278
00279 increment = xmin < xcoords[0] ? 1 : 0;
00280 if (xmax > xcoords[segments])
00281 increment++;
00282 if (increment == 0)
00283 return;
00284 xstarts = (int *) alloc_mem ((segments + 1 + increment) * sizeof (int));
00285 quads =
00286 (QUAD_COEFFS *) alloc_mem ((segments + increment) * sizeof (QUAD_COEFFS));
00287 if (xmin < xcoords[0]) {
00288 xstarts[0] = xmin;
00289 quads[0].a = 0;
00290 quads[0].b = gradient;
00291 quads[0].c = y (xcoords[0]) - quads[0].b * xcoords[0];
00292 dest_segment = 1;
00293 }
00294 else
00295 dest_segment = 0;
00296 for (segment = 0; segment < segments; segment++) {
00297 xstarts[dest_segment] = xcoords[segment];
00298 quads[dest_segment] = quadratics[segment];
00299 dest_segment++;
00300 }
00301 xstarts[dest_segment] = xcoords[segment];
00302 if (xmax > xcoords[segments]) {
00303 quads[dest_segment].a = 0;
00304 quads[dest_segment].b = gradient;
00305 quads[dest_segment].c = y (xcoords[segments])
00306 - quads[dest_segment].b * xcoords[segments];
00307 dest_segment++;
00308 xstarts[dest_segment] = xmax + 1;
00309 }
00310 segments = dest_segment;
00311 free_mem(xcoords);
00312 free_mem(quadratics);
00313 xcoords = (INT32 *) xstarts;
00314 quadratics = quads;
00315 }
00316
00317
00321 #ifndef GRAPHICS_DISABLED
00322 void QSPLINE::plot(
00323 WINDOW window,
00324 COLOUR colour
00325 ) const {
00326 INT32 segment;
00327 INT16 step;
00328 double increment;
00329 double x;
00330
00331 line_color_index(window, colour);
00332 for (segment = 0; segment < segments; segment++) {
00333 increment =
00334 (double) (xcoords[segment + 1] -
00335 xcoords[segment]) / QSPLINE_PRECISION;
00336 x = xcoords[segment];
00337 for (step = 0; step <= QSPLINE_PRECISION; step++) {
00338 if (segment == 0 && step == 0)
00339 move2d (window, x, quadratics[segment].y (x));
00340 else
00341 draw2d (window, x, quadratics[segment].y (x));
00342 x += increment;
00343 }
00344 }
00345 }
00346 #endif