Image.h

00001 //   Read the documentation to learn more about C++ code generator
00002 //   versioning.
00003 //      This is version 2.1 release dated Nov 2008
00004 //      Astrophysics Science Division,
00005 //      NASA/ Goddard Space Flight Center
00006 //      HEASARC
00007 //      http://heasarc.gsfc.nasa.gov
00008 //      e-mail: ccfits@legacy.gsfc.nasa.gov
00009 //
00010 //      Original author: Ben Dorman
00011 
00012 #ifndef IMAGE_H
00013 #define IMAGE_H 1
00014 
00015 // functional
00016 #include <functional>
00017 // vector
00018 #include <vector>
00019 // valarray
00020 #include <valarray>
00021 // numeric
00022 #include <numeric>
00023 #ifdef _MSC_VER
00024 #include "MSconfig.h" //form std::min
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         //      Read data reads the image if readFlag is true and
00046         //      optional keywords if supplied. Thus, with no arguments,
00047         //      readData() does nothing.
00048         const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls);
00049         //      Read data reads the image if readFlag is true and
00050         //      optional keywords if supplied. Thus, with no arguments,
00051         //      readData() does nothing.
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         //      Read data reads the image if readFlag is true and
00054         //      optional keywords if supplied. Thus, with no arguments,
00055         //      readData() does nothing.
00056         void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue = 0);
00057         //      Read data reads the image if readFlag is true and
00058         //      optional keywords if supplied. Thus, with no arguments,
00059         //      readData() does nothing.
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       // Additional Public Declarations
00070 
00071     protected:
00072       // Additional Protected Declarations
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       // Additional Private Declarations
00080 
00081     private: //## implementation
00082       // Data Members for Class Attributes
00083         bool m_isRead;
00084 
00085       // Data Members for Associations
00086         std::valarray< T > m_image;
00087 
00088       // Additional Implementation Declarations
00089 
00090   };
00091 
00092   // Parameterized Class CCfits::Image 
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   // Parameterized Class CCfits::Image 
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       // all stack allocated.
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                 // truncate to valid array size if too much data asked for.
00179                 // note first is 1-based index)
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                 // initialize m_image to nullValue. resize if necessary.
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      //size_t imageSize(std::accumulate(&naxes[0],&naxes[N],init,std::multiplies<T>() ));
00250      // m_isRead = (startPoint  == 0 && imageSize == static_cast<unsigned long>(nElements)); 
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         // input vectors' size equality will be verified in prepareForSubset.
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     // naxes, firstVertex, lastVertex, and stride must all be the same size.
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     // nPoints refers to the subset of m_image INCLUDING the zero'ed elements 
00363     // resulting from the stride parameter.  
00364     // subSizeWithStride refers to the same subset, not counting the zeros.
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     // Trying to avoid at all costs an assignment between 2 valarrays of 
00378     // different sizes when m_image gets set below.
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   // Additional Declarations
00440 
00441 } // namespace CCfits
00442 
00443 
00444 #endif

Generated on Thu Nov 6 10:44:39 2008 for CCfits by  doxygen 1.5.4