CCfits
2.4
|
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 #ifdef _MSC_VER 00021 #include "MSconfig.h" //form std::min 00022 #endif 00023 #include "CCfits.h" 00024 #include "FitsError.h" 00025 #include "FITSUtil.h" 00026 00027 00028 namespace CCfits { 00029 00030 00031 00032 template <typename T> 00033 class Image 00034 { 00035 00036 public: 00037 Image(const Image< T > &right); 00038 Image (const std::valarray<T>& imageArray = std::valarray<T>()); 00039 ~Image(); 00040 Image< T > & operator=(const Image< T > &right); 00041 00042 // Read data reads the image if readFlag is true and 00043 // optional keywords if supplied. Thus, with no arguments, 00044 // readData() does nothing. 00045 const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls); 00046 // Read data reads the image if readFlag is true and 00047 // optional keywords if supplied. Thus, with no arguments, 00048 // readData() does nothing. 00049 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); 00050 // Read data reads the image if readFlag is true and 00051 // optional keywords if supplied. Thus, with no arguments, 00052 // readData() does nothing. 00053 void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue = 0); 00054 // Read data reads the image if readFlag is true and 00055 // optional keywords if supplied. Thus, with no arguments, 00056 // readData() does nothing. 00057 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); 00058 void writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::valarray<T>& inData, const std::vector<long>& naxes); 00059 bool isRead () const; 00060 void isRead (bool value); 00061 const std::valarray< T >& image () const; 00062 void setImage (const std::valarray< T >& value); 00063 const T image (size_t index) const; 00064 void setImage (size_t index, T value); 00065 00066 // Additional Public Declarations 00067 00068 protected: 00069 // Additional Protected Declarations 00070 00071 private: 00072 std::valarray<T>& image (); 00073 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); 00074 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); 00075 00076 // Additional Private Declarations 00077 00078 private: //## implementation 00079 // Data Members for Class Attributes 00080 bool m_isRead; 00081 00082 // Data Members for Associations 00083 std::valarray< T > m_image; 00084 00085 // Additional Implementation Declarations 00086 00087 }; 00088 00089 // Parameterized Class CCfits::Image 00090 00091 template <typename T> 00092 inline bool Image<T>::isRead () const 00093 { 00094 return m_isRead; 00095 } 00096 00097 template <typename T> 00098 inline void Image<T>::isRead (bool value) 00099 { 00100 m_isRead = value; 00101 } 00102 00103 template <typename T> 00104 inline const std::valarray< T >& Image<T>::image () const 00105 { 00106 return m_image; 00107 } 00108 00109 template <typename T> 00110 inline void Image<T>::setImage (const std::valarray< T >& value) 00111 { 00112 m_image.resize(value.size()); 00113 m_image = value; 00114 } 00115 00116 template <typename T> 00117 inline const T Image<T>::image (size_t index) const 00118 { 00119 return m_image[index]; 00120 } 00121 00122 template <typename T> 00123 inline void Image<T>::setImage (size_t index, T value) 00124 { 00125 m_image[index] = value; 00126 } 00127 00128 // Parameterized Class CCfits::Image 00129 00130 template <typename T> 00131 Image<T>::Image(const Image<T> &right) 00132 : m_isRead(right.m_isRead), 00133 m_image(right.m_image) 00134 { 00135 } 00136 00137 template <typename T> 00138 Image<T>::Image (const std::valarray<T>& imageArray) 00139 : m_isRead(false), 00140 m_image(imageArray) 00141 { 00142 } 00143 00144 00145 template <typename T> 00146 Image<T>::~Image() 00147 { 00148 } 00149 00150 00151 template <typename T> 00152 Image<T> & Image<T>::operator=(const Image<T> &right) 00153 { 00154 // all stack allocated. 00155 m_isRead = right.m_isRead; 00156 m_image.resize(right.m_image.size()); 00157 m_image = right.m_image; 00158 return *this; 00159 } 00160 00161 00162 template <typename T> 00163 const std::valarray<T>& Image<T>::readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls) 00164 { 00165 const size_t N(naxes.size()); 00166 if (N > 0) 00167 { 00168 int status(0); 00169 int any (0); 00170 FITSUtil::MatchType<T> imageType; 00171 unsigned long init(1); 00172 unsigned long nelements(std::accumulate(naxes.begin(),naxes.end(),init, 00173 std::multiplies<long>())); 00174 00175 // truncate to valid array size if too much data asked for. 00176 // note first is 1-based index) 00177 long elementsToRead(std::min(static_cast<unsigned long>(nElements), 00178 nelements - first + 1)); 00179 if ( elementsToRead < nElements) 00180 { 00181 std::cerr << 00182 "***CCfits Warning: data request exceeds image size, truncating\n"; 00183 } 00184 FITSUtil::FitsNullValue<T> null; 00185 // initialize m_image to nullValue. resize if necessary. 00186 if (m_image.size() != static_cast<size_t>(elementsToRead)) 00187 { 00188 m_image.resize(elementsToRead,null()); 00189 } 00190 if (fits_read_img(fPtr,imageType(),first,elementsToRead, 00191 nullValue,&m_image[0],&any,&status) != 0) throw FitsError(status); 00192 00193 nulls = (any != 0); 00194 m_isRead = (first == 1 && nelements == static_cast<unsigned long>(nElements)); 00195 } 00196 else 00197 { 00198 m_isRead = true; 00199 m_image.resize(0); 00200 } 00201 return m_image; 00202 } 00203 00204 template <typename T> 00205 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) 00206 { 00207 00208 00209 00210 FITSUtil::CVarray<long> carray; 00211 00212 00213 int any(0); 00214 int status(0); 00215 const size_t N(naxes.size()); 00216 00217 size_t arraySize(1); 00218 00219 for (size_t j = 0; j < N; ++j) 00220 { 00221 arraySize *= (lastVertex[j] - firstVertex[j] + 1); 00222 } 00223 00224 FITSUtil::auto_array_ptr<long> pFpixel(carray(firstVertex)); 00225 FITSUtil::auto_array_ptr<long> pLpixel(carray(lastVertex)); 00226 FITSUtil::auto_array_ptr<long> pStride(carray(stride)); 00227 00228 FITSUtil::MatchType<T> imageType; 00229 00230 size_t n(m_image.size()); 00231 if (n != arraySize) m_image.resize(arraySize); 00232 if (fits_read_subset(fPtr,imageType(), 00233 pFpixel.get(),pLpixel.get(), 00234 pStride.get(),nullValue,&m_image[0],&any,&status) != 0) 00235 { 00236 throw FitsError(status); 00237 00238 } 00239 00240 nulls = (any != 0); 00241 return m_image; 00242 } 00243 00244 template <typename T> 00245 void Image<T>::writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue) 00246 { 00247 00248 00249 int status(0); 00250 size_t init(1); 00251 size_t totalSize= static_cast<size_t>(std::accumulate(naxes.begin(),naxes.end(),init,std::multiplies<long>() )); 00252 FITSUtil::FitsNullValue<T> null; 00253 if (m_image.size() != totalSize) m_image.resize(totalSize,null()); 00254 FITSUtil::CAarray<T> convert; 00255 FITSUtil::auto_array_ptr<T> pArray(convert(inData)); 00256 T* array = pArray.get(); 00257 00258 00259 FITSUtil::MatchType<T> imageType; 00260 long type(imageType()); 00261 00262 if (fits_write_imgnull(fPtr,type,first,nElements,array, 00263 nullValue,&status) || fits_flush_file(fPtr,&status) != 0) 00264 { 00265 throw FitsError(status); 00266 00267 } 00268 00269 00270 00271 m_image[std::slice(first-1,nElements,1)] = inData; 00272 } 00273 00274 template <typename T> 00275 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) 00276 { 00277 // input vectors' size equality will be verified in prepareForSubset. 00278 const size_t nDim = naxes.size(); 00279 FITSUtil::auto_array_ptr<long> pFPixel(new long[nDim]); 00280 FITSUtil::auto_array_ptr<long> pLPixel(new long[nDim]); 00281 std::valarray<T> subset; 00282 prepareForSubset(naxes,firstVertex,lastVertex,stride,inData,subset); 00283 00284 long* fPixel = pFPixel.get(); 00285 long* lPixel = pLPixel.get(); 00286 for (size_t i=0; i<nDim; ++i) 00287 { 00288 fPixel[i] = firstVertex[i]; 00289 lPixel[i] = lastVertex[i]; 00290 } 00291 00292 FITSUtil::CAarray<T> convert; 00293 FITSUtil::auto_array_ptr<T> pArray(convert(subset)); 00294 T* array = pArray.get(); 00295 FITSUtil::MatchType<T> imageType; 00296 int status(0); 00297 00298 if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status) 00299 || fits_flush_file(fPtr,&status) != 0) throw FitsError(status); 00300 } 00301 00302 template <typename T> 00303 std::valarray<T>& Image<T>::image () 00304 { 00305 00306 return m_image; 00307 } 00308 00309 template <typename T> 00310 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) 00311 { 00312 00313 // naxes, firstVertex, lastVertex, and stride must all be the same size. 00314 const size_t N = naxes.size(); 00315 if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size()) 00316 { 00317 string errMsg("*** CCfits Error: Image write function requires that naxes, firstVertex,"); 00318 errMsg += " \nlastVertex, and stride vectors all be the same size.\n"; 00319 bool silent = false; 00320 throw FitsException(errMsg, silent); 00321 } 00322 for (size_t i=0; i<N; ++i) 00323 { 00324 if (naxes[i] < 1) 00325 { 00326 bool silent = false; 00327 throw FitsException("*** CCfits Error: Invalid naxes value sent to image write function.\n", silent); 00328 } 00329 string rangeErrMsg("*** CCfits Error: Out-of-range value sent to image write function in arg: "); 00330 if (firstVertex[i] < 1 || firstVertex[i] > naxes[i]) 00331 { 00332 bool silent = false; 00333 rangeErrMsg += "firstVertex\n"; 00334 throw FitsException(rangeErrMsg, silent); 00335 } 00336 if (lastVertex[i] < firstVertex[i] || lastVertex[i] > naxes[i]) 00337 { 00338 bool silent = false; 00339 rangeErrMsg += "lastVertex\n"; 00340 throw FitsException(rangeErrMsg, silent); 00341 } 00342 if (stride[i] < 1) 00343 { 00344 bool silent = false; 00345 rangeErrMsg += "stride\n"; 00346 throw FitsException(rangeErrMsg, silent); 00347 } 00348 } 00349 00350 // nPoints refers to the subset of m_image INCLUDING the zero'ed elements 00351 // resulting from the stride parameter. 00352 // subSizeWithStride refers to the same subset, not counting the zeros. 00353 size_t subSizeWithStride = 1; 00354 size_t nPoints = 1; 00355 std::vector<size_t> subIncr(N); 00356 for (size_t i=0; i<N; ++i) 00357 { 00358 subIncr[i] = nPoints; 00359 nPoints *= static_cast<size_t>(1+lastVertex[i]-firstVertex[i]); 00360 subSizeWithStride *= static_cast<size_t>(1+(lastVertex[i]-firstVertex[i])/stride[i]); 00361 } 00362 FITSUtil::FitsNullValue<T> null; 00363 subset.resize(nPoints, null()); 00364 00365 // Trying to avoid at all costs an assignment between 2 valarrays of 00366 // different sizes when m_image gets set below. 00367 if (subSizeWithStride != inData.size()) 00368 { 00369 bool silent = false; 00370 string errMsg("*** CCfits Error: Data array size is not consistent with the values"); 00371 errMsg += "\n in range and stride vectors sent to the image write function.\n"; 00372 throw FitsException(errMsg, silent); 00373 } 00374 00375 size_t startPoint = 0; 00376 size_t dimMult = 1; 00377 std::vector<size_t> incr(N); 00378 for (size_t j = 0; j < N; ++j) 00379 { 00380 startPoint += dimMult*(firstVertex[j]-1); 00381 incr[j] = dimMult; 00382 dimMult *= static_cast<size_t>(naxes[j]); 00383 } 00384 const size_t imageSize = dimMult; 00385 m_image.resize(imageSize,null()); 00386 00387 size_t inDataPos = 0; 00388 size_t iSub = 0; 00389 loop(N-1, firstVertex, lastVertex, stride, startPoint, incr, inData, inDataPos, subIncr, subset, iSub); 00390 } 00391 00392 template <typename T> 00393 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) 00394 { 00395 size_t start = static_cast<size_t>(firstVertex[iDim]); 00396 size_t stop = static_cast<size_t>(lastVertex[iDim]); 00397 size_t skip = static_cast<size_t>(stride[iDim]); 00398 if (iDim == 0) 00399 { 00400 size_t length = stop - start + 1; 00401 for (size_t i=0; i<length; i+=skip) 00402 { 00403 m_image[i+iPos] = inData[iDat]; 00404 subset[i+iSub] = inData[iDat++]; 00405 } 00406 } 00407 else 00408 { 00409 size_t jump = incr[iDim]*skip; 00410 size_t subJump = subIncr[iDim]*skip; 00411 for (size_t i=start; i<=stop; i+=skip) 00412 { 00413 loop(iDim-1, firstVertex, lastVertex, stride, iPos, incr, inData, iDat, subIncr, subset, iSub); 00414 iPos += jump; 00415 iSub += subJump; 00416 } 00417 } 00418 } 00419 00420 template <typename T> 00421 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) 00422 { 00423 std::vector<long> stride(firstVertex.size(), 1); 00424 writeImage(fPtr, firstVertex, lastVertex, stride, inData, naxes); 00425 } 00426 00427 // Additional Declarations 00428 00429 } // namespace CCfits 00430 00431 00432 #endif