CCfits  2.4
Image.h
00001 //      Astrophysics Science Division,
00002 //      NASA/ Goddard Space Flight Center
00003 //      HEASARC
00004 //      http://heasarc.gsfc.nasa.gov
00005 //      e-mail: ccfits@legacy.gsfc.nasa.gov
00006 //
00007 //      Original author: Ben Dorman
00008 
00009 #ifndef IMAGE_H
00010 #define IMAGE_H 1
00011 
00012 // functional
00013 #include <functional>
00014 // valarray
00015 #include <valarray>
00016 // vector
00017 #include <vector>
00018 // numeric
00019 #include <numeric>
00020 #ifdef _MSC_VER
00021 #include "MSconfig.h" //form std::min
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         //      Read data reads the image if readFlag is true and
00043         //      optional keywords if supplied. Thus, with no arguments,
00044         //      readData() does nothing.
00045         const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls);
00046         //      Read data reads the image if readFlag is true and
00047         //      optional keywords if supplied. Thus, with no arguments,
00048         //      readData() does nothing.
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         //      Read data reads the image if readFlag is true and
00051         //      optional keywords if supplied. Thus, with no arguments,
00052         //      readData() does nothing.
00053         void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue = 0);
00054         //      Read data reads the image if readFlag is true and
00055         //      optional keywords if supplied. Thus, with no arguments,
00056         //      readData() does nothing.
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       // Additional Public Declarations
00067 
00068     protected:
00069       // Additional Protected Declarations
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       // Additional Private Declarations
00077 
00078     private: //## implementation
00079       // Data Members for Class Attributes
00080         bool m_isRead;
00081 
00082       // Data Members for Associations
00083         std::valarray< T > m_image;
00084 
00085       // Additional Implementation Declarations
00086 
00087   };
00088 
00089   // Parameterized Class CCfits::Image 
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   // Parameterized Class CCfits::Image 
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       // all stack allocated.
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                 // truncate to valid array size if too much data asked for.
00176                 // note first is 1-based index)
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                 // initialize m_image to nullValue. resize if necessary.
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         // input vectors' size equality will be verified in prepareForSubset.
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     // naxes, firstVertex, lastVertex, and stride must all be the same size.
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     // nPoints refers to the subset of m_image INCLUDING the zero'ed elements 
00351     // resulting from the stride parameter.  
00352     // subSizeWithStride refers to the same subset, not counting the zeros.
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     // Trying to avoid at all costs an assignment between 2 valarrays of 
00366     // different sizes when m_image gets set below.
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   // Additional Declarations
00428 
00429 } // namespace CCfits
00430 
00431 
00432 #endif