← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to imagery/i.rectify/crs.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.rectify/crs.c
1
 
 
2
 
/***************************************************************************/
3
 
 
4
 
/***************************************************************************/
5
 
/*
6
 
   CRS.C - Center for Remote Sensing rectification routines
7
 
 
8
 
   Written By: Brian J. Buckley
9
 
 
10
 
   At: The Center for Remote Sensing
11
 
   Michigan State University
12
 
   302 Berkey Hall
13
 
   East Lansing, MI  48824
14
 
   (517)353-7195
15
 
 
16
 
   Written: 12/19/91
17
 
 
18
 
   Last Update: 12/26/91 Brian J. Buckley
19
 
   Last Update:  1/24/92 Brian J. Buckley
20
 
   Added printout of trnfile. Triggered by BDEBUG.
21
 
   Last Update:  1/27/92 Brian J. Buckley
22
 
   Fixed bug so that only the active control points were used.
23
 
 */
24
 
 
25
 
/***************************************************************************/
26
 
 
27
 
/***************************************************************************/
28
 
 
29
 
#include <stdio.h>
30
 
#include <string.h>
31
 
#include <stdlib.h>
32
 
#include <math.h>
33
 
#include <limits.h>
34
 
#include <grass/gis.h>
35
 
 
36
 
/*
37
 
   #define MSDOS 1
38
 
 */
39
 
 
40
 
/*
41
 
   #define BDEBUG
42
 
 */
43
 
#ifdef BDEBUG
44
 
FILE *fp;
45
 
#endif
46
 
 
47
 
#ifdef MSDOS
48
 
 
49
 
#include "stdlib.h"
50
 
#include "dummy.h"
51
 
 
52
 
typedef double DOUBLE;
53
 
 
54
 
#else
55
 
 
56
 
#include <grass/imagery.h>
57
 
 
58
 
typedef double DOUBLE;
59
 
 
60
 
#endif
61
 
 
62
 
#include "crs.h"
63
 
 
64
 
/* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS.  THESE FUNCTIONS EXPECT
65
 
   SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
66
 
 
67
 
struct MATRIX
68
 
{
69
 
    int n;                      /* SIZE OF THIS MATRIX (N x N) */
70
 
    DOUBLE *v;
71
 
};
72
 
 
73
 
/* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
74
 
 
75
 
#define M(row,col) m->v[(((row)-1)*(m->n))+(col)-1]
76
 
 
77
 
/***************************************************************************/
78
 
/*
79
 
 */
80
 
 
81
 
/***************************************************************************/
82
 
 
83
 
#define MSUCCESS     1          /* SUCCESS */
84
 
#define MNPTERR      0          /* NOT ENOUGH POINTS */
85
 
#define MUNSOLVABLE -1          /* NOT SOLVABLE */
86
 
#define MMEMERR     -2          /* NOT ENOUGH MEMORY */
87
 
#define MPARMERR    -3          /* PARAMETER ERROR */
88
 
#define MINTERR     -4          /* INTERNAL ERROR */
89
 
 
90
 
/***************************************************************************/
91
 
/*
92
 
   FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
93
 
 */
94
 
 
95
 
/***************************************************************************/
96
 
 
97
 
#ifdef MSDOS
98
 
 
99
 
static int calccoef(struct Control_Points *, double *, double *, int);
100
 
static int calcls(struct Control_Points *, struct MATRIX *,
101
 
                  DOUBLE *, DOUBLE *, double *, double *);
102
 
static int exactdet(struct Control_Points *, struct MATRIX *,
103
 
                    DOUBLE *, DOUBLE *, double *, double *);
104
 
static int solvemat(struct MATRIX *, DOUBLE *, DOUBLE *, double *, double *);
105
 
static DOUBLE term(int, double, double);
106
 
 
107
 
#ifdef BDEBUG
108
 
static int checkgeoref(struct Control_Points *, double *, double *, int, int,
109
 
                       FILE * fp);
