← Back to branch summary

~ubuntu-branches/ubuntu/hardy/qgis/hardy

« back to all changes in this revision

Viewing changes to plugins/georeferencer/qgsleastsquares.cpp

  • Committer: Bazaar Package Importer
  • Author(s): William Grant
  • Date: 2007-05-06 13:42:32 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20070506134232-pyli6t388w5asd8x
Tags: 0.8.0-3ubuntu1
* Merge from Debian unstable. Remaining Ubuntu changes:
  - debian/rules, debian/qgis.install, debian/qgis.dirs debian/qgis.desktop:
    Add and install .desktop.
* debian/qgis.desktop: Remove Applications category; it's not real.
* Modify Maintainer value to match Debian-Maintainer-Field Spec

Show diffs side-by-side

added added

removed removed

Lines of Context:
plugins/georeferencer/qgsleastsquares.cpp
1
 
#include <cmath>
2
 
#include <stdexcept>
3
 
 
4
 
#include <gsl/gsl_linalg.h>
5
 
 
6
 
#include "qgsleastsquares.h"
7
 
 
8
 
 
9
 
void QgsLeastSquares::linear(std::vector<QgsPoint> mapCoords, 
10
 
                             std::vector<QgsPoint> pixelCoords,
11
 
                             QgsPoint& origin, double& pixelSize) {
12
 
  int n = mapCoords.size();
13
 
  if (n < 2) {
14
 
    throw std::domain_error("Fit to a linear transform requires at "
15
 
                            "least 2 points.");
16
 
  }
17
 
 
18
 
  double sumPx(0), sumPy(0), sumPx2(0), sumPy2(0), sumPxMx(0), sumPyMy(0),
19
 
    sumMx(0), sumMy(0);
20
 
  for (int i = 0; i < n; ++i) {
21
 
    sumPx += pixelCoords[i].x();
22
 
    sumPy += pixelCoords[i].y();
23
 
    sumPx2 += std::pow(pixelCoords[i].x(), 2);
24
 
    sumPy2 += std::pow(pixelCoords[i].y(), 2);
25
 
    sumPxMx += pixelCoords[i].x() * mapCoords[i].x();
26
 
    sumPyMy += pixelCoords[i].y() * mapCoords[i].y();
27
 
    sumMx += mapCoords[i].x();
28
 
    sumMy += mapCoords[i].y();
29
 
  }
30
 
  
31
 
  double deltaX = n * sumPx2 - std::pow(sumPx, 2);
32
 
  double deltaY = n * sumPy2 - std::pow(sumPy, 2);
33
 
  
34
 
  double aX = (sumPx2 * sumMx - sumPx * sumPxMx) / deltaX;
35
 
  double aY = (sumPy2 * sumMy - sumPy * sumPyMy) / deltaY;
36
 
  double bX = (n * sumPxMx - sumPx * sumMx) / deltaX;
37
 
  double bY = (n * sumPyMy - sumPy * sumMy) / deltaY;
38
 
  
39
 
  origin.setX(aX);
40
 
  origin.setY(aY);
41
 
  pixelSize = (std::abs(bX) + std::abs(bY)) / 2;
42
 
}
43
 
 
44
 
  
45
 
