Image.h

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

Generated on Fri Nov 3 17:09:05 2006 for CCfits by  doxygen 1.4.7