110
 
#endif
111
 
 
112
 
#else
113
 
 
114
 
static int calccoef();
115
 
static int calcls();
116
 
static int exactdet();
117
 
static int solvemat();
118
 
static DOUBLE term();
119
 
 
120
 
#ifdef BDEBUG
121
 
static int checkgeoref();
122
 
#endif
123
 
 
124
 
#endif
125
 
 
126
 
/***************************************************************************/
127
 
/*
128
 
   USE THIS TRANSFORMATION FUNCTION IF YOU WANT TO DO ARRAYS.
129
 
 */
130
 
 
131
 
/***************************************************************************/
132
 
 
133
 
#ifdef BETTERGEOREF
134
 
 
135
 
extern void CRS_georef(int, DOUBLE *, DOUBLE *, DOUBLE *, DOUBLE *, int);
136
 
 
137
 
void CRS_georef2(int order,     /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
138
 
                                   ORDER USED TO CALCULATE THE COEFFICIENTS */
139
 
                 double E[],    /* EASTING COEFFICIENTS */
140
 
                 double N[],    /* NORTHING COEFFICIENTS */
141
 
                 double e[],    /* EASTINGS TO BE TRANSFORMED */
142
 
                 double n[],    /* NORTHINGS TO BE TRANSFORMED */
143
 
                 int numpts     /* NUMBER OF POINTS TO BE TRANSFORMED */
144
 
    )
145
 
{
146
 
    DOUBLE e3, e2n, e2, en2, en, e1, n3, n2, n1;
147
 
    int i;
148
 
 
149
 
    switch (order) {
150
 
    case 3:
151
 
 
152
 
        for (i = 0; i < numpts; i++) {
153
 
            e1 = e[i];
154
 
            n1 = n[i];
155
 
            e2 = e1 * e1;
156
 
            en = e1 * n1;
157
 
            n2 = n1 * n1;
158
 
            e3 = e1 * e2;
159
 
            e2n = e2 * n1;
160
 
            en2 = e1 * n2;
161
 
            n3 = n1 * n2;
162
 
 
163
 
            e[i] = E[0] +
164
 
                E[1] * e1 + E[2] * n1 +
165
 
                E[3] * e2 + E[4] * en + E[5] * n2 +
166
 
                E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
167
 
            n[i] = N[0] +
168
 
                N[1] * e1 + N[2] * n1 +
169
 
                N[3] * e2 + N[4] * en + N[5] * n2 +
170
 
                N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
171
 
        }
172
 
        break;
173
 
 
174
 
    case 2:
175
 
 
176
 
        for (i = 0; i < numpts; i++) {
177
 
            e1 = e[i];
178
 
            n1 = n[i];
179
 
            e2 = e1 * e1;
180
 
            n2 = n1 * n1;
181
 
            en = e1 * n1;
182
 
 
183
 
            e[i] = E[0] + E[1] * e1 + E[2] * n1 +
184
 
                E[3] * e2 + E[4] * en + E[5] * n2;
185
 
            n[i] = N[0] + N[1] * e1 + N[2] * n1 +
186
 
                N[3] * e2 + N[4] * en + N[5] * n2;
187
 
        }
188
 
        break;
189
 
 
190
 
    case 1:
191
 
 
192
 
        for (i = 0; i < numpts; i++) {
193
 
            e1 = e[i];
194
 
            n1 = n[i];
195
 
            e[i] = E[0] + E[1] * e1 + E[2] * n1;
196
 
            n[i] = N[0] + N[1] * e1 + N[2] * n1;
197
 
        }
198
 
        break;
199
 
    }
200
 
}
201
 
 
202
 
#endif
203
 
 
204
 
/***************************************************************************/
205
 
/*
206
 
   TRANSFORM A SINGLE COORDINATE PAIR.
207
 
 */
208
 
 
209
 
/***************************************************************************/
210
 
 
211
 
