00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef IMAGE_H
00013 #define IMAGE_H 1
00014
00015
00016 #include <vector>
00017
00018 #include <valarray>
00019
00020 #include <functional>
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[0],&naxes[N],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[0],&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 const size_t N(naxes.size());
00262 size_t init(1);
00263 size_t totalSize= static_cast<size_t>(std::accumulate(&naxes[0],&naxes[N],init,std::multiplies<long>() ));
00264 FITSUtil::FitsNullValue<T> null;
00265 if (m_image.size() != totalSize) m_image.resize(totalSize,null());
00266 FITSUtil::CAarray<T> convert;
00267 FITSUtil::auto_array_ptr<T> pArray(convert(inData));
00268 T* array = pArray.get();
00269
00270
00271 FITSUtil::MatchType<T> imageType;
00272 long type(imageType());
00273
00274 if (fits_write_imgnull(fPtr,type,first,nElements,array,
00275 nullValue,&status) || fits_flush_file(fPtr,&status) != 0)
00276 {
00277 throw FitsError(status);
00278
00279 }
00280
00281
00282
00283 m_image[std::slice(first-1,nElements,1)] = inData;
00284 }
00285
00286 template <typename T>
00287 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)
00288 {
00289
00290 const size_t nDim = naxes.size();
00291 FITSUtil::auto_array_ptr<long> pFPixel(new long[nDim]);
00292 FITSUtil::auto_array_ptr<long> pLPixel(new long[nDim]);
00293 std::valarray<T> subset;
00294 prepareForSubset(naxes,firstVertex,lastVertex,stride,inData,subset);
00295
00296 long* fPixel = pFPixel.get();
00297 long* lPixel = pLPixel.get();
00298 for (size_t i=0; i<nDim; ++i)
00299 {
00300 fPixel[i] = firstVertex[i];
00301 lPixel[i] = lastVertex[i];
00302 }
00303
00304 FITSUtil::CAarray<T> convert;
00305 FITSUtil::auto_array_ptr<T> pArray(convert(subset));
00306 T* array = pArray.get();
00307 FITSUtil::MatchType<T> imageType;
00308 int status(0);
00309
00310 if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status)
00311 || fits_flush_file(fPtr,&status) != 0) throw FitsError(status);
00312 }
00313
00314 template <typename T>
00315 std::valarray<T>& Image<T>::image ()
00316 {
00317
00318 return m_image;
00319 }
00320
00321 template <typename T>
00322 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)
00323 {
00324
00325
00326 const size_t N = naxes.size();
00327 if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size())
00328 {
00329 string errMsg("*** CCfits Error: Image write function requires that naxes, firstVertex,");
00330 errMsg += " \nlastVertex, and stride vectors all be the same size.\n";
00331 bool silent = false;
00332 throw FitsException(errMsg, silent);
00333 }
00334 for (size_t i=0; i<N; ++i)
00335 {
00336 if (naxes[i] < 1)
00337 {
00338 bool silent = false;
00339 throw FitsException("*** CCfits Error: Invalid naxes value sent to image write function.\n", silent);
00340 }
00341 string rangeErrMsg("*** CCfits Error: Out-of-range value sent to image write function in arg: ");
00342 if (firstVertex[i] < 1 || firstVertex[i] > naxes[i])
00343 {
00344 bool silent = false;
00345 rangeErrMsg += "firstVertex\n";
00346 throw FitsException(rangeErrMsg, silent);
00347 }
00348 if (lastVertex[i] < firstVertex[i] || lastVertex[i] > naxes[i])
00349 {
00350 bool silent = false;
00351 rangeErrMsg += "lastVertex\n";
00352 throw FitsException(rangeErrMsg, silent);
00353 }
00354 if (stride[i] < 1)
00355 {
00356 bool silent = false;
00357 rangeErrMsg += "stride\n";
00358 throw FitsException(rangeErrMsg, silent);
00359 }
00360 }
00361
00362
00363
00364
00365 size_t subSizeWithStride = 1;
00366 size_t nPoints = 1;
00367 std::vector<size_t> subIncr(N);
00368 for (size_t i=0; i<N; ++i)
00369 {
00370 subIncr[i] = nPoints;
00371 nPoints *= static_cast<size_t>(1+lastVertex[i]-firstVertex[i]);
00372 subSizeWithStride *= static_cast<size_t>(1+(lastVertex[i]-firstVertex[i])/stride[i]);
00373 }
00374 FITSUtil::FitsNullValue<T> null;
00375 subset.resize(nPoints, null());
00376
00377
00378
00379 if (subSizeWithStride != inData.size())
00380 {
00381 bool silent = false;
00382 string errMsg("*** CCfits Error: Data array size is not consistent with the values");
00383 errMsg += "\n in range and stride vectors sent to the image write function.\n";
00384 throw FitsException(errMsg, silent);
00385 }
00386
00387 size_t startPoint = 0;
00388 size_t dimMult = 1;
00389 std::vector<size_t> incr(N);
00390 for (size_t j = 0; j < N; ++j)
00391 {
00392 startPoint += dimMult*(firstVertex[j]-1);
00393 incr[j] = dimMult;
00394 dimMult *= static_cast<size_t>(naxes[j]);
00395 }
00396 const size_t imageSize = dimMult;
00397 m_image.resize(imageSize,null());
00398
00399 size_t inDataPos = 0;
00400 size_t iSub = 0;
00401 loop(N-1, firstVertex, lastVertex, stride, startPoint, incr, inData, inDataPos, subIncr, subset, iSub);
00402 }
00403
00404 template <typename T>
00405 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)
00406 {
00407 size_t start = static_cast<size_t>(firstVertex[iDim]);
00408 size_t stop = static_cast<size_t>(lastVertex[iDim]);
00409 size_t skip = static_cast<size_t>(stride[iDim]);
00410 if (iDim == 0)
00411 {
00412 size_t length = stop - start + 1;
00413 for (size_t i=0; i<length; i+=skip)
00414 {
00415 m_image[i+iPos] = inData[iDat];
00416 subset[i+iSub] = inData[iDat++];
00417 }
00418 }
00419 else
00420 {
00421 size_t jump = incr[iDim]*skip;
00422 size_t subJump = subIncr[iDim]*skip;
00423 for (size_t i=start; i<=stop; i+=skip)
00424 {
00425 loop(iDim-1, firstVertex, lastVertex, stride, iPos, incr, inData, iDat, subIncr, subset, iSub);
00426 iPos += jump;
00427 iSub += subJump;
00428 }
00429 }
00430 }
00431
00432 template <typename T>
00433 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)
00434 {
00435 std::vector<long> stride(firstVertex.size(), 1);
00436 writeImage(fPtr, firstVertex, lastVertex, stride, inData, naxes);
00437 }
00438
00439
00440
00441 }
00442
00443
00444 #endif