CCfits  2.4
Image.h
1 // Astrophysics Science Division,
2 // NASA/ Goddard Space Flight Center
3 // HEASARC
4 // http://heasarc.gsfc.nasa.gov
5 // e-mail: ccfits@legacy.gsfc.nasa.gov
6 //
7 // Original author: Ben Dorman
8 
9 #ifndef IMAGE_H
10 #define IMAGE_H 1
11 
12 // functional
13 #include <functional>
14 // valarray
15 #include <valarray>
16 // vector
17 #include <vector>
18 // numeric
19 #include <numeric>
20 #ifdef _MSC_VER
21 #include "MSconfig.h" //form std::min
22 #endif
23 #include "CCfits.h"
24 #include "FitsError.h"
25 #include "FITSUtil.h"
26 
27 
28 namespace CCfits {
29 
30 
31 
32  template <typename T>
33  class Image
34  {
35 
36  public:
37  Image(const Image< T > &right);
38  Image (const std::valarray<T>& imageArray = std::valarray<T>());
39  ~Image();
40  Image< T > & operator=(const Image< T > &right);
41 
42  // Read data reads the image if readFlag is true and
43  // optional keywords if supplied. Thus, with no arguments,
44  // readData() does nothing.
45  const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls);
46  // Read data reads the image if readFlag is true and
47  // optional keywords if supplied. Thus, with no arguments,
48  // readData() does nothing.
49  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);
50  // Read data reads the image if readFlag is true and
51  // optional keywords if supplied. Thus, with no arguments,
52  // readData() does nothing.
53  void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue = 0);
54  // Read data reads the image if readFlag is true and
55  // optional keywords if supplied. Thus, with no arguments,
56  // readData() does nothing.
57  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);
58  void writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::valarray<T>& inData, const std::vector<long>& naxes);
59  bool isRead () const;
60  void isRead (bool value);
61  const std::valarray< T >& image () const;
62  void setImage (const std::valarray< T >& value);
63  const T image (size_t index) const;
64  void setImage (size_t index, T value);
65 
66  // Additional Public Declarations
67 
68  protected:
69  // Additional Protected Declarations
70 
71  private:
72  std::valarray<T>& image ();
73  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);
74  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);
75 
76  // Additional Private Declarations
77 
78  private: //## implementation
79  // Data Members for Class Attributes
80  bool m_isRead;
81 
82  // Data Members for Associations
83  std::valarray< T > m_image;
84 
85  // Additional Implementation Declarations
86 
87  };
88 
89  // Parameterized Class CCfits::Image
90 
91  template <typename T>
92  inline bool Image<T>::isRead () const
93  {
94  return m_isRead;
95  }
96 
97  template <typename T>
98  inline void Image<T>::isRead (bool value)
99  {
100  m_isRead = value;
101  }
102 
103  template <typename T>
104  inline const std::valarray< T >& Image<T>::image () const
105  {
106  return m_image;
107  }
108 
109  template <typename T>
110  inline void Image<T>::setImage (const std::valarray< T >& value)
111  {
112  m_image.resize(value.size());
113  m_image = value;
114  }
115 
116  template <typename T>
117  inline const T Image<T>::image (size_t index) const
118  {
119  return m_image[index];
120  }
121 
122  template <typename T>
123  inline void Image<T>::setImage (size_t index, T value)
124  {
125  m_image[index] = value;
126  }
127 
128  // Parameterized Class CCfits::Image
129 
130  template <typename T>
131  Image<T>::Image(const Image<T> &right)
132  : m_isRead(right.m_isRead),
133  m_image(right.m_image)
134  {
135  }
136 
137  template <typename T>
138  Image<T>::Image (const std::valarray<T>& imageArray)
139  : m_isRead(false),
140  m_image(imageArray)
141  {
142  }
143 
144 
145  template <typename T>
146  Image<T>::~Image()
147  {
148  }
149 
150 
151  template <typename T>
152  Image<T> & Image<T>::operator=(const Image<T> &right)
153  {
154  // all stack allocated.
155  m_isRead = right.m_isRead;
156  m_image.resize(right.m_image.size());
157  m_image = right.m_image;
158  return *this;
159  }
160 
161 
162  template <typename T>
163  const std::valarray<T>& Image<T>::readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls)
164  {
165  const size_t N(naxes.size());
166  if (N > 0)
167  {
168  int status(0);
169  int any (0);
170  FITSUtil::MatchType<T> imageType;
171  unsigned long init(1);
172  unsigned long nelements(std::accumulate(naxes.begin(),naxes.end(),init,
173  std::multiplies<long>()));
174 
175  // truncate to valid array size if too much data asked for.
176  // note first is 1-based index)
177  long elementsToRead(std::min(static_cast<unsigned long>(nElements),
178  nelements - first + 1));
179  if ( elementsToRead < nElements)
180  {
181  std::cerr <<
182  "***CCfits Warning: data request exceeds image size, truncating\n";
183  }
184  FITSUtil::FitsNullValue<T> null;
185  // initialize m_image to nullValue. resize if necessary.
186  if (m_image.size() != static_cast<size_t>(elementsToRead))
187  {
188  m_image.resize(elementsToRead,null());
189  }
190  if (fits_read_img(fPtr,imageType(),first,elementsToRead,
191  nullValue,&m_image[0],&any,&status) != 0) throw FitsError(status);
192 
193  nulls = (any != 0);
194  m_isRead = (first == 1 && nelements == static_cast<unsigned long>(nElements));
195  }
196  else
197  {
198  m_isRead = true;
199  m_image.resize(0);
200  }
201  return m_image;
202  }
203 
204  template <typename T>
205  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)
206  {
207 
208 
209 
210  FITSUtil::CVarray<long> carray;
211 
212 
213  int any(0);
214  int status(0);
215  const size_t N(naxes.size());
216 
217  size_t arraySize(1);
218 
219  for (size_t j = 0; j < N; ++j)
220  {
221  arraySize *= (lastVertex[j] - firstVertex[j] + 1);
222  }
223 
224  FITSUtil::auto_array_ptr<long> pFpixel(carray(firstVertex));
225  FITSUtil::auto_array_ptr<long> pLpixel(carray(lastVertex));
226  FITSUtil::auto_array_ptr<long> pStride(carray(stride));
227 
228  FITSUtil::MatchType<T> imageType;
229 
230  size_t n(m_image.size());
231  if (n != arraySize) m_image.resize(arraySize);
232  if (fits_read_subset(fPtr,imageType(),
233  pFpixel.get(),pLpixel.get(),
234  pStride.get(),nullValue,&m_image[0],&any,&status) != 0)
235  {
236  throw FitsError(status);
237 
238  }
239 
240  nulls = (any != 0);
241  return m_image;
242  }
243 
244  template <typename T>
245  void Image<T>::writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, T* nullValue)
246  {
247 
248 
249  int status(0);
250  size_t init(1);
251  size_t totalSize= static_cast<size_t>(std::accumulate(naxes.begin(),naxes.end(),init,std::multiplies<long>() ));
252  FITSUtil::FitsNullValue<T> null;
253  if (m_image.size() != totalSize) m_image.resize(totalSize,null());
254  FITSUtil::CAarray<T> convert;
255  FITSUtil::auto_array_ptr<T> pArray(convert(inData));
256  T* array = pArray.get();
257 
258 
259  FITSUtil::MatchType<T> imageType;
260  long type(imageType());
261 
262  if (fits_write_imgnull(fPtr,type,first,nElements,array,
263  nullValue,&status) || fits_flush_file(fPtr,&status) != 0)
264  {
265  throw FitsError(status);
266 
267  }
268 
269 
270 
271  m_image[std::slice(first-1,nElements,1)] = inData;
272  }
273 
274  template <typename T>
275  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)
276  {
277  // input vectors' size equality will be verified in prepareForSubset.
278  const size_t nDim = naxes.size();
279  FITSUtil::auto_array_ptr<long> pFPixel(new long[nDim]);
280  FITSUtil::auto_array_ptr<long> pLPixel(new long[nDim]);
281  std::valarray<T> subset;
282  prepareForSubset(naxes,firstVertex,lastVertex,stride,inData,subset);
283 
284  long* fPixel = pFPixel.get();
285  long* lPixel = pLPixel.get();
286  for (size_t i=0; i<nDim; ++i)
287  {
288  fPixel[i] = firstVertex[i];
289  lPixel[i] = lastVertex[i];
290  }
291 
292  FITSUtil::CAarray<T> convert;
293  FITSUtil::auto_array_ptr<T> pArray(convert(subset));
294  T* array = pArray.get();
295  FITSUtil::MatchType<T> imageType;
296  int status(0);
297 
298  if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status)
299  || fits_flush_file(fPtr,&status) != 0) throw FitsError(status);
300  }
301 
302  template <typename T>
303  std::valarray<T>& Image<T>::image ()
304  {
305 
306  return m_image;
307  }
308 
309  template <typename T>
310  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)
311  {
312 
313  // naxes, firstVertex, lastVertex, and stride must all be the same size.
314  const size_t N = naxes.size();
315  if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size())
316  {
317  string errMsg("*** CCfits Error: Image write function requires that naxes, firstVertex,");
318  errMsg += " \nlastVertex, and stride vectors all be the same size.\n";
319  bool silent = false;
320  throw FitsException(errMsg, silent);
321  }
322  for (size_t i=0; i<N; ++i)
323  {
324  if (naxes[i] < 1)
325  {
326  bool silent = false;
327  throw FitsException("*** CCfits Error: Invalid naxes value sent to image write function.\n", silent);
328  }
329  string rangeErrMsg("*** CCfits Error: Out-of-range value sent to image write function in arg: ");
330  if (firstVertex[i] < 1 || firstVertex[i] > naxes[i])
331  {
332  bool silent = false;
333  rangeErrMsg += "firstVertex\n";
334  throw FitsException(rangeErrMsg, silent);
335  }
336  if (lastVertex[i] < firstVertex[i] || lastVertex[i] > naxes[i])
337  {
338  bool silent = false;
339  rangeErrMsg += "lastVertex\n";
340  throw FitsException(rangeErrMsg, silent);
341  }
342  if (stride[i] < 1)
343  {
344  bool silent = false;
345  rangeErrMsg += "stride\n";
346  throw FitsException(rangeErrMsg, silent);
347  }
348  }
349 
350  // nPoints refers to the subset of m_image INCLUDING the zero'ed elements
351  // resulting from the stride parameter.
352  // subSizeWithStride refers to the same subset, not counting the zeros.
353  size_t subSizeWithStride = 1;
354  size_t nPoints = 1;
355  std::vector<size_t> subIncr(N);
356  for (size_t i=0; i<N; ++i)
357  {
358  subIncr[i] = nPoints;
359  nPoints *= static_cast<size_t>(1+lastVertex[i]-firstVertex[i]);
360  subSizeWithStride *= static_cast<size_t>(1+(lastVertex[i]-firstVertex[i])/stride[i]);
361  }
362  FITSUtil::FitsNullValue<T> null;
363  subset.resize(nPoints, null());
364 
365  // Trying to avoid at all costs an assignment between 2 valarrays of
366  // different sizes when m_image gets set below.
367  if (subSizeWithStride != inData.size())
368  {
369  bool silent = false;
370  string errMsg("*** CCfits Error: Data array size is not consistent with the values");
371  errMsg += "\n in range and stride vectors sent to the image write function.\n";
372  throw FitsException(errMsg, silent);
373  }
374 
375  size_t startPoint = 0;
376  size_t dimMult = 1;
377  std::vector<size_t> incr(N);
378  for (size_t j = 0; j < N; ++j)
379  {
380  startPoint += dimMult*(firstVertex[j]-1);
381  incr[j] = dimMult;
382  dimMult *= static_cast<size_t>(naxes[j]);
383  }
384  const size_t imageSize = dimMult;
385  m_image.resize(imageSize,null());
386 
387  size_t inDataPos = 0;
388  size_t iSub = 0;
389  loop(N-1, firstVertex, lastVertex, stride, startPoint, incr, inData, inDataPos, subIncr, subset, iSub);
390  }
391 
392  template <typename T>
393  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)
394  {
395  size_t start = static_cast<size_t>(firstVertex[iDim]);
396  size_t stop = static_cast<size_t>(lastVertex[iDim]);
397  size_t skip = static_cast<size_t>(stride[iDim]);
398  if (iDim == 0)
399  {
400  size_t length = stop - start + 1;
401  for (size_t i=0; i<length; i+=skip)
402  {
403  m_image[i+iPos] = inData[iDat];
404  subset[i+iSub] = inData[iDat++];
405  }
406  }
407  else
408  {
409  size_t jump = incr[iDim]*skip;
410  size_t subJump = subIncr[iDim]*skip;
411  for (size_t i=start; i<=stop; i+=skip)
412  {
413  loop(iDim-1, firstVertex, lastVertex, stride, iPos, incr, inData, iDat, subIncr, subset, iSub);
414  iPos += jump;
415  iSub += subJump;
416  }
417  }
418  }
419 
420  template <typename T>
421  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)
422  {
423  std::vector<long> stride(firstVertex.size(), 1);
424  writeImage(fPtr, firstVertex, lastVertex, stride, inData, naxes);
425  }
426 
427  // Additional Declarations
428 
429 } // namespace CCfits
430 
431 
432 #endif