int CRS_georef(double e1,       /* EASTINGS TO BE TRANSFORMED */
212
 
               double n1,       /* NORTHINGS TO BE TRANSFORMED */
213
 
               double *e,       /* EASTINGS TO BE TRANSFORMED */
214
 
               double *n,       /* NORTHINGS TO BE TRANSFORMED */
215
 
               double E[],      /* EASTING COEFFICIENTS */
216
 
               double N[],      /* NORTHING COEFFICIENTS */
217
 
               int order        /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
218
 
                                   ORDER USED TO CALCULATE THE COEFFICIENTS */
219
 
    )
220
 
{
221
 
    DOUBLE e3, e2n, en2, n3, e2, en, n2;
222
 
 
223
 
    switch (order) {
224
 
    case 1:
225
 
 
226
 
        *e = E[0] + E[1] * e1 + E[2] * n1;
227
 
        *n = N[0] + N[1] * e1 + N[2] * n1;
228
 
        break;
229
 
 
230
 
    case 2:
231
 
 
232
 
        e2 = e1 * e1;
233
 
        n2 = n1 * n1;
234
 
        en = e1 * n1;
235
 
 
236
 
        *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en + E[5] * n2;
237
 
        *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en + N[5] * n2;
238
 
        break;
239
 
 
240
 
    case 3:
241
 
 
242
 
        e2 = e1 * e1;
243
 
        en = e1 * n1;
244
 
        n2 = n1 * n1;
245
 
        e3 = e1 * e2;
246
 
        e2n = e2 * n1;
247
 
        en2 = e1 * n2;
248
 
        n3 = n1 * n2;
249
 
 
250
 
        *e = E[0] +
251
 
            E[1] * e1 + E[2] * n1 +
252
 
            E[3] * e2 + E[4] * en + E[5] * n2 +
253
 
            E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
254
 
        *n = N[0] +
255
 
            N[1] * e1 + N[2] * n1 +
256
 
            N[3] * e2 + N[4] * en + N[5] * n2 +
257
 
            N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
258
 
        break;
259
 
 
260
 
    default:
261
 
 
262
 
        return (MPARMERR);
263
 
        break;
264
 
    }
265
 
 
266
 
    return (MSUCCESS);
267
 
}
268
 
 
269
 
/***************************************************************************/
270
 
/*
271
 
   COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
272
 
 */
273
 
 
274
 
/***************************************************************************/
275
 
 
276
 
int
277
 
CRS_compute_georef_equations(struct Control_Points *cp, double E12[],
278
 
                             double N12[], double E21[], double N21[],
279
 
                             int order)
280
 
{
281
 
    double *tempptr;
282
 
    int status;
283
 
 
284
 
    if (order < 1 || order > MAXORDER)
285
 
        return (MPARMERR);
286
 
 
287
 
#ifdef BDEBUG
288
 
    fp = fopen("error.dat", "w");
289
 
    if (fp == NULL)
290
 
        return (-1);
291
 
#endif
292
 
 
293
 
    /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
294
 
 
295
 
    status = calccoef(cp, E12, N12, order);
296
 
    if (status != MSUCCESS)
297
 
        return (status);
298
 
 
299
 
#ifdef BDEBUG
300
 
    checkgeoref(cp, E12, N12, order, 1, fp);
301
 
#endif
302
 
 
303
 
    /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
304
 
 
305
 
    tempptr = cp->e1;
306
 
    cp->e1 = cp->e2;
307
 
    cp->e2 = tempptr;
308
 
    tempptr = cp->n1;
309
 
    cp->n1 = cp->n2;
310
 
    cp->n2 = tempptr;
311
 
 
312
 
    /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
313
 
 
314
 
    status = calccoef(cp, E21, N21, order);
315
 
 
316
 
#ifdef BDEBUG
317
 
    checkgeoref(cp, E21, N21, order, 0, fp);
318
 
    fclose(fp);
319
 
#endif
320
 
 
321
 
    /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
322
 
 
323
 
    tempptr = cp->e1;
324
 
    cp->e1 = cp->e2;
325
 
    cp->e2 = tempptr;
326
 
    tempptr = cp->n1;
327
 
    cp->n1 = cp->n2;
328
 
    cp->n2 = tempptr;
329
 
 
330
 
    return (status);
331
 
}
332
 
 
333
 
