USGS

Isis 3.0 Object Programmers' Reference

 Home
GruenTypes.h
Go to the documentation of this file.
1 #ifndef GruenTypes_h
2 #define GruenTypes_h
3 
26 #include <cmath>
27 #include <iostream>
28 #include <sstream>
29 #include <iomanip>
30 
31 #include <tnt/tnt_array1d.h>
32 #include <tnt/tnt_array2d.h>
33 #include <tnt/tnt_array2d_utils.h>
34 
35 #include "Affine.h"
36 #include "Chip.h"
37 #include "Constants.h"
38 #include "IException.h"
39 #include "PvlKeyword.h"
40 #include "SpecialPixel.h"
41 
42 
43 namespace Isis {
44 
45  typedef Affine::AMatrix GMatrix;
46  typedef TNT::Array1D<double> GVector;
47 
48  // Constraints
49  enum { NCONSTR = 8 };
50 
69  class Coordinate {
70  public:
71  Coordinate() : m_y(Null), m_x(Null) { }
72  Coordinate(double y, double x) : m_y(y), m_x(x) { }
73  Coordinate(const Chip &chip) : m_y(chip.CubeLine()),
74  m_x(chip.CubeSample()){ }
75  ~Coordinate() { }
76 
78  void setLatLon(const double &latitude, const double &longitude) {
79  m_y = latitude;
80  m_x = longitude;
81  return;
82  }
83 
85  void setLineSamp(const double &line, const double &sample) {
86  m_y = line;
87  m_x = sample;
88  return;
89  }
90 
92  Coordinate &operator+=(const Coordinate &other) {
93  if ( isValid() && other.isValid()) {
94  m_y += other.m_y;
95  m_x += other.m_x;
96  }
97  else {
98  m_y = m_x = Null;
99  }
100  return (*this);
101  }
102 
104  Coordinate &operator-=(const Coordinate &other) {
105  if ( isValid() && other.isValid() ) {
106  m_y -= other.m_y;
107  m_x -= other.m_x;
108  }
109  else {
110  m_y = m_x = Null;
111  }
112  return (*this);
113  }
114 
116  double getDistance(const Coordinate &pntA = Coordinate(0.0, 0.0)) const {
117  double yd = pntA.m_y - m_y;
118  double xd = pntA.m_x - m_x;
119  return (std::sqrt((xd * xd) + (yd * yd)));
120  }
121 
123  inline bool isValid() const {
124  return ( !(IsSpecial(m_x) || IsSpecial(m_y)));
125  }
126 
127  inline double getLatitude() const { return (m_y); }
128  inline double getLongitude() const { return (m_x); }
129  inline double getLine() const { return (m_y); }
130  inline double getSample() const { return (m_x); }
131 
132  double m_y; // Overloaded as Latitude or line
133  double m_x; // Overloaded as Longitude or sample
134  };
135 
136 
146  inline Coordinate operator+(const Coordinate &A, const Coordinate &B) {
147  if ( A.isValid() && B.isValid() ) {
148  return (Coordinate(A.m_y+B.m_y, A.m_x+B.m_x));
149  }
150  else {
151  return (Coordinate());
152  }
153  }
154 
164  inline Coordinate operator-(const Coordinate &A, const Coordinate &B) {
165  if ( A.isValid() && B.isValid() ) {
166  return (Coordinate(A.m_y-B.m_y, A.m_x-B.m_x));
167  }
168  else {
169  return (Coordinate());
170  }
171  }
172 
187  class PointPair {
188  public:
189  PointPair() : m_left(), m_right() { }
190  PointPair(const double &line, const double &sample) : m_left(line, sample),
191  m_right() { }
192  PointPair(const Coordinate &left,
193  const Coordinate &right = Coordinate()) : m_left(left),
194  m_right(right) { }
195 
197  inline bool isValid() const {
198  return (m_left.isValid() && m_right.isValid());
199  }
200  inline const Coordinate &getLeft() const { return (m_left); }
201  inline const Coordinate &getRight() const { return (m_right); }
202 
203  inline double getLine() const { return (getLeftLine()); }
204  inline double getSample() const { return (getLeftSample()); }
205  inline double getLeftLine() const { return (m_left.getLine()); }
206  inline double getLeftSample() const { return (m_left.getSample()); }
207  inline double getRightLine() const { return (m_right.getLine()); }
208  inline double getRightSample() const { return (m_right.getSample()); }
209 
210  Coordinate m_left; // Left image line/sample
211  Coordinate m_right; // Right image line/sample
212  };
213 
214 
222  class Radiometric {
223  public:
224  Radiometric() : m_shift(0.0), m_gain(0.0) { }
225  Radiometric(const double &shift, const double &gain) :
226  m_shift(shift), m_gain(gain) { }
227 
228  inline double Shift() const { return (m_shift); }
229  inline double Gain() const { return (m_gain); }
230 
231 
233  Radiometric &operator+=(const Radiometric &B) {
234  m_shift += B.m_shift;
235  m_gain += B.m_gain;
236  return (*this);
237  }
238 
239  double m_shift; // Radiometric shift
240  double m_gain; // Radiometric gain
241  };
242 
244  inline Radiometric operator+(const Radiometric &A, const Radiometric &B) {
245  return (Radiometric(A.m_shift+B.m_shift, A.m_gain+B.m_gain));
246  }
247 
248 
258  class AffineRadio {
259  public:
260  AffineRadio() : m_affine(Affine::getIdentity()), m_radio() {}
261  AffineRadio(const GMatrix &A) : m_affine(A), m_radio() { }
262  AffineRadio(const GMatrix &M, const double &shift,
263  const double &gain) : m_affine(M), m_radio(shift, gain) { }
264  AffineRadio(const GVector &alpha) : m_affine(), m_radio() {
265  clone(alpha); // Clone from an alpha matrix
266  }
267  AffineRadio(const Radiometric &radio) : m_affine(Affine::getIdentity()),
268  m_radio(radio) { }
269  ~AffineRadio() { }
270 
272  AffineRadio &operator+=(const AffineRadio &other) {
273  m_affine = m_affine + (other.m_affine - Affine::getIdentity());
274  m_radio += other.m_radio;
275  return (*this);
276  }
277 
279  void Translate(const Coordinate &offset) {
280  GMatrix trans = Affine::getIdentity();
281  trans[0][2] = offset.getSample();
282  trans[1][2] = offset.getLine();
283  m_affine = TNT::matmult(trans, m_affine);
284  return;
285  }
286 
288  Coordinate getPoint(const Coordinate &location) const {
289  double x = m_affine[0][0] * location.getSample() +
290  m_affine[0][1] * location.getLine() + m_affine[0][2];
291  double y = m_affine[1][0] * location.getSample() +
292  m_affine[1][1] * location.getLine() + m_affine[1][2];
293  return (Coordinate(y, x));
294  }
295 
296  GMatrix m_affine; // Affine transform
297  Radiometric m_radio; // Radiometric gain and shift
298 
299  private:
301  void clone(const GVector &alpha) {
302  if ( alpha.dim1() != 8 ) {
303  QString mess = "Alpha array for AffineRadio must have 8 elements "
304  " but has " + toString(alpha.dim1());
305  throw IException(IException::Programmer, mess, _FILEINFO_);
306  }
307  m_affine = Affine::getIdentity();
308  m_affine[0][0] += alpha[1];
309  m_affine[0][1] += alpha[2];
310  m_affine[0][2] += alpha[0];
311 
312  m_affine[1][0] += alpha[4];
313  m_affine[1][1] += alpha[5];
314  m_affine[1][2] += alpha[3];
315 
316  m_radio.m_shift = alpha[6];
317  m_radio.m_gain = alpha[7];
318  }
319  };
320 
321 
335  struct AffineTolerance {
336  public:
337  AffineTolerance() : m_transTol(0.1), m_scaleTol(0.5), m_shearTol(0.5) { }
338  AffineTolerance(const double &transTol, const double &scaleTol,
339  const double &shearTol) : m_transTol(transTol),
340  m_scaleTol(scaleTol),
341  m_shearTol(shearTol) { }
342  ~AffineTolerance() { }
343 
344  double m_transTol; // Affine translation tolerance
345  double m_scaleTol; // Affine scale tolerance
346  double m_shearTol; // Affine shear tolerance
347  };
348 
349 
364  class Threshold {
365  public:
366  Threshold() : m_thresh(6, 0.0) { }
367  Threshold(const Chip &chip, const AffineTolerance &tolerance) : m_thresh(6) {
368  m_thresh[0] = tolerance.m_scaleTol / (((double)(chip.Samples() - 1)) / 2.0);
369  m_thresh[1] = tolerance.m_shearTol / (((double)(chip.Lines() - 1)) / 2.0);
370  m_thresh[2] = tolerance.m_transTol;
371 
372  m_thresh[3] = tolerance.m_shearTol / (((double)(chip.Samples() - 1)) / 2.0);
373  m_thresh[4] = tolerance.m_scaleTol / (((double)(chip.Lines() - 1)) / 2.0);
374  m_thresh[5] = tolerance.m_transTol;
375  }
376  ~Threshold() { }
377 
379  bool hasConverged(const AffineRadio &affine) const {
380  GMatrix Malpha = affine.m_affine - Affine::getIdentity();
381  const double *alpha = Malpha[0];
382  for ( int i = 0 ; i < m_thresh.dim1() ; i++ ) {
383  if ( std::fabs(alpha[i]) >= m_thresh[i] ) return (false);
384  }
385  return (true); // If all values are below threshold
386  }
387 
388  GVector m_thresh;
389  };
390 
391 
399  struct Analysis {
400  Analysis() : m_npts(0), m_variance(0.0), m_sevals(),
401  m_kmat(), m_status(-1) {
402  for ( int i = 0 ; i < 2 ; i++ ) {
403  m_sevals[i] = 999.0;
404  m_kmat[i] = 999.0;
405  }
406  }
407  ~Analysis() { }
408 
409  inline bool isValid() const { return (m_status == 0); }
410 
412  inline double getEigen() const {
413  double eigen = std::sqrt(m_sevals[0] * m_sevals[0] +
414  m_sevals[1] * m_sevals[1]);
415  return (eigen);
416  }
417 
419  inline void setZeroState() {
420  for ( int i = 0 ; i < 2 ; i++ ) {
421  m_sevals[i] = 0.0;
422  m_kmat[i] = 0.0;
423  }
424  m_status = 0;
425  return;
426  }
427 
428  BigInt m_npts;
429  double m_variance;
430  double m_sevals[2]; // 2 sorted eigenvalues
431  double m_kmat[2]; // Sample/Line uncertainty
432  int m_status; // Status
433  };
434 
435 
449  class MatchPoint {
450  public:
451  MatchPoint() : m_point(), m_affine(), m_analysis(), m_status(-1) { }
452  MatchPoint(const AffineRadio &radio) : m_point(), m_affine(radio),
453  m_analysis(), m_status(-1) { }
454  MatchPoint(const PointPair &point) : m_point(point), m_affine(),
455  m_analysis(), m_status(-1) { }
456  ~MatchPoint() { }
457 
458  inline int getStatus() const { return (m_status); }
459  const MatchPoint &setStatus(int status) {
460  m_status = status;
461  return (*this);
462  }
463 
464  inline bool isValid() const { return (m_status == 0); }
465  inline double getEigen() const { return (m_analysis.getEigen()); }
466 
468  Coordinate getAffinePoint(const Coordinate &coord = Coordinate(0.0, 0.0))
469  const {
470  return (m_affine.getPoint(coord));
471  }
472 
473  PointPair m_point; // Pattern (left) and search (right) points
474  AffineRadio m_affine; // Resulting Affine transform
475  Analysis m_analysis; // Error analysis of registration
476  int m_nIters; // Number of iterations required to match
477  int m_status; // Status - good is 0
478  };
479 
480 } // namespace Isis
481 #endif