Isis 3.0 Application Source Code Reference |
Home |
00001 #if !defined(HiCalUtil_h) 00002 #define HiCalUtil_h 00003 /** 00004 * @file 00005 * $Revision: 3878 $ 00006 * $Date: 2012-01-18 10:23:13 -0700 (Wed, 18 Jan 2012) $ 00007 * 00008 * Unless noted otherwise, the portions of Isis written by the USGS are 00009 * public domain. See individual third-party library and package descriptions 00010 * for intellectual property information, user agreements, and related 00011 * information. 00012 * 00013 * Although Isis has been used by the USGS, no warranty, expressed or 00014 * implied, is made by the USGS as to the accuracy and functioning of such 00015 * software and related material nor shall the fact of distribution 00016 * constitute any such warranty, and no responsibility is assumed by the 00017 * USGS in connection therewith. 00018 * 00019 * For additional information, launch 00020 * $ISISROOT/doc//documents/Disclaimers/Disclaimers.html 00021 * in a browser or see the Privacy & Disclaimers page on the Isis website, 00022 * http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on 00023 * http://www.usgs.gov/privacy.html. 00024 */ 00025 00026 #include <cmath> 00027 #include <string> 00028 #include <vector> 00029 #include <iostream> 00030 #include <sstream> 00031 00032 #include "HiCalTypes.h" 00033 #include "DbProfile.h" 00034 #include "iString.h" 00035 #include "Statistics.h" 00036 #include "PvlKeyword.h" 00037 #include "NumericalApproximation.h" 00038 #include "Filename.h" 00039 #include "Pvl.h" 00040 #include "iException.h" 00041 00042 namespace Isis { 00043 00044 template <typename T> inline T MIN(const T &A, const T &B) { 00045 if ( A < B ) { return (A); } 00046 else { return (B); } 00047 } 00048 00049 template <typename T> inline T MAX(const T &A, const T &B) { 00050 if ( A > B ) { return (A); } 00051 else { return (B); } 00052 } 00053 00054 /** 00055 * @brief Counts number of valid pixels in vector 00056 * @param v Vector to inspect 00057 * @return int Number valid pixels in vector 00058 */ 00059 inline int ValidCount(const HiVector &v) { 00060 int n = 0; 00061 for (int i = 0 ; i < v.dim() ; i++) { 00062 if (!IsSpecial(v[i])) n++; 00063 } 00064 return (n); 00065 } 00066 00067 /** 00068 * @brief Counts number of invalid pixels in vector 00069 * @param v Vector to inspect 00070 * @return int Number invalid (special) pixels in vector 00071 */ 00072 inline int InValidCount(const HiVector &v) { 00073 int n = 0; 00074 for (int i = 0 ; i < v.dim() ; i++) { 00075 if (IsSpecial(v[i])) n++; 00076 } 00077 return (n); 00078 } 00079 00080 /** 00081 * @brief Convert HiRISE Cpmm number to Ccd number 00082 * 00083 * @param cpmm Cpmm number 00084 */ 00085 inline int CpmmToCcd(int cpmm) { 00086 const int cpmm2ccd[] = {0,1,2,3,12,4,10,11,5,13,6,7,8,9}; 00087 if ( (cpmm < 0) || (cpmm >= (int)(sizeof(cpmm2ccd)/sizeof(int))) ) { 00088 std::string mess = "CpmmToCdd: Bad CpmmNumber (" + iString(cpmm) + ")"; 00089 throw iException::Message(iException::User, mess.c_str(), _FILEINFO_); 00090 } 00091 return (cpmm2ccd[cpmm]); 00092 } 00093 00094 00095 /** 00096 * @brief Convert HiRISE Ccd number to string filter name 00097 * 00098 * @param ccd Ccd number of device 00099 */ 00100 inline std::string CcdToFilter(int ccd) { 00101 if ( (ccd < 0) || (ccd > 13) ) { 00102 std::string mess = "CcdToFilter: Bad Ccd Number (" + iString(ccd) + ")"; 00103 throw iException::Message(iException::User, mess.c_str(), _FILEINFO_); 00104 } 00105 00106 std::string filter; 00107 if ( ccd <= 9 ) { filter = "RED"; } 00108 else if (ccd <= 11) { filter = "IR"; } 00109 else { filter = "BG"; } 00110 return (filter); 00111 } 00112 00113 /** 00114 * @brief Crop specified lines from a buffer 00115 * 00116 * This function extracts lines from a buffer and returns a new buffer with the 00117 * extracted lines. 00118 * 00119 * @param m Buffer to extract lines from 00120 * @param sline Starting line number (first line is 0) 00121 * @param eline Last line to extract 00122 * 00123 * @return HiMatrix Buffer containing cropped lines 00124 */ 00125 inline HiMatrix cropLines(const HiMatrix &m, int sline, int eline) { 00126 int nlines(eline-sline+1); 00127 HiMatrix mcrop(nlines, m.dim2()); 00128 for (int l = 0 ; l < nlines ; l++) { 00129 for (int s = 0 ; s < m.dim2() ; s++) { 00130 mcrop[l][s] = m[l+sline][s]; 00131 } 00132 } 00133 return (mcrop); 00134 } 00135 00136 /** 00137 * @brief Crop specified samples from a buffer 00138 * 00139 * This function extracts samples from a buffer and returns a new buffer with 00140 * the extracted samples. 00141 * 00142 * @param m Buffer to extract samples from 00143 * @param ssamp Startling sample (first sample 0) 00144 * @param esamp Ending sample to extract 00145 * 00146 * @return HiMatrix Buffer with cropped samples 00147 */ 00148 inline HiMatrix cropSamples(const HiMatrix &m, int ssamp, int esamp) { 00149 int nsamps(esamp-ssamp+1); 00150 HiMatrix mcrop(m.dim1(), nsamps); 00151 for (int l = 0 ; l < m.dim1() ; l++) { 00152 for (int s = 0 ; s < nsamps ; s++) { 00153 mcrop[l][s] = m[l][s+ssamp]; 00154 } 00155 } 00156 return (mcrop); 00157 } 00158 00159 /** 00160 * @brief Reduces by averaging specified lines from a buffer 00161 * 00162 * This function produces a vector from a 2-D buffer of averaged lines at each 00163 * vertical sample location 00164 * 00165 * @param m Buffer to reduce 00166 * @param sline Starting line number (first line is 0) 00167 * @param eline Last line to average (-1 means use all lines) 00168 * 00169 * @return HiVector Buffer containing averaged lines 00170 */ 00171 inline HiVector averageLines(const HiMatrix &m, int sline = 0, int eline = -1) { 00172 eline = (eline == -1) ? m.dim1() - 1 : eline; 00173 HiVector v(m.dim2()); 00174 for (int s = 0 ; s < m.dim2() ; s++) { 00175 Statistics stats; 00176 for (int l = sline ; l <= eline ; l++) { 00177 stats.AddData(&m[l][s], 1); 00178 } 00179 v[s] = stats.Average(); 00180 } 00181 return (v); 00182 } 00183 00184 /** 00185 * @brief Reduces by averaging specified samples from a buffer 00186 * 00187 * This function produces a vector from a 2-D buffer of averaged samples at each 00188 * horizontal line location 00189 * 00190 * @param m Buffer to reduce 00191 * @param ssamp Starting sample number (first sample is 0) 00192 * @param esamp Last sample to average (-1 means use all samples) 00193 * 00194 * @return HiVector Buffer containing averaged samples 00195 */ 00196 inline HiVector averageSamples(const HiMatrix &m, int ssamp = 0, 00197 int esamp = -1) { 00198 esamp = (esamp == -1) ? m.dim2() - 1 : esamp; 00199 HiVector v(m.dim1()); 00200 for (int l = 0 ; l < m.dim1() ; l++) { 00201 Statistics stats; 00202 for (int s = ssamp ; s <= esamp ; s++) { 00203 stats.AddData(&m[l][s], 1); 00204 } 00205 v[l] = stats.Average(); 00206 } 00207 return (v); 00208 } 00209 00210 /** 00211 * @brief Find a keyword in a profile using default for non-existant keywords 00212 * 00213 * This template function will extract keyword values from a profile and provide 00214 * the values at the indicated index. If the keyword does not exist or the 00215 * indicated value at index, the provided default will be used instead. 00216 * 00217 */ 00218 template <typename T> 00219 T ConfKey(const DbProfile &conf, const std::string &keyname, const T &defval, 00220 int index = 0) { 00221 if (!conf.exists(keyname)) { return (defval); } 00222 if (conf.count(keyname) < index) { return (defval); } 00223 iString iValue(conf.value(keyname, index)); 00224 T value = iValue; // This makes it work with a string? 00225 return (value); 00226 } 00227 00228 /** 00229 * @brief Helper function to convert values to Integers 00230 * 00231 * @param T Type of value to convert 00232 * @param value Value to convert 00233 * 00234 * @return int Converted value 00235 */ 00236 template <typename T> int ToInteger(const T &value) { 00237 return (iString(value).Trim(" \r\t\n").ToInteger()); 00238 } 00239 00240 /** 00241 * @brief Helper function to convert values to doubles 00242 * 00243 * @param T Type of value to convert 00244 * @param value Value to convert 00245 * 00246 * @return double Converted value 00247 */ 00248 template <typename T> double ToDouble(const T &value) { 00249 return (iString(value).Trim(" \r\t\n").ToDouble()); 00250 } 00251 00252 /** 00253 * @brief Helper function to convert values to strings 00254 * 00255 * @param T Type of value to convert 00256 * @param value Value to convert 00257 * 00258 * @return string Converted value 00259 */ 00260 template <typename T> std::string ToString(const T &value) { 00261 return (iString(value).Trim(" \r\t\n")); 00262 } 00263 00264 /** 00265 * @brief Shortened string equality test 00266 * 00267 * @param v1 First value 00268 * @param v2 Second value 00269 * 00270 * @return bool True if they are equal w/o regard to case 00271 */ 00272 inline bool IsEqual(const std::string &v1, const std::string &v2 = "TRUE") { 00273 return (iString::UpCase(v1) == iString::UpCase(v2)); 00274 } 00275 00276 /** 00277 * @brief Determines if the keyword value is the expected value 00278 * 00279 * This function checks the existance of a keyword in a profile, extracts the 00280 * first value and tests it for equivelance to the expected value. Note this 00281 * test is case insensitive. 00282 * 00283 * @param prof Profile to find the expected keyword in 00284 * @param key Name of keyword in profile to test 00285 * @param value Value to test for in keyword 00286 * 00287 * @return bool Returns true only if the keyword exists in the profile and it is 00288 * the value expected. 00289 */ 00290 inline bool IsTrueValue(const DbProfile &prof, const std::string &key, 00291 const std::string &value = "TRUE") { 00292 if ( prof.exists(key) ) { 00293 return (IsEqual(prof(key), value)); 00294 } 00295 return (false); 00296 } 00297 00298 /** 00299 * @brief Checks profile flag to skip the current Module 00300 * 00301 * This function looks for the keyword \b Debug::SkipModule and checks its 00302 * value. True is returned if the value is TRUE (case insensentive). 00303 * 00304 * @param prof Module profile from config file 00305 * 00306 * @return bool True if the value of the Debug::SkipModule keyword is TRUE, 00307 * otherwise it returns false for all other values. 00308 */ 00309 inline bool SkipModule(const DbProfile &prof) { 00310 return (IsTrueValue(prof, "Debug::SkipModule", "TRUE")); 00311 } 00312 00313 00314 inline HiMatrix appendLines(const HiMatrix &top, const HiMatrix &bottom) { 00315 // Ensure same number samples 00316 if (top.dim2() != bottom.dim2()) { 00317 std::ostringstream mess; 00318 mess << "Top buffer samples (" << top.dim2() 00319 << ") do not match bottom buffer samples (" << bottom.dim2() 00320 << ")"; 00321 throw iException::Message(iException::User, mess.str(), _FILEINFO_); 00322 } 00323 00324 int nlines(top.dim1()+bottom.dim1()); 00325 HiMatrix mat(nlines, top.dim2()); 00326 for (int lt = 0 ; lt < top.dim1() ; lt++) { 00327 for (int s = 0 ; s < top.dim2() ; s++) { 00328 mat[lt][s] = top[lt][s]; 00329 } 00330 } 00331 00332 int topl = top.dim1(); 00333 for (int lb = 0 ; lb < bottom.dim1() ; lb++) { 00334 for (int s = 0 ; s < top.dim2() ; s++) { 00335 mat[topl+lb][s] = bottom[lb][s]; 00336 } 00337 } 00338 return (mat); 00339 } 00340 00341 inline HiMatrix appendSamples(const HiMatrix &left, const HiMatrix &right) { 00342 // Ensure same number samples 00343 if (right.dim1() != right.dim1()) { 00344 std::ostringstream mess; 00345 mess << "Left buffer lines (" << left.dim1() 00346 << ") do not match right buffer lines (" << right.dim1() 00347 << ")"; 00348 throw iException::Message(iException::User, mess.str(), _FILEINFO_); 00349 } 00350 00351 int nsamps(left.dim2()+right.dim2()); 00352 HiMatrix mat(left.dim1(), nsamps); 00353 for (int ll = 0 ; ll < left.dim1() ; ll++) { 00354 for (int s = 0 ; s < left.dim2() ; s++) { 00355 mat[ll][s] = left[ll][s]; 00356 } 00357 } 00358 00359 int lefts = left.dim2(); 00360 for (int lr = 0 ; lr < right.dim1() ; lr++) { 00361 for (int s = 0 ; s < right.dim2() ; s++) { 00362 mat[lr][lefts+s] = right[lr][s]; 00363 } 00364 } 00365 return (mat); 00366 } 00367 00368 /** 00369 * @brief Compute HiRISE line times 00370 * 00371 * This class will compute the time (in seconds) for a HiRISE observation line 00372 * based upon the binning mode and line number. It is assumed that the first 00373 * line will be time 0, but that is up to the caller to keep that straight. 00374 * 00375 * @author ????-??-?? Unknown 00376 * 00377 * @internal 00378 */ 00379 class HiLineTimeEqn { 00380 public: 00381 HiLineTimeEqn() : _bin(1), _ltime(1.0) { } 00382 HiLineTimeEqn(int bin, double ltime) : _bin(bin), _ltime(ltime) { } 00383 virtual ~HiLineTimeEqn() { } 00384 00385 void setLineTime(double ltime) { _ltime = ltime; } 00386 void setBin(int bin) { _bin = bin; } 00387 double Time(const double line) const { 00388 return (line * (_bin * _ltime * 1.0E-6)); 00389 } 00390 double operator()(const double line) const { return (Time(line)); } 00391 00392 private: 00393 double _bin; 00394 double _ltime; 00395 }; 00396 00397 /** 00398 * @brief Implements (classic) HiRISE temperature equation 00399 * 00400 * This function computes the dark current temperature and returns the results 00401 * in electrons/sec/pixel. 00402 * 00403 * @param temperature Temperature (typically of the FPGA) 00404 * @param napcm2 Dark current for silicone diodes (nano-ampres/cm^2) 00405 * @param px Pixel size in microns 00406 * 00407 * @return double Dark current temperature (electrons/sec/pixel) 00408 */ 00409 inline double HiTempEqn(const double temperature, const double napcm2 = 2.0, 00410 const double px = 12.0) { 00411 double temp = temperature + 273.0; 00412 double eg = 1.1557 - (7.021E-4 * temp * temp) / (1108.0 + temp); 00413 const double K = 1.38E-23; 00414 const double Q = 1.6E-19; 00415 return (napcm2*(px*px)*2.55E7*std::pow(temp,1.5)*std::exp(-eg*Q/2.0/K/temp)); 00416 } 00417 00418 /** 00419 * @brief Rebins a vector to a different size 00420 * 00421 * This function can rebin to both larger and smaller sizes. It fits the data 00422 * to a cubic spline and then computes the value at the rebin pixel index. One 00423 * advantage to this approach is that on input, special pixels are ignored - on 00424 * output there will never be special pixels unless there are not enough points 00425 * to conpute the cubic spline on which case this function throws an exception. 00426 * 00427 * @param v Input vector to rebin 00428 * @param n Size of the new output vector 00429 * 00430 * @return HiVector 00431 * @history 2008-11-05 Jeannie Walldren Replaced references to 00432 * DataInterp class with NumericalApproximation. 00433 */ 00434 inline HiVector rebin(const HiVector &v, int n) { 00435 if ( n == v.dim() ) { return (v); } 00436 NumericalApproximation nterp(NumericalApproximation::CubicNatural); 00437 double mag((double) v.dim()/ (double) n); 00438 00439 for ( int i = 0 ; i < v.dim() ; i++ ) { 00440 double x((double) i); 00441 if ( !IsSpecial(v[i]) ) { 00442 nterp.AddData(x,v[i]); 00443 } 00444 } 00445 00446 // Compute the spline and fill the output vector 00447 HiVector vout(n); 00448 for ( int j = 0 ; j < n ; j++ ) { 00449 double x = (double) j * mag; 00450 vout[j] = nterp.Evaluate(x, NumericalApproximation::NearestEndpoint); 00451 } 00452 return (vout); 00453 } 00454 00455 /** 00456 * @brief Deletes HiRISE specific BLOBS from cube file 00457 * 00458 * Ths function removes only the HiRISE specific 00459 * 00460 * @param label Input label associated with file from which to remove the HiRISE 00461 * blobs 00462 */ 00463 inline void RemoveHiBlobs(Pvl &label) { 00464 for ( int blob = 0 ; blob < label.Objects() ; blob++ ) { 00465 if ( label.Object(blob).IsNamed("Table") ) { 00466 std::string name = label.Object(blob)["Name"][0]; 00467 if ( iString::Equal(name,"HiRISE Calibration Ancillary") || 00468 iString::Equal(name,"HiRISE Calibration Image") || 00469 iString::Equal(name,"HiRISE Ancillary") ) { 00470 label.DeleteObject(blob); 00471 blob--; 00472 } 00473 } 00474 } 00475 return; 00476 } 00477 00478 /** 00479 * @brief Return the statistics for a vector of data 00480 * 00481 * The default statistic returned is the median of the dataset but can be 00482 * changed with a compile time change. The vector passed in will be sorted so 00483 * that the median can be determined. If the vector has an even number of 00484 * elements in it, the average of the two center values will be returned. The 00485 * array is assumed to be clean data, no ISIS special pixels. 00486 * 00487 * @param data Vector containing data compute the statistic. It will be 00488 * altered (sorted) upon return to the caller. 00489 * 00490 * @return double The median (default) of the data 00491 */ 00492 inline double GainLineStat(std::vector<double> &data) { 00493 00494 // Check for special conditions 00495 if (data.size() == 0) return (Null); 00496 if (data.size() == 1) return (data[0]); 00497 00498 #if defined(USE_AVERAGE) 00499 Statistics stats; 00500 stats.AddData(&data[0], data.size()); 00501 return (stat.Average()); 00502 #else 00503 std::sort(data.begin(), data.end()); 00504 int meanIndex = data.size()/2; 00505 if ((data.size() % 2) == 1) return data[meanIndex]; 00506 return ((data[meanIndex-1]+data[meanIndex])/2.0); 00507 #endif 00508 } 00509 00510 }; // namespace Isis 00511 #endif