/***************************************************************************/
334
 
/*
335
 
   COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
336
 
 */
337
 
 
338
 
/***************************************************************************/
339
 
 
340
 
static int
341
 
calccoef(struct Control_Points *cp, double E[], double N[], int order)
342
 
{
343
 
    struct MATRIX m;
344
 
    DOUBLE *a;
345
 
    DOUBLE *b;
346
 
    int numactive;              /* NUMBER OF ACTIVE CONTROL POINTS */
347
 
    int status, i;
348
 
 
349
 
    /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
350
 
 
351
 
    for (i = numactive = 0; i < cp->count; i++) {
352
 
        if (cp->status[i] > 0)
353
 
            numactive++;
354
 
    }
355
 
 
356
 
    /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
357
 
       A TRANSFORMATION OF THIS ORDER */
358
 
 
359
 
    m.n = ((order + 1) * (order + 2)) / 2;
360
 
 
361
 
    if (numactive < m.n)
362
 
        return (MNPTERR);
363
 
 
364
 
    /* INITIALIZE MATRIX */
365
 
 
366
 
    m.v = (DOUBLE *) G_calloc(m.n * m.n, sizeof(DOUBLE));
367
 
    if (m.v == NULL) {
368
 
        return (MMEMERR);
369
 
    }
370
 
    a = (DOUBLE *) G_calloc(m.n, sizeof(DOUBLE));
371
 
    if (a == NULL) {
372
 
        G_free((char *)m.v);
373
 
        return (MMEMERR);
374
 
    }
375
 
    b = (DOUBLE *) G_calloc(m.n, sizeof(DOUBLE));
376
 
    if (b == NULL) {
377
 
        G_free((char *)m.v);
378
 
        G_free((char *)a);
379
 
        return (MMEMERR);
380
 
    }
381
 
 
382
 
    if (numactive == m.n)
383
 
        status = exactdet(cp, &m, a, b, E, N);
384
 
    else
385
 
        status = calcls(cp, &m, a, b, E, N);
386
 
 
387
 
    G_free((char *)m.v);
388
 
    G_free((char *)a);
389
 
    G_free((char *)b);
390
 
 
391
 
    return (status);
392
 
}
393
 
 
394
 
/***************************************************************************/
395
 
/*
396
 
   CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
397
 
   NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
398
 
 */
399
 
 
400
 
/***************************************************************************/
401
 
 
402
 
static int exactdet(struct Control_Points *cp, struct MATRIX *m, DOUBLE a[], DOUBLE b[], double E[],    /* EASTING COEFFICIENTS */
403
 
                    double N[]  /* NORTHING COEFFICIENTS */
404
 
    )
405
 