void QgsLeastSquares::helmert(std::vector<QgsPoint> mapCoords, 
46
 
                              std::vector<QgsPoint> pixelCoords,
47
 
                              QgsPoint& origin, double& pixelSize, 
48
 
                              double& rotation) {
49
 
  int n = mapCoords.size();
50
 
  if (n < 2) {
51
 
    throw std::domain_error("Fit to a Helmert transform requires at "
52
 
                            "least 2 points.");
53
 
  }
54
 
  
55
 
  double A = 0, B = 0, C = 0, D = 0, E = 0, F = 0, G = 0, H = 0, I = 0, J = 0;
56
 
  for (int i = 0; i < n; ++i) {
57
 
    A += pixelCoords[i].x();
58
 
    B += pixelCoords[i].y();
59
 
    C += mapCoords[i].x();
60
 
    D += mapCoords[i].y();
61
 
    E += mapCoords[i].x() * pixelCoords[i].x();
62
 
    F += mapCoords[i].y() * pixelCoords[i].y();
63
 
    G += std::pow(pixelCoords[i].x(), 2);
64
 
    H += std::pow(pixelCoords[i].y(), 2);
65
 
    I += mapCoords[i].x() * pixelCoords[i].y();
66
 
    J += pixelCoords[i].x() * mapCoords[i].y();
67
 
  }
68
 
  
69
 
  /* The least squares fit for the parameters { a, b, x0, y0 } is the solution
70
 
     to the matrix equation Mx = b, where M and b is given below. I *think*
71
 
     that this is correct but I derived it myself late at night. Look at 
72
 
     helmert.jpg if you suspect bugs. */
73
 
  
74
 
  double MData[] = { A,   -B,    n,    0,
75
 
                     B,    A,    0,    n,
76
 
                     G+H,  0,    A,    B,                    
77
 
                     0,    G+H, -B,    A };
78
 
 
79
 
  double bData[] = { C,    D,    E+F,  J-I };
80
 
  
81
 
  // we want to solve the equation M*x = b, where x = [a b x0 y0]
82
 
  gsl_matrix_view M = gsl_matrix_view_array(MData, 4, 4);
83
 
  gsl_vector_view b = gsl_vector_view_array(bData, 4);
84
 
  gsl_vector* x = gsl_vector_alloc(4);
85
 
  gsl_permutation* p = gsl_permutation_alloc(4);
86
 
  int s;
87
 
  gsl_linalg_LU_decomp(&M.matrix, p, &s);
88
 
  gsl_linalg_LU_solve(&M.matrix, p, &b.vector, x);
89
 
  gsl_permutation_free(p);
90
 
  
91
 
  origin.setX(gsl_vector_get(x, 2));
92
 
  origin.setY(gsl_vector_get(x, 3));
93
 
  pixelSize = std::sqrt(std::pow(gsl_vector_get(x, 0), 2) +
94
 
                        std::pow(gsl_vector_get(x, 1), 2));
95
 
  rotation = std::atan2(gsl_vector_get(x, 1), gsl_vector_get(x, 0));    
96
 
}
97
 
 
98
 
 
99
 
void QgsLeastSquares::affine(std::vector<QgsPoint> mapCoords,
100
 
                             std::vector<QgsPoint> pixelCoords) {
101
 
  int n = mapCoords.size();
102
 
  if (n < 4) {
103
 
    throw std::domain_error("Fit to an affine transform requires at "
104
 
                            "least 4 points.");
105
 
  }
106
 
  
107
 
  double A = 0, B = 0, C = 0, D = 0, E = 0, F = 0, 
108
 
    G = 0, H = 0, I = 0, J = 0, K = 0;
109
 
  for (int i = 0; i < n; ++i) {
110
 
    A += pixelCoords[i].x();
111
 
    B += pixelCoords[i].y();
112
 
    C += mapCoords[i].x();
113
 
    D += mapCoords[i].y();
114
 
    E += std::pow(pixelCoords[i].x(), 2);
115
 
    F += std::pow(pixelCoords[i].y(), 2);
116
 
    G += pixelCoords[i].x() * pixelCoords[i].y();
117
 
    H += pixelCoords[i].x() * mapCoords[i].x();
118
 
    I += pixelCoords[i].y() * mapCoords[i].y();
119
 
    J += pixelCoords[i].x() * mapCoords[i].y();
120
 
    K += mapCoords[i].x() * pixelCoords[i].y();
121
 
  }
122
 
 
123
 
  /* The least squares fit for the parameters { a, b, c, d, x0, y0 } is the 
124
 
     solution to the matrix equation Mx = b, where M and b is given below. 
125
 
     I *think* that this is correct but I derived it myself late at night. 
126
 
     Look at affine.jpg if you suspect bugs. */
127
 
  
128
 
  double MData[] = { A,    B,    0,    0,    n,    0,
129
 
                     0,    0,    A,    B,    0,    n,
130
 
                     E,    G,    0,    0,    A,    0,
131
 
                     G,    F,    0,    0,    B,    0,
132
 
                     0,    0,    E,    G,    0,    A,
133
 
                     0,    0,    G,    F,    0,    B };
134
 
 
135
 
  double bData[] = { C,    D,    H,    K,    J,    I };
136
 
  
137
 
  // we want to solve the equation M*x = b, where x = [a b c d x0 y0]
138
 
  gsl_matrix_view M = gsl_matrix_view_array(MData, 6, 6);
139
 
  gsl_vector_view b = gsl_vector_view_array(bData, 6);
140
 
  gsl_vector* x = gsl_vector_alloc(6);
141
 
  gsl_permutation* p = gsl_permutation_alloc(6);
142
 
  int s;
143
 
  gsl_linalg_LU_decomp(&M.matrix, p, &s);
144
 
  gsl_linalg_LU_solve(&M.matrix, p, &b.vector, x);
145
 
  gsl_permutation_free(p);
146
 
 
147
 
}
148