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