{
406
 
    int pntnow, currow, j;
407
 
 
408
 
    currow = 1;
409
 
    for (pntnow = 0; pntnow < cp->count; pntnow++) {
410
 
        if (cp->status[pntnow] > 0) {
411
 
            /* POPULATE MATRIX M */
412
 
 
413
 
#ifdef BDEBUG
414
 
            fprintf(fp, "%2d ", pntnow + 1);
415
 
#endif
416
 
 
417
 
            for (j = 1; j <= m->n; j++) {
418
 
                M(currow, j) = term(j, cp->e1[pntnow], cp->n1[pntnow]);
419
 
#ifdef BDEBUG
420
 
                fprintf(fp, "%+14.7le ", M(currow, j));
421
 
                if (j == 5)
422
 
                    fprintf(fp, "\n   ");
423
 
#endif
424
 
            }
425
 
#ifdef BDEBUG
426
 
            fprintf(fp, "\n");
427
 
#endif
428
 
 
429
 
            /* POPULATE MATRIX A AND B */
430
 
 
431
 
            a[currow - 1] = cp->e2[pntnow];
432
 
            b[currow - 1] = cp->n2[pntnow];
433
 
#ifdef BDEBUG
434
 
            fprintf(fp, "   %+14.7le ", a[currow - 1]);
435
 
            fprintf(fp, "%+14.7le\n", b[currow - 1]);
436
 
#endif
437
 
 
438
 
            currow++;
439
 
        }
440
 
#ifdef BDEBUG
441
 
        else {
442
 
            fprintf(fp, "%2d UNUSED\n", pntnow + 1);
443
 
        }
444
 
#endif
445
 
    }
446
 
 
447
 
    if (currow - 1 != m->n)
448
 
        return (MINTERR);
449
 
 
450
 
    return (solvemat(m, a, b, E, N));
451
 
}
452
 
 
453
 
/***************************************************************************/
454
 
/*
455
 
   CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
456
 
   NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.  THIS
457
 
   ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
458
 
 */
459
 
 
460
 
/***************************************************************************/
461
 
 
462
 
static int calcls(struct Control_Points *cp, struct MATRIX *m, DOUBLE a[], DOUBLE b[], double E[],      /* EASTING COEFFICIENTS */
463
 
                  double N[]    /* NORTHING COEFFICIENTS */
464
 
    )
465
 
{
466
 
    int i, j, n, numactive = 0;
467
 
 
468
 
    /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
469
 
 
470
 
    for (i = 1; i <= m->n; i++) {
471
 
        for (j = i; j <= m->n; j++)
472
 
            M(i, j) = 0.0;
473
 
        a[i - 1] = b[i - 1] = 0.0;
474
 
    }
475
 
 
476
 
    /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
477
 
       THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
478
 
 
479
 
    for (n = 0; n < cp->count; n++) {
480
 
        if (cp->status[n] > 0) {
481
 
            numactive++;
482
 
            for (i = 1; i <= m->n; i++) {
483
 
                for (j = i; j <= m->n; j++)
484
 
                    M(i, j) +=
485
 
                        term(i, cp->e1[n], cp->n1[n]) * term(j, cp->e1[n],
486
 
                                                             cp->n1[n]);
487
 
 
488
 
                a[i - 1] += cp->e2[n] * term(i, cp->e1[n], cp->n1[n]);
489
 
                b[i - 1] += cp->n2[n] * term(i, cp->e1[n], cp->n1[n]);
490
 
            }
491
 
        }
492
 
    }
493
 
 
494
 
    if (numactive <= m->n)
495
 
        return (MINTERR);
496
 
 
497
 
    /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
498
 
 
499
 
    for (i = 2; i <= m->n; i++) {
500
 
        for (j = 1; j < i; j++)
501
 
            M(i, j) = M(j, i);
502
 
    }
503
 
 
504
 
    return (solvemat(m, a, b, E, N));
505
 
}
506
 
 
507
 
/***************************************************************************/
508
 
/*
509
 
   CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
510
 
 
511
 
   ORDER\TERM   1    2    3    4    5    6    7    8    9   10
512
 
   1        e0n0 e1n0 e0n1
513
 
   2        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
514
 
   3        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
515
 
 */
516
 
 
517
 
/***************************************************************************/
518
 
 
519
 
static DOUBLE term(int term, double e, double n)
520
 
