← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to imagery/i.ortho.photo/lib/georef.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.ortho.photo/lib/georef.c
 
1
#include <stdlib.h>
 
2
#include <signal.h>
 
3
#include "orthophoto.h"
 
4
 
 
5
/* TODO: replace with imagery lib I_compute_georef_equations() */
 
6
 
 
7
static int floating_exception;
 
8
static void catch(int);
 
9
static double determinant(double, double, double, double, double,
 
10
                          double, double, double, double);
 
11
 
 
12
/* find coefficients A,B,C for e2 = A + B*e1 + C*n1
 
13
 * also compute the reverse equations
 
14
 *
 
15
 * return 0 if no points
 
16
 *       -1 if not solvable
 
17
 *        1 if ok
 
18
 *
 
19
 * method is least squares.
 
20
 * the least squares problem reduces to solving the following
 
21
 * system of equations for A,B,C
 
22
 *
 
23
 *   s0*A + s1*B + s2*C = x0
 
24
 *   s1*A + s3*B + s4*C = x1
 
25
 *   s2*A + s4*B + s5*C = x2
 
26
 *
 
27
 * use Cramer's rule
 
28
 *
 
29
 *     | x0 s1 s2 |      | s0 x0 s2 |      | s0 s1 x0 |
 
30
 *     | x1 s3 s4 |      | s1 x1 s4 |      | s1 s3 x1 |
 
31
 *     | x2 s4 s5 |      | s2 x2 s5 |      | s2 s4 x2 |
 
32
 * A = ------------  B = ------------  C = ------------ 
 
33
 *     | s0 s1 s2 |      | s0 s1 s2 |      | s0 s1 s2 |
 
34
 *     | s1 s3 s4 |      | s1 s3 s4 |      | s1 s3 s4 |
 
35
 *     | s2 s4 s5 |      | s2 s4 s5 |      | s2 s4 s5 |
 
36
 *
 
37
 */
 
38
 
 
39
int
 
40
I_compute_ref_equations(struct Ortho_Photo_Points *cp, double E12[3],
 
41
                        double N12[3], double E21[3], double N21[3])
 
42
{
 
43
    double s0, s1, s2, s3, s4, s5;
 
44
    double x0, x1, x2;
 
45
    double det;
 
46
    void (*sigfpe) ();
 
47
    int i;
 
48
 
 
49
 
 
50
    s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
 
51
    for (i = 0; i < cp->count; i++) {
 
52
        if (cp->status[i] <= 0)
 
53
            continue;
 
54
        s0 += 1.0;
 
55
        s1 += cp->e1[i];
 
56
        s2 += cp->n1[i];
 
57
        s3 += cp->e1[i] * cp->e1[i];
 
58
        s4 += cp->e1[i] * cp->n1[i];
 
59
        s5 += cp->n1[i] * cp->n1[i];
 
60
    }
 
61
    if (s0 < 0.5)
 
62
        return 0;
 
63
 
 
64
    floating_exception = 0;
 
65
    sigfpe = signal(SIGFPE, catch);
 
66
 
 
67
    /* eastings */
 
68
    x0 = x1 = x2 = 0.0;
 
69
    for (i = 0; i < cp->count; i++) {
 
70
        if (cp->status[i] <= 0)
 
71
            continue;
 
72
        x0 += cp->e2[i];
 
73
        x1 += cp->e1[i] * cp->e2[i];
 
74
        x2 += cp->n1[i] * cp->e2[i];
 
75
    }
 
76
 
 
77
    det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
 
78
    if (det == 0.0) {
 
79
        signal(SIGFPE, sigfpe);
 
80
        return -1;
 
81
    }
 
82
    E12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
 
83
    E12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
 
84
    E12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
 
85
 
 
86
    /* northings */
 
87
    x0 = x1 = x2 = 0.0;
 
88
    for (i = 0; i < cp->count; i++) {
 
89
        if (cp->status[i] <= 0)
 
90
            continue;
 
91
        x0 += cp->n2[i];
 
92
        x1 += cp->e1[i] * cp->n2[i];
 
93
        x2 += cp->n1[i] * cp->n2[i];
 
94
    }
 
95
 
 
96
    det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
 
97
    if (det == 0.0) {
 
98
        signal(SIGFPE, sigfpe);
 
99
        return -1;
 
100
    }
 
101
    N12[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
 
102
    N12[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
 
103
    N12[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
 
104
 
 
105
    /* the inverse equations */
 
106
 
 
107
    s0 = s1 = s2 = s3 = s4 = s5 = 0.0;
 
108
    for (i = 0; i < cp->count; i++) {
 
109
        if (cp->status[i] <= 0)
 
110
            continue;
 
111
        s0 += 1.0;
 
112
        s1 += cp->e2[i];
 
113
        s2 += cp->n2[i];
 
114
        s3 += cp->e2[i] * cp->e2[i];
 
115
        s4 += cp->e2[i] * cp->n2[i];
 
116
        s5 += cp->n2[i] * cp->n2[i];
 
117
    }
 
118
 
 
119
    /* eastings */
 
120
    x0 = x1 = x2 = 0.0;
 
121
    for (i = 0; i < cp->count; i++) {
 
122
        if (cp->status[i] <= 0)
 
123
            continue;
 
124
        x0 += cp->e1[i];
 
125
        x1 += cp->e2[i] * cp->e1[i];
 
126
        x2 += cp->n2[i] * cp->e1[i];
 
127
    }
 
128
 
 
129
    det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
 
130
    if (det == 0.0) {
 
131
        signal(SIGFPE, sigfpe);
 
132
        return -1;
 
133
    }
 
134
    E21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
 
135
    E21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
 
136
    E21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
 
137
 
 
138
    /* northings */
 
139
    x0 = x1 = x2 = 0.0;
 
140
    for (i = 0; i < cp->count; i++) {
 
141
        if (cp->status[i] <= 0)
 
142
            continue;
 
143
        x0 += cp->n1[i];
 
144
        x1 += cp->e2[i] * cp->n1[i];
 
145
        x2 += cp->n2[i] * cp->n1[i];
 
146
    }
 
147
 
 
148
    det = determinant(s0, s1, s2, s1, s3, s4, s2, s4, s5);
 
149
    if (det == 0.0) {
 
150
        signal(SIGFPE, sigfpe);
 
151
        return -1;
 
152
    }
 
153
    N21[0] = determinant(x0, s1, s2, x1, s3, s4, x2, s4, s5) / det;
 
154
    N21[1] = determinant(s0, x0, s2, s1, x1, s4, s2, x2, s5) / det;
 
155
    N21[2] = determinant(s0, s1, x0, s1, s3, x1, s2, s4, x2) / det;
 
156
 
 
157
    signal(SIGFPE, sigfpe);
 
158
    return floating_exception ? -1 : 1;
 
159
}
 
160
 
 
161
static double determinant(double a, double b, double c, double d,
 
162
                          double e, double f, double g, double h, double i)
 
163
{
 
164
    /* compute determinant of 3x3 matrix
 
165
     *     | a b c |
 
166
     *     | d e f |
 
167
     *     | g h i |
 
168
     */
 
169
    return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
 
170
}
 
171
 
 
172
static void catch(int n)
 
173
{
 
174
    floating_exception = 1;
 
175
    signal(n, catch);
 
176
}