USGS

Isis 3.0 Object Programmers' Reference

Home

OverlapNormalization.cpp
Go to the documentation of this file.
1 
23 #include "OverlapNormalization.h"
24 
25 #include <iomanip>
26 
27 #include "BasisFunction.h"
28 #include "IException.h"
29 #include "LeastSquares.h"
30 #include "Statistics.h"
31 
32 using namespace std;
33 
34 namespace Isis {
35 
45  OverlapNormalization::OverlapNormalization(std::vector<Statistics *> statsList) {
46  m_gainFunction = NULL;
47  m_gainLsq = NULL;
48  m_offsetFunction = NULL;
49  m_offsetLsq = NULL;
50 
51  m_statsList = statsList;
52  m_gainFunction = new BasisFunction("BasisFunction",
53  statsList.size(), statsList.size());
54  m_gainLsq = new LeastSquares(*m_gainFunction);
55  m_offsetFunction = new BasisFunction("BasisFunction",
56  statsList.size(), statsList.size());
57  m_offsetLsq = new LeastSquares(*m_offsetFunction);
58 
59  m_gains.resize(statsList.size());
60  m_offsets.resize(statsList.size());
61  for (unsigned int i = 0; i < statsList.size(); i++) {
62  m_gains[i] = 1.0;
63  m_offsets[i] = 0.0;
64  }
65  m_solved = false;
66  }
67 
72  OverlapNormalization::~OverlapNormalization() {
73  if (m_gainFunction != NULL) delete m_gainFunction;
74  if (m_offsetFunction != NULL) delete m_offsetFunction;
75  if (m_gainLsq != NULL) delete m_gainLsq;
76  if (m_offsetLsq != NULL) delete m_offsetLsq;
77  for (unsigned int i = 0; i < m_statsList.size(); i++) delete m_statsList[i];
78  };
79 
106  OverlapNormalization::AddStatus OverlapNormalization::AddOverlap(
107  const Statistics &area1, const unsigned index1,
108  const Statistics &area2, const unsigned index2,
109  double weight) {
110  if (index1 >= m_statsList.size()) {
111  string msg = "The index 1 is outside the bounds of the list.";
112  throw IException(IException::Programmer, msg, _FILEINFO_);
113  }
114  if (index2 >= m_statsList.size()) {
115  string msg = "The index 2 is outside the bounds of the list.";
116  throw IException(IException::Programmer, msg, _FILEINFO_);
117  }
118 
119  // If there is no overlapping area, then the overlap is invalid
120  if (area1.ValidPixels() == 0 || area2.ValidPixels() == 0) {
121  return NoOverlap;
122  }
123 
124  // The weight must be a positive real number
125  if (weight <= 0.0) {
126  string msg = "All weights must be positive real numbers.";
127  throw IException(IException::Programmer, msg, _FILEINFO_);
128  }
129 
131  o.area1 = area1;
132  o.area2 = area2;
133  o.index1 = index1;
134  o.index2 = index2;
135 
136  double avg1 = area1.Average();
137  double avg2 = area2.Average();
138 
139  // Averages must not be 0 to avoid messing up least squares
140  if (avg1 == 0 || avg2 == 0) return NoContrast;
141 
142  m_overlapList.push_back(o);
143  m_deltas.push_back(avg2 - avg1);
144  m_weights.push_back(weight);
145  m_solved = false;
146  return Success;
147  }
148 
163  void OverlapNormalization::Solve(SolutionType type,
164  LeastSquares::SolveMethod method) {
165  // Make sure that there is at least one overlap
166  if (m_overlapList.size() == 0) {
167  string msg = "None of the input images overlap";
168  throw IException(IException::User, msg, _FILEINFO_);
169  }
170 
171  // Make sure the number of valid overlaps + hold images is greater than the
172  // number of input images (otherwise the least squares equation will be
173  // unsolvable due to having more unknowns than knowns)
174  if (m_overlapList.size() + m_idHoldList.size() < m_statsList.size()) {
175  string msg = "Unable to normalize overlaps. The number of overlaps and "
176  "holds must be greater than the number of input images";
177  throw IException(IException::User, msg, _FILEINFO_);
178  }
179 
180  if ( method == LeastSquares::SPARSE ) {
181  m_offsetLsq = NULL;
182  m_gainLsq = NULL;
183  delete m_offsetLsq;
184  delete m_gainLsq;
185  int sparseMatrixRows = m_overlapList.size() + m_idHoldList.size();
186  int sparseMatrixCols = m_offsetFunction->Coefficients();
187  m_offsetLsq = new LeastSquares(*m_offsetFunction, true, sparseMatrixRows, sparseMatrixCols, true);
188  sparseMatrixCols = m_gainFunction->Coefficients();
189  m_gainLsq = new LeastSquares(*m_gainFunction, true, sparseMatrixRows, sparseMatrixCols, true);
190  const std::vector<double> alphaWeight(sparseMatrixCols, 1/1000.0);
191  m_offsetLsq->SetParameterWeights( alphaWeight );
192  m_gainLsq->SetParameterWeights( alphaWeight );
193  }
194 
195  // Calculate offsets
196  if (type != Gains && type != GainsWithoutNormalization) {
197  // Add knowns to least squares for each overlap
198  for (int overlap = 0; overlap < (int)m_overlapList.size(); overlap++) {
199  Overlap curOverlap = m_overlapList[overlap];
200  int id1 = curOverlap.index1;
201  int id2 = curOverlap.index2;
202 
203  vector<double> input;
204  input.resize(m_statsList.size());
205  for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
206  input[id1] = 1.0;
207  input[id2] = -1.0;
208 
209  m_offsetLsq->AddKnown(input, m_deltas[overlap], m_weights[overlap]);
210  }
211 
212  // Add a known to the least squares for each hold image
213  for (int h = 0; h < (int)m_idHoldList.size(); h++) {
214  int hold = m_idHoldList[h];
215 
216  vector<double> input;
217  input.resize(m_statsList.size());
218  for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
219  input[hold] = 1.0;
220  m_offsetLsq->AddKnown(input, 0.0, 1e30);
221  }
222 
223  // Solve the least squares and get the offset coefficients to apply to the
224  // images
225  m_offsets.resize(m_statsList.size());
226  m_offsetLsq->Solve(method);
227  for (int i = 0; i < m_offsetFunction->Coefficients(); i++) {
228  m_offsets[i] = m_offsetFunction->Coefficient(i);
229  }
230  }
231 
232  // Calculate Gains
233  if (type != Offsets) {
234  // Add knowns to least squares for each overlap
235  for (int overlap = 0; overlap < (int)m_overlapList.size(); overlap++) {
236  Overlap curOverlap = m_overlapList[overlap];
237  int id1 = curOverlap.index1;
238  int id2 = curOverlap.index2;
239 
240  vector<double> input;
241  input.resize(m_statsList.size());
242  for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
243  input[id1] = 1.0;
244  input[id2] = -1.0;
245 
246  double tanp;
247 
248  if (type != GainsWithoutNormalization) {
249  if (curOverlap.area1.StandardDeviation() == 0.0) {
250  tanp = 0.0; // Set gain to 1.0
251  }
252  else {
253  tanp = curOverlap.area2.StandardDeviation()
254  / curOverlap.area1.StandardDeviation();
255  }
256  }
257  else {
258  if (curOverlap.area1.Average() == 0.0) {
259  tanp = 0.0;
260  }
261  else {
262  tanp = curOverlap.area2.Average() / curOverlap.area1.Average();
263  }
264  }
265 
266  if (tanp > 0.0) {
267  m_gainLsq->AddKnown(input, log(tanp), m_weights[overlap]);
268  }
269  else {
270  m_gainLsq->AddKnown(input, 0.0, 1e30); // Set gain to 1.0
271  }
272  }
273 
274  // Add a known to the least squares for each hold image
275  for (int h = 0; h < (int)m_idHoldList.size(); h++) {
276  int hold = m_idHoldList[h];
277 
278  vector<double> input;
279  input.resize(m_statsList.size());
280  for (int i = 0; i < (int)input.size(); i++) input[i] = 0.0;
281  input[hold] = 1.0;
282  m_gainLsq->AddKnown(input, 0.0, 1e30);
283  }
284 
285  // Solve the least squares and get the gain coefficients to apply to the
286  // images
287  m_gains.resize(m_statsList.size());
288  m_gainLsq->Solve(method);
289  for (int i = 0; i < m_gainFunction->Coefficients(); i++) {
290  m_gains[i] = exp(m_gainFunction->Coefficient(i));
291  }
292  }
293 
294  m_solved = true;
295  }
296 
309  double OverlapNormalization::Average(const unsigned index) const {
310  if (index >= m_statsList.size()) {
311  string msg = "The index was out of bounds for the list of statistics.";
312  throw IException(IException::Programmer, msg, _FILEINFO_);
313  }
314 
315  return m_statsList[index]->Average();
316  }
317 
330  double OverlapNormalization::Gain(const unsigned index) const {
331  if (index >= m_statsList.size()) {
332  string msg = "The index was out of bounds for the list of statistics.";
333  throw IException(IException::Programmer, msg, _FILEINFO_);
334  }
335 
336  return m_gains[index];
337  }
338 
351  double OverlapNormalization::Offset(const unsigned index) const {
352  if (index >= m_statsList.size()) {
353  string msg = "The index was out of bounds for the list of statistics.";
354  throw IException(IException::Programmer, msg, _FILEINFO_);
355  }
356 
357  return m_offsets[index];
358  }
359 
374  double OverlapNormalization::Evaluate(double dn, unsigned index) const {
375  if (!m_solved) {
376  string msg = "The least squares equation has not been successfully ";
377  msg += "solved yet.";
378  throw IException(IException::Programmer, msg, _FILEINFO_);
379  }
380 
381  if (Isis::IsSpecial(dn)) return dn;
382  return (dn - Average(index)) * Gain(index) + Average(index) + Offset(index);
383  }
384 }