{
521
 
    switch (term) {
522
 
    case 1:
523
 
        return ((DOUBLE) 1.0);
524
 
    case 2:
525
 
        return ((DOUBLE) e);
526
 
    case 3:
527
 
        return ((DOUBLE) n);
528
 
    case 4:
529
 
        return ((DOUBLE) (e * e));
530
 
    case 5:
531
 
        return ((DOUBLE) (e * n));
532
 
    case 6:
533
 
        return ((DOUBLE) (n * n));
534
 
    case 7:
535
 
        return ((DOUBLE) (e * e * e));
536
 
    case 8:
537
 
        return ((DOUBLE) (e * e * n));
538
 
    case 9:
539
 
        return ((DOUBLE) (e * n * n));
540
 
    case 10:
541
 
        return ((DOUBLE) (n * n * n));
542
 
    }
543
 
    return ((DOUBLE) 0.0);
544
 
}
545
 
 
546
 
/***************************************************************************/
547
 
/*
548
 
   SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
549
 
   GAUSSIAN ELIMINATION METHOD.
550
 
 
551
 
   | M11 M12 ... M1n | | E0   |   | a0   |
552
 
   | M21 M22 ... M2n | | E1   | = | a1   |
553
 
   |  .   .   .   .  | | .    |   | .    |
554
 
   | Mn1 Mn2 ... Mnn | | En-1 |   | an-1 |
555
 
 
556
 
   and
557
 
 
558
 
   | M11 M12 ... M1n | | N0   |   | b0   |
559
 
   | M21 M22 ... M2n | | N1   | = | b1   |
560
 
   |  .   .   .   .  | | .    |   | .    |
561
 
   | Mn1 Mn2 ... Mnn | | Nn-1 |   | bn-1 |
562
 
 */
563
 
 
564
 
/***************************************************************************/
565
 
 
566
 
static int solvemat(struct MATRIX *m,
567
 
                    DOUBLE a[], DOUBLE b[], double E[], double N[])
568
 
{
569
 
    int i, j, i2, j2, imark;
570
 
    DOUBLE factor, temp;
571
 
    DOUBLE pivot;               /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
572
 
 
573
 
    for (i = 1; i <= m->n; i++) {
574
 
        j = i;
575
 
 
576
 
        /* find row with largest magnitude value for pivot value */
577
 
 
578
 
        pivot = M(i, j);
579
 
        imark = i;
580
 
        for (i2 = i + 1; i2 <= m->n; i2++) {
581
 
            temp = fabs(M(i2, j));
582
 
            if (temp > fabs(pivot)) {
583
 
                pivot = M(i2, j);
584
 
                imark = i2;
585
 
            }
586
 
        }
587
 
 
588
 
        /* if the pivot is very small then the points are nearly co-linear */
589
 
        /* co-linear points result in an undefined matrix, and nearly */
590
 
        /* co-linear points results in a solution with rounding error */
591
 
 
592
 
        if (pivot == 0.0)
593
 
            return (MUNSOLVABLE);
594
 
 
595
 
        /* if row with highest pivot is not the current row, switch them */
596
 
 
597
 
        if (imark != i) {
598
 
            for (j2 = 1; j2 <= m->n; j2++) {
599
 
                temp = M(imark, j2);
600
 
                M(imark, j2) = M(i, j2);
601
 
                M(i, j2) = temp;
602
 
            }
603
 
 
604
 
            temp = a[imark - 1];
605
 
            a[imark - 1] = a[i - 1];
606
 
            a[i - 1] = temp;
607
 
 
608
 
            temp = b[imark - 1];
609
 
            b[imark - 1] = b[i - 1];
610
 
            b[i - 1] = temp;
611
 
        }
612
 
 
613
 
        /* compute zeros above and below the pivot, and compute
614
 
           values for the rest of the row as well */
615
 
 
616
 
        for (i2 = 1; i2 <= m->n; i2++) {
617
 
            if (i2 != i) {
618
 
                factor = M(i2, j) / pivot;
619
 
                for (j2 = j; j2 <= m->n; j2++)
620
 
                    M(i2, j2) -= factor * M(i, j2);
621
 
                a[i2 - 1] -= factor * a[i - 1];
622
 
                b[i2 - 1] -= factor * b[i - 1];
623
 
            }
624
 
        }
625
 
    }
626
 
 
627
 
    /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
628
 
       COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
629
 
 
630
 
    for (i = 1; i <= m->n; i++) {
631
 
        E[i - 1] = a[i - 1] / M(i, i);
632
 
        N[i - 1] = b[i - 1] / M(i, i);
633
 
    }
634
 
 
635
 
    return (MSUCCESS);
636
 
}
637
 
 
638
 
/***************************************************************************/
639
 
/*
640
 
 */
641
 
 
642
 
/***************************************************************************/
643
 
 
644
 
#ifdef BDEBUG
645
 
 
646
 
static int checkgeoref(struct Control_Points *cp,
647
 
                       double E[], double N[], int order, int forward,
648
 
                       FILE * fp)
649
 
{
650
 
    DOUBLE xrms, yrms, dx, dy, dx2, dy2, totaldist, dist;
651
 
    double tempx, tempy;
652
 
    int i, n, numactive;
653
 
 
654
 
    n = ((order + 1) * (order + 2)) / 2;
655
 
 
656
 
    if (forward)
657
 
        fprintf(fp, "FORWARD:\n");
658
 
    else
659
 
        fprintf(fp, "BACKWARD:\n");
660
 
 
661
 
    fprintf(fp, "%d order\n", order);
662
 
    for (i = 0; i < n; i++)
663
 
        fprintf(fp, "%+.17E %+.17E\n", E[i], N[i]);
664
 
 
665
 
    xrms = yrms = dx2 = dy2 = totaldist = 0.0;
666
 
    numactive = 0;
667
 
    for (i = 0; i < cp->count; i++) {
668
 
        fprintf(fp, "\nCONTROL POINT: %d\n", i + 1);
669
 
 
670
 
        fprintf(fp, "%20s: %+.20lE %+.20lE\n", "ORIGINAL POINT",
671
 
                cp->e1[i], cp->n1[i]);
672
 
        fprintf(fp, "%20s: %+.20lE %+.20lE\n", "DESIRED POINT",
673
 
                cp->e2[i], cp->n2[i]);
674
 
 
675
 
        if (cp->status[i] > 0) {
676
 
            numactive++;
677
 
            CRS_georef(cp->e1[i], cp->n1[i], &tempx, &tempy, E, N, order);
678
 
 
679
 
            fprintf(fp, "%20s: %+.20lE %+.20lE\n", "CALCULATED POINT", tempx,
680
 
                    tempy);
681
 
            dx = tempx - cp->e2[i];
682
 
            dy = tempy - cp->n2[i];
683
 
            fprintf(fp, "%20s: %+.20lE %+.20lE\n", "RESIDUAL ERROR", dx, dy);
684
 
            dx2 = dx * dx;
685
 
            dy2 = dy * dy;
686
 
            dist = sqrt(dx2 + dy2);
687
 
            fprintf(fp, "%20s: %+.20lE\n", "DISTANCE (RMS) ERROR", dist);
688
 
 
689
 
            xrms += dx2;
690
 
            yrms += dy2;
691
 
 
692
 
            totaldist += dist;
693
 
        }
694
 
        else
695
 
            fprintf(fp, "NOT USED\n");
696
 
    }
697
 
    xrms = sqrt(xrms / (DOUBLE) numactive);
698
 
    yrms = sqrt(yrms / (DOUBLE) numactive);
699
 
 
700
 
    fprintf(fp, "\n%20s: %+.20lE %+.20lE\n", "RMS ERROR", xrms, yrms);
701
 
 
702
 
    fprintf(fp, "\n%20s: %+.20lE\n", "TOTAL RMS ERROR",
703
 
            sqrt(xrms * xrms + yrms * yrms));
704
 
 
705
 
    fprintf(fp, "\n%20s: %+.20lE\n", "AVG. DISTANCE ERROR",
706
 
            totaldist / numactive);
707
 
 
708
 
    return (0);
709
 
}
710
 
 
711
 
#endif