CCfits
2.5
|
00001 // Astrophysics Science Division, 00002 // NASA/ Goddard Space Flight Center 00003 // HEASARC 00004 // http://heasarc.gsfc.nasa.gov 00005 // e-mail: ccfits@legacy.gsfc.nasa.gov 00006 // 00007 // Original author: Ben Dorman 00008 00009 #ifndef IMAGE_H 00010 #define IMAGE_H 1 00011 00012 // functional 00013 #include <functional> 00014 // valarray 00015 #include <valarray> 00016 // vector 00017 #include <vector> 00018 // numeric 00019 #include <numeric> 00020 #include <sstream> 00021 00022 #ifdef _MSC_VER 00023 #include "MSconfig.h" //form std::min 00024 #endif 00025 #include "CCfits.h" 00026 #include "FitsError.h" 00027 #include "FITSUtil.h" 00028 00029 00030 namespace CCfits { 00031 00032 00033 00034 template <typename T> 00035 class Image 00036 { 00037 00038 public: 00039 Image(const Image< T > &right); 00040 Image (const std::valarray<T>& imageArray = std::valarray<T>()); 00041 ~Image(); 00042 Image< T > & operator=(const Image< T > &right); 00043 00044 const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls); 00045 const std::valarray<T>& readImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, T* nullValue, const std::vector<long>& naxes, bool& nulls); 00046 // If write operation causes an expansion of the image's outer-most dimension, newNaxisN will be set to the new value. Else it will be 0. 00047 void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN, T* nullValue = 0); 00048 void writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN); 00049 void writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN); 00050 bool isRead () const; 00051 // This allows higher level classes to notify Image that a user-input 00052 // scaling value has changed. Image can then decide how this 00053 // should affect reading from cache vs. disk. 00054 void scalingHasChanged(); 00055 // Give the user (via higher level classes) a way to explicitly set the m_isRead flag 00056 // to false, thus providing a fail-safe override of reading from the cache. 00057 void resetRead(); 00058 const std::valarray< T >& image () const; 00059 00060 // Additional Public Declarations 00061 00062 protected: 00063 // Additional Protected Declarations 00064 00065 private: 00066 std::valarray<T>& image (); 00067 void prepareForSubset (const std::vector<long>& naxes, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, const std::valarray<T>& inData, std::valarray<T>& subset); 00068 void loop (size_t iDim, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, size_t iPos, const std::vector<size_t>& incr, const std::valarray<T>& inData, size_t& iDat, const std::vector<size_t>& subIncr, std::valarray<T>& subset, size_t iSub); 00069 bool isNullValChanged(T* newNull) const; 00070 void setLastNullInfo(T* newNull); 00071 00072 // Additional Private Declarations 00073 00074 private: //## implementation 00075 // Data Members for Class Attributes 00076 // When m_isRead = true, assume m_fullImageCache contains the full image from the file. 00077 bool m_isRead; 00078 00079 // Information regarding the usage of null values for the 00080 // most recent read operation. 00081 bool m_usingNullVal; 00082 T m_lastNullVal; 00083 00084 // Data Members for Associations 00085 std::valarray< T > m_fullImageCache; 00086 std::valarray<T> m_currentRead; 00087 00088 // Additional Implementation Declarations 00089 00090 }; 00091 00092 // Parameterized Class CCfits::Image 00093 00094 template <typename T> 00095 inline bool Image<T>::isRead () const 00096 { 00097 return m_isRead; 00098 } 00099 00100 template <typename T> 00101 inline const std::valarray< T >& Image<T>::image () const 00102 { 00103 return m_fullImageCache; 00104 } 00105 00106 00107 // Parameterized Class CCfits::Image 00108 00109 template <typename T> 00110 Image<T>::Image(const Image<T> &right) 00111 : m_isRead(right.m_isRead), 00112 m_usingNullVal(right.m_usingNullVal), 00113 m_lastNullVal(right.m_lastNullVal), 00114 m_fullImageCache(right.m_fullImageCache), 00115 m_currentRead(right.m_currentRead) 00116 { 00117 } 00118 00119 template <typename T> 00120 Image<T>::Image (const std::valarray<T>& imageArray) 00121 : m_isRead(false), 00122 m_usingNullVal(false), 00123 m_lastNullVal(0), 00124 m_fullImageCache(imageArray), 00125 m_currentRead() 00126 { 00127 } 00128 00129 00130 template <typename T> 00131 Image<T>::~Image() 00132 { 00133 } 00134 00135 00136 template <typename T> 00137 Image<T> & Image<T>::operator=(const Image<T> &right) 00138 { 00139 // all stack allocated. 00140 m_isRead = right.m_isRead; 00141 m_usingNullVal = right.m_usingNullVal, 00142 m_lastNullVal = right.m_lastNullVal, 00143 m_fullImageCache.resize(right.m_fullImageCache.size()); 00144 m_fullImageCache = right.m_fullImageCache; 00145 m_currentRead.resize(right.m_currentRead.size()); 00146 m_currentRead = right.m_currentRead; 00147 return *this; 00148 } 00149 00150 00151 template <typename T> 00152 const std::valarray<T>& Image<T>::readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls) 00153 { 00154 if (!naxes.size()) 00155 { 00156 m_currentRead.resize(0); 00157 return m_currentRead; 00158 } 00159 unsigned long init(1); 00160 unsigned long nTotalElements(std::accumulate(naxes.begin(),naxes.end(),init, 00161 std::multiplies<long>())); 00162 00163 if (first <= 0) 00164 { 00165 string errMsg("*** CCfits Error: For image read, lowest allowed value for first element is 1.\n"); 00166 bool silent = false; 00167 throw FitsException(errMsg, silent); 00168 } 00169 // 0-based index for slice 00170 unsigned long start = (unsigned long)first - 1; 00171 if (start >= nTotalElements) 00172 { 00173 string errMsg("*** CCfits Error: For image read, starting element is out of range.\n"); 00174 bool silent = false; 00175 throw FitsException(errMsg, silent); 00176 } 00177 if (nElements < 0) 00178 { 00179 string errMsg("*** CCfits Error: Negative nElements value specified for image read.\n"); 00180 bool silent = false; 00181 throw FitsException(errMsg, silent); 00182 } 00183 const unsigned long elementsRequested = (unsigned long)nElements; 00184 00185 int status(0); 00186 int any (0); 00187 FITSUtil::MatchType<T> imageType; 00188 00189 // truncate to valid array size if too much data asked for. 00190 unsigned long elementsToRead(std::min(elementsRequested, 00191 nTotalElements - start)); 00192 if ( elementsToRead < elementsRequested) 00193 { 00194 std::cerr << 00195 "***CCfits Warning: data request exceeds image size, truncating\n"; 00196 } 00197 const bool isFullRead = (elementsToRead == nTotalElements); 00198 const bool isDifferentNull = isNullValChanged(nullValue); 00199 if (!m_isRead || isDifferentNull) 00200 { 00201 // Must perform a read from disk. 00202 m_isRead = false; 00203 if (isFullRead) 00204 { 00205 m_fullImageCache.resize(elementsToRead); 00206 if (fits_read_img(fPtr,imageType(),first,elementsToRead, 00207 nullValue,&m_fullImageCache[0],&any,&status) != 0) throw FitsError(status); 00208 m_isRead = true; 00209 // For this case only, we'll pass m_fullImageCache back up (to be 00210 // copied into user-supplied array). This spares having to do 00211 // what may be a very large copy into m_currentRead. 00212 } 00213 else 00214 { 00215 m_fullImageCache.resize(0); 00216 m_currentRead.resize(elementsToRead); 00217 if (fits_read_img(fPtr,imageType(),first,elementsToRead, 00218 nullValue,&m_currentRead[0],&any,&status) != 0) throw FitsError(status); 00219 } 00220 nulls = (any != 0); 00221 setLastNullInfo(nullValue); 00222 } 00223 else 00224 { 00225 if (!isFullRead) 00226 { 00227 m_currentRead.resize((size_t)elementsToRead); 00228 // This may be a costly copy, though should still be faster 00229 // than disk read. 00230 m_currentRead = m_fullImageCache[std::slice((size_t)start, (size_t)elementsToRead,1)]; 00231 } 00232 } 00233 if (isFullRead) 00234 return m_fullImageCache; 00235 00236 return m_currentRead; 00237 00238 } 00239 00240 template <typename T> 00241 const std::valarray<T>& Image<T>::readImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, T* nullValue, const std::vector<long>& naxes, bool& nulls) 00242 { 00243 const size_t N = naxes.size(); 00244 if (!N) 00245 { 00246 m_currentRead.resize(0); 00247 return m_currentRead; 00248 } 00249 if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size()) 00250 { 00251 string errMsg("*** CCfits Error: Image read function requires that naxes, firstVertex,"); 00252 errMsg += " \nlastVertex, and stride vectors all be the same size.\n"; 00253 bool silent = false; 00254 throw FitsException(errMsg, silent); 00255 } 00256 00257 FITSUtil::CVarray<long> carray; 00258 int any(0); 00259 int status(0); 00260 long requestedSize=1; 00261 long nTotalSize=1; 00262 00263 for (size_t j = 0; j < N; ++j) 00264 { 00265 // Intentional truncation during division. 00266 requestedSize *= ((lastVertex[j] - firstVertex[j])/stride[j] + 1); 00267 nTotalSize *= naxes[j]; 00268 if (firstVertex[j] < 1 || lastVertex[j] > naxes[j]) 00269 { 00270 string errMsg("*** CCfits Error: Out-of-bounds vertex value.\n"); 00271 bool silent=false; 00272 throw FitsException(errMsg,silent); 00273 } 00274 if (firstVertex[j] > lastVertex[j]) 00275 { 00276 string errMsg("*** CCfits Error: firstVertex values must not be larger than corresponding lastVertex values.\n"); 00277 bool silent = false; 00278 throw FitsException(errMsg,silent); 00279 } 00280 } 00281 const bool isFullRead = (requestedSize == nTotalSize); 00282 const bool isDifferentNull = isNullValChanged(nullValue); 00283 if (!m_isRead || isDifferentNull) 00284 { 00285 // Must perform a read from disk. 00286 FITSUtil::auto_array_ptr<long> pFpixel(carray(firstVertex)); 00287 FITSUtil::auto_array_ptr<long> pLpixel(carray(lastVertex)); 00288 FITSUtil::auto_array_ptr<long> pStride(carray(stride)); 00289 00290 FITSUtil::MatchType<T> imageType; 00291 m_isRead = false; 00292 if (isFullRead) 00293 { 00294 m_fullImageCache.resize(requestedSize); 00295 if (fits_read_subset(fPtr,imageType(), 00296 pFpixel.get(),pLpixel.get(), 00297 pStride.get(),nullValue,&m_fullImageCache[0],&any,&status) != 0) 00298 throw FitsError(status); 00299 m_isRead = true; 00300 } 00301 else 00302 { 00303 m_currentRead.resize(requestedSize); 00304 if (fits_read_subset(fPtr,imageType(), 00305 pFpixel.get(),pLpixel.get(), 00306 pStride.get(),nullValue,&m_currentRead[0],&any,&status) != 0) 00307 throw FitsError(status); 00308 } 00309 nulls = (any != 0); 00310 setLastNullInfo(nullValue); 00311 } 00312 else 00313 { 00314 if (!isFullRead) 00315 { 00316 // Must convert firstVertex,lastVertex,stride to gslice parameters. 00317 // Note that in cfitsio, the NAXIS1 dimension varies the fastest 00318 // when laid out in an array in memory (ie. Fortran style). Therefore NAXISn 00319 // ordering must be reversed to C style before passing to gslice. 00320 size_t startPos=0; 00321 std::valarray<size_t> gsliceLength(size_t(0),N); 00322 std::valarray<size_t> gsliceStride(size_t(0),N); 00323 00324 std::vector<long> naxesProducts(N); 00325 long accum=1; 00326 for (size_t i=0; i<N; ++i) 00327 { 00328 naxesProducts[i] = accum; 00329 accum *= naxes[i]; 00330 } 00331 00332 for (size_t i=0; i<N; ++i) 00333 { 00334 startPos += static_cast<size_t>((firstVertex[i]-1)*naxesProducts[i]); 00335 // Here's where we reverse the order: 00336 const size_t gsPos = N-1-i; 00337 // Division truncation is intentional. 00338 gsliceLength[gsPos] = static_cast<size_t>((1 + (lastVertex[i]-firstVertex[i])/stride[i])); 00339 gsliceStride[gsPos] = static_cast<size_t>(stride[i]*naxesProducts[i]); 00340 } 00341 m_currentRead.resize(requestedSize); 00342 m_currentRead = m_fullImageCache[std::gslice(startPos, gsliceLength, gsliceStride)]; 00343 } 00344 } 00345 if (isFullRead) 00346 return m_fullImageCache; 00347 return m_currentRead; 00348 } 00349 00350 template <typename T> 00351 void Image<T>::writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN, T* nullValue) 00352 { 00353 int status(0); 00354 if (first < 1 || nElements < 1) 00355 { 00356 string errMsg("*** CCfits Error: first and nElements values must be > 0\n"); 00357 bool silent = false; 00358 throw FitsException(errMsg, silent); 00359 } 00360 FITSUtil::CAarray<T> convert; 00361 FITSUtil::auto_array_ptr<T> pArray(convert(inData)); 00362 T* array = pArray.get(); 00363 00364 m_isRead = false; 00365 newNaxisN = 0; 00366 00367 FITSUtil::MatchType<T> imageType; 00368 long type(imageType()); 00369 00370 if (fits_write_imgnull(fPtr,type,first,nElements,array, 00371 nullValue,&status)!= 0) 00372 { 00373 throw FitsError(status); 00374 } 00375 const size_t nDim=naxes.size(); 00376 long origTotSize=1; 00377 for (size_t i=0; i<nDim; ++i) 00378 origTotSize *= naxes[i]; 00379 const long highestOutputElem = first + nElements - 1; 00380 if (highestOutputElem > origTotSize) 00381 { 00382 // NAXIS(nDIM) may have increased. 00383 std::ostringstream oss; 00384 oss <<"NAXIS" << nDim; 00385 string keyname(oss.str()); 00386 long newVal = 1 + (highestOutputElem-1)/(origTotSize/naxes[nDim-1]); 00387 if (newVal != naxes[nDim-1]) 00388 { 00389 if (fits_update_key(fPtr,TLONG,(char *)keyname.c_str(),&newVal,0,&status) != 0) 00390 { 00391 throw FitsError(status); 00392 } 00393 newNaxisN = newVal; 00394 } 00395 } 00396 if (fits_flush_file(fPtr,&status) != 0) 00397 throw FitsError(status); 00398 00399 } 00400 00401 template <typename T> 00402 void Image<T>::writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN) 00403 { 00404 // input vectors' size equality will be verified in prepareForSubset. 00405 const size_t nDim = naxes.size(); 00406 FITSUtil::auto_array_ptr<long> pFPixel(new long[nDim]); 00407 FITSUtil::auto_array_ptr<long> pLPixel(new long[nDim]); 00408 std::valarray<T> subset; 00409 m_isRead = false; 00410 newNaxisN = 0; 00411 prepareForSubset(naxes,firstVertex,lastVertex,stride,inData,subset); 00412 00413 long* fPixel = pFPixel.get(); 00414 long* lPixel = pLPixel.get(); 00415 for (size_t i=0; i<nDim; ++i) 00416 { 00417 fPixel[i] = firstVertex[i]; 00418 lPixel[i] = lastVertex[i]; 00419 } 00420 00421 FITSUtil::CAarray<T> convert; 00422 FITSUtil::auto_array_ptr<T> pArray(convert(subset)); 00423 T* array = pArray.get(); 00424 FITSUtil::MatchType<T> imageType; 00425 int status(0); 00426 00427 if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status) ) 00428 throw FitsError(status); 00429 00430 if (lPixel[nDim-1] > naxes[nDim-1]) 00431 { 00432 std::ostringstream oss; 00433 oss << "NAXIS" << nDim; 00434 string keyname(oss.str()); 00435 long newVal = lPixel[nDim-1]; 00436 if (fits_update_key(fPtr,TLONG,(char *)keyname.c_str(),&newVal,0,&status) != 0) 00437 { 00438 throw FitsError(status); 00439 } 00440 newNaxisN = lPixel[nDim-1]; 00441 } 00442 if (fits_flush_file(fPtr,&status) != 0) 00443 throw FitsError(status); 00444 00445 } 00446 00447 template <typename T> 00448 std::valarray<T>& Image<T>::image () 00449 { 00450 00451 return m_fullImageCache; 00452 } 00453 00454 template <typename T> 00455 void Image<T>::prepareForSubset (const std::vector<long>& naxes, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, const std::valarray<T>& inData, std::valarray<T>& subset) 00456 { 00457 00458 // naxes, firstVertex, lastVertex, and stride must all be the same size. 00459 const size_t N = naxes.size(); 00460 if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size()) 00461 { 00462 string errMsg("*** CCfits Error: Image write function requires that naxes, firstVertex,"); 00463 errMsg += " \nlastVertex, and stride vectors all be the same size.\n"; 00464 bool silent = false; 00465 throw FitsException(errMsg, silent); 00466 } 00467 for (size_t i=0; i<N; ++i) 00468 { 00469 if (naxes[i] < 1) 00470 { 00471 bool silent = false; 00472 throw FitsException("*** CCfits Error: Invalid naxes value sent to image write function.\n", silent); 00473 } 00474 string rangeErrMsg("*** CCfits Error: Out-of-range value sent to image write function in arg: "); 00475 if (firstVertex[i] < 1 || (firstVertex[i] > naxes[i] && i != N-1)) 00476 { 00477 bool silent = false; 00478 rangeErrMsg += "firstVertex\n"; 00479 throw FitsException(rangeErrMsg, silent); 00480 } 00481 if (lastVertex[i] < firstVertex[i] || (lastVertex[i] > naxes[i] && i != N-1)) 00482 { 00483 bool silent = false; 00484 rangeErrMsg += "lastVertex\n"; 00485 throw FitsException(rangeErrMsg, silent); 00486 } 00487 if (stride[i] < 1) 00488 { 00489 bool silent = false; 00490 rangeErrMsg += "stride\n"; 00491 throw FitsException(rangeErrMsg, silent); 00492 } 00493 } 00494 00495 // nPoints refers to the subset of the image INCLUDING the zero'ed elements 00496 // resulting from the stride parameter. 00497 // subSizeWithStride refers to the same subset, not counting the zeros. 00498 size_t subSizeWithStride = 1; 00499 size_t nPoints = 1; 00500 std::vector<size_t> subIncr(N); 00501 for (size_t i=0; i<N; ++i) 00502 { 00503 subIncr[i] = nPoints; 00504 nPoints *= static_cast<size_t>(1+lastVertex[i]-firstVertex[i]); 00505 subSizeWithStride *= static_cast<size_t>(1+(lastVertex[i]-firstVertex[i])/stride[i]); 00506 } 00507 subset.resize(nPoints, 0); 00508 00509 if (subSizeWithStride != inData.size()) 00510 { 00511 bool silent = false; 00512 string errMsg("*** CCfits Error: Data array size is not consistent with the values"); 00513 errMsg += "\n in range and stride vectors sent to the image write function.\n"; 00514 throw FitsException(errMsg, silent); 00515 } 00516 00517 size_t startPoint = 0; 00518 size_t dimMult = 1; 00519 std::vector<size_t> incr(N); 00520 for (size_t j = 0; j < N; ++j) 00521 { 00522 startPoint += dimMult*(firstVertex[j]-1); 00523 incr[j] = dimMult; 00524 dimMult *= static_cast<size_t>(naxes[j]); 00525 } 00526 00527 size_t inDataPos = 0; 00528 size_t iSub = 0; 00529 loop(N-1, firstVertex, lastVertex, stride, startPoint, incr, inData, inDataPos, subIncr, subset, iSub); 00530 } 00531 00532 template <typename T> 00533 void Image<T>::loop (size_t iDim, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::vector<long>& stride, size_t iPos, const std::vector<size_t>& incr, const std::valarray<T>& inData, size_t& iDat, const std::vector<size_t>& subIncr, std::valarray<T>& subset, size_t iSub) 00534 { 00535 size_t start = static_cast<size_t>(firstVertex[iDim]); 00536 size_t stop = static_cast<size_t>(lastVertex[iDim]); 00537 size_t skip = static_cast<size_t>(stride[iDim]); 00538 if (iDim == 0) 00539 { 00540 size_t length = stop - start + 1; 00541 for (size_t i=0; i<length; i+=skip) 00542 { 00543 subset[i+iSub] = inData[iDat++]; 00544 } 00545 } 00546 else 00547 { 00548 size_t jump = incr[iDim]*skip; 00549 size_t subJump = subIncr[iDim]*skip; 00550 for (size_t i=start; i<=stop; i+=skip) 00551 { 00552 loop(iDim-1, firstVertex, lastVertex, stride, iPos, incr, inData, iDat, subIncr, subset, iSub); 00553 iPos += jump; 00554 iSub += subJump; 00555 } 00556 } 00557 } 00558 00559 template <typename T> 00560 bool Image<T>::isNullValChanged(T* newNull) const 00561 { 00562 bool isChanged = false; 00563 if (m_usingNullVal) 00564 { 00565 // If m_usingNullVal is true, we can assume m_lastNullVal != 0. 00566 if (newNull) 00567 { 00568 T newVal = *newNull; 00569 if (newVal != m_lastNullVal) 00570 isChanged = true; 00571 } 00572 else 00573 isChanged = true; 00574 } 00575 else 00576 { 00577 if (newNull && (*newNull != 0)) 00578 isChanged = true; 00579 } 00580 00581 return isChanged; 00582 } 00583 00584 template <typename T> 00585 void Image<T>::setLastNullInfo(T* newNull) 00586 { 00587 if (!newNull || *newNull == 0) 00588 { 00589 m_usingNullVal = false; 00590 m_lastNullVal = 0; 00591 } 00592 else 00593 { 00594 m_usingNullVal = true; 00595 m_lastNullVal = *newNull; 00596 } 00597 } 00598 00599 template <typename T> 00600 void Image<T>::writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN) 00601 { 00602 std::vector<long> stride(firstVertex.size(), 1); 00603 writeImage(fPtr, firstVertex, lastVertex, stride, inData, naxes, newNaxisN); 00604 } 00605 00606 template <typename T> 00607 void Image<T>::scalingHasChanged() 00608 { 00609 m_isRead = false; 00610 } 00611 00612 template <typename T> 00613 void Image<T>::resetRead() 00614 { 00615 m_isRead = false; 00616 } 00617 00618 // Additional Declarations 00619 00620 } // namespace CCfits 00621 00622 00623 #endif