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