00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef IMAGE_H
00010 #define IMAGE_H 1
00011
00012
00013 #include <functional>
00014
00015 #include <valarray>
00016
00017 #include <vector>
00018
00019 #include <numeric>
00020 #ifdef _MSC_VER
00021 #include "MSconfig.h"
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
00043
00044
00045 const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls);
00046
00047
00048
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
00051
00052
00053 void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue = 0);
00054
00055
00056
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
00067
00068 protected:
00069
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
00077
00078 private:
00079
00080 bool m_isRead;
00081
00082
00083 std::valarray< T > m_image;
00084
00085
00086
00087 };
00088
00089
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
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
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
00176
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
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
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
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
00351
00352
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
00366
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
00428
00429 }
00430
00431
00432 #endif