CCfits  2.6
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 #include <sstream>
21 
22 #ifdef _MSC_VER
23 #include "MSconfig.h" //form std::min
24 #endif
25 #include "CCfits.h"
26 #include "FitsError.h"
27 #include "FITSUtil.h"
28 
29 
30 namespace CCfits {
31 
32 
33 
34  template <typename T>
35  class Image
36  {
37 
38  public:
39  Image(const Image< T > &right);
40  Image (const std::valarray<T>& imageArray = std::valarray<T>());
41  ~Image();
42  Image< T > & operator=(const Image< T > &right);
43 
44  const std::valarray<T>& readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls);
45  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);
46  // If write operation causes an expansion of the image's outer-most dimension, newNaxisN will be set to the new value. Else it will be 0.
47  void writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN, T* nullValue = 0);
48  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, long& newNaxisN);
49  void writeImage (fitsfile* fPtr, const std::vector<long>& firstVertex, const std::vector<long>& lastVertex, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN);
50  bool isRead () const;
51  // This allows higher level classes to notify Image that a user-input
52  // scaling value has changed. Image can then decide how this
53  // should affect reading from cache vs. disk.
54  void scalingHasChanged();
55  // Give the user (via higher level classes) a way to explicitly set the m_isRead flag
56  // to false, thus providing a fail-safe override of reading from the cache.
57  void resetRead();
58  const std::valarray< T >& image () const;
59 
60  // Additional Public Declarations
61 
62  protected:
63  // Additional Protected Declarations
64 
65  private:
66  std::valarray<T>& image ();
67  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);
68  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);
69  bool isNullValChanged(T* newNull) const;
70  void setLastNullInfo(T* newNull);
71 
72  // Additional Private Declarations
73 
74  private: //## implementation
75  // Data Members for Class Attributes
76  // When m_isRead = true, assume m_fullImageCache contains the full image from the file.
77  bool m_isRead;
78 
79  // Information regarding the usage of null values for the
80  // most recent read operation.
81  bool m_usingNullVal;
82  T m_lastNullVal;
83 
84  // Data Members for Associations
85  std::valarray< T > m_fullImageCache;
86  std::valarray<T> m_currentRead;
87 
88  // Additional Implementation Declarations
89 
90  };
91 
92  // Parameterized Class CCfits::Image
93 
94  template <typename T>
95  inline bool Image<T>::isRead () const
96  {
97  return m_isRead;
98  }
99 
100  template <typename T>
101  inline const std::valarray< T >& Image<T>::image () const
102  {
103  return m_fullImageCache;
104  }
105 
106 
107  // Parameterized Class CCfits::Image
108 
109  template <typename T>
110  Image<T>::Image(const Image<T> &right)
111  : m_isRead(right.m_isRead),
112  m_usingNullVal(right.m_usingNullVal),
113  m_lastNullVal(right.m_lastNullVal),
114  m_fullImageCache(right.m_fullImageCache),
115  m_currentRead(right.m_currentRead)
116  {
117  }
118 
119  template <typename T>
120  Image<T>::Image (const std::valarray<T>& imageArray)
121  : m_isRead(false),
122  m_usingNullVal(false),
123  m_lastNullVal(0),
124  m_fullImageCache(imageArray),
125  m_currentRead()
126  {
127  }
128 
129 
130  template <typename T>
131  Image<T>::~Image()
132  {
133  }
134 
135 
136  template <typename T>
137  Image<T> & Image<T>::operator=(const Image<T> &right)
138  {
139  // all stack allocated.
140  m_isRead = right.m_isRead;
141  m_usingNullVal = right.m_usingNullVal,
142  m_lastNullVal = right.m_lastNullVal,
143  m_fullImageCache.resize(right.m_fullImageCache.size());
144  m_fullImageCache = right.m_fullImageCache;
145  m_currentRead.resize(right.m_currentRead.size());
146  m_currentRead = right.m_currentRead;
147  return *this;
148  }
149 
150 
151  template <typename T>
152  const std::valarray<T>& Image<T>::readImage (fitsfile* fPtr, long first, long nElements, T* nullValue, const std::vector<long>& naxes, bool& nulls)
153  {
154  if (!naxes.size())
155  {
156  m_currentRead.resize(0);
157  return m_currentRead;
158  }
159  unsigned long init(1);
160  unsigned long nTotalElements(std::accumulate(naxes.begin(),naxes.end(),init,
161  std::multiplies<long>()));
162 
163  if (first <= 0)
164  {
165  string errMsg("*** CCfits Error: For image read, lowest allowed value for first element is 1.\n");
166  bool silent = false;
167  throw FitsException(errMsg, silent);
168  }
169  // 0-based index for slice
170  unsigned long start = (unsigned long)first - 1;
171  if (start >= nTotalElements)
172  {
173  string errMsg("*** CCfits Error: For image read, starting element is out of range.\n");
174  bool silent = false;
175  throw FitsException(errMsg, silent);
176  }
177  if (nElements < 0)
178  {
179  string errMsg("*** CCfits Error: Negative nElements value specified for image read.\n");
180  bool silent = false;
181  throw FitsException(errMsg, silent);
182  }
183  const unsigned long elementsRequested = (unsigned long)nElements;
184 
185  int status(0);
186  int any (0);
187  FITSUtil::MatchType<T> imageType;
188 
189  // truncate to valid array size if too much data asked for.
190  unsigned long elementsToRead(std::min(elementsRequested,
191  nTotalElements - start));
192  if ( elementsToRead < elementsRequested)
193  {
194  std::cerr <<
195  "***CCfits Warning: data request exceeds image size, truncating\n";
196  }
197  const bool isFullRead = (elementsToRead == nTotalElements);
198  const bool isDifferentNull = isNullValChanged(nullValue);
199  if (!m_isRead || isDifferentNull)
200  {
201  // Must perform a read from disk.
202  m_isRead = false;
203  if (isFullRead)
204  {
205  m_fullImageCache.resize(elementsToRead);
206  if (fits_read_img(fPtr,imageType(),first,elementsToRead,
207  nullValue,&m_fullImageCache[0],&any,&status) != 0) throw FitsError(status);
208  m_isRead = true;
209  // For this case only, we'll pass m_fullImageCache back up (to be
210  // copied into user-supplied array). This spares having to do
211  // what may be a very large copy into m_currentRead.
212  }
213  else
214  {
215  m_fullImageCache.resize(0);
216  m_currentRead.resize(elementsToRead);
217  if (fits_read_img(fPtr,imageType(),first,elementsToRead,
218  nullValue,&m_currentRead[0],&any,&status) != 0) throw FitsError(status);
219  }
220  nulls = (any != 0);
221  setLastNullInfo(nullValue);
222  }
223  else
224  {
225  if (!isFullRead)
226  {
227  m_currentRead.resize((size_t)elementsToRead);
228  // This may be a costly copy, though should still be faster
229  // than disk read.
230  m_currentRead = m_fullImageCache[std::slice((size_t)start, (size_t)elementsToRead,1)];
231  }
232  }
233  if (isFullRead)
234  return m_fullImageCache;
235 
236  return m_currentRead;
237 
238  }
239 
240  template <typename T>
241  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)
242  {
243  const size_t N = naxes.size();
244  if (!N)
245  {
246  m_currentRead.resize(0);
247  return m_currentRead;
248  }
249  if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size())
250  {
251  string errMsg("*** CCfits Error: Image read function requires that naxes, firstVertex,");
252  errMsg += " \nlastVertex, and stride vectors all be the same size.\n";
253  bool silent = false;
254  throw FitsException(errMsg, silent);
255  }
256 
257  FITSUtil::CVarray<long> carray;
258  int any(0);
259  int status(0);
260  long requestedSize=1;
261  long nTotalSize=1;
262 
263  for (size_t j = 0; j < N; ++j)
264  {
265  // Intentional truncation during division.
266  requestedSize *= ((lastVertex[j] - firstVertex[j])/stride[j] + 1);
267  nTotalSize *= naxes[j];
268  if (firstVertex[j] < 1 || lastVertex[j] > naxes[j])
269  {
270  string errMsg("*** CCfits Error: Out-of-bounds vertex value.\n");
271  bool silent=false;
272  throw FitsException(errMsg,silent);
273  }
274  if (firstVertex[j] > lastVertex[j])
275  {
276  string errMsg("*** CCfits Error: firstVertex values must not be larger than corresponding lastVertex values.\n");
277  bool silent = false;
278  throw FitsException(errMsg,silent);
279  }
280  }
281  const bool isFullRead = (requestedSize == nTotalSize);
282  const bool isDifferentNull = isNullValChanged(nullValue);
283  if (!m_isRead || isDifferentNull)
284  {
285  // Must perform a read from disk.
286  FITSUtil::auto_array_ptr<long> pFpixel(carray(firstVertex));
287  FITSUtil::auto_array_ptr<long> pLpixel(carray(lastVertex));
288  FITSUtil::auto_array_ptr<long> pStride(carray(stride));
289 
290  FITSUtil::MatchType<T> imageType;
291  m_isRead = false;
292  if (isFullRead)
293  {
294  m_fullImageCache.resize(requestedSize);
295  if (fits_read_subset(fPtr,imageType(),
296  pFpixel.get(),pLpixel.get(),
297  pStride.get(),nullValue,&m_fullImageCache[0],&any,&status) != 0)
298  throw FitsError(status);
299  m_isRead = true;
300  }
301  else
302  {
303  m_currentRead.resize(requestedSize);
304  if (fits_read_subset(fPtr,imageType(),
305  pFpixel.get(),pLpixel.get(),
306  pStride.get(),nullValue,&m_currentRead[0],&any,&status) != 0)
307  throw FitsError(status);
308  }
309  nulls = (any != 0);
310  setLastNullInfo(nullValue);
311  }
312  else
313  {
314  if (!isFullRead)
315  {
316  // Must convert firstVertex,lastVertex,stride to gslice parameters.
317  // Note that in cfitsio, the NAXIS1 dimension varies the fastest
318  // when laid out in an array in memory (ie. Fortran style). Therefore NAXISn
319  // ordering must be reversed to C style before passing to gslice.
320  size_t startPos=0;
321  std::valarray<size_t> gsliceLength(size_t(0),N);
322  std::valarray<size_t> gsliceStride(size_t(0),N);
323 
324  std::vector<long> naxesProducts(N);
325  long accum=1;
326  for (size_t i=0; i<N; ++i)
327  {
328  naxesProducts[i] = accum;
329  accum *= naxes[i];
330  }
331 
332  for (size_t i=0; i<N; ++i)
333  {
334  startPos += static_cast<size_t>((firstVertex[i]-1)*naxesProducts[i]);
335  // Here's where we reverse the order:
336  const size_t gsPos = N-1-i;
337  // Division truncation is intentional.
338  gsliceLength[gsPos] = static_cast<size_t>((1 + (lastVertex[i]-firstVertex[i])/stride[i]));
339  gsliceStride[gsPos] = static_cast<size_t>(stride[i]*naxesProducts[i]);
340  }
341  m_currentRead.resize(requestedSize);
342  m_currentRead = m_fullImageCache[std::gslice(startPos, gsliceLength, gsliceStride)];
343  }
344  }
345  if (isFullRead)
346  return m_fullImageCache;
347  return m_currentRead;
348  }
349 
350  template <typename T>
351  void Image<T>::writeImage (fitsfile* fPtr, long first, long nElements, const std::valarray<T>& inData, const std::vector<long>& naxes, long& newNaxisN, T* nullValue)
352  {
353  int status(0);
354  if (first < 1 || nElements < 1)
355  {
356  string errMsg("*** CCfits Error: first and nElements values must be > 0\n");
357  bool silent = false;
358  throw FitsException(errMsg, silent);
359  }
360  FITSUtil::CAarray<T> convert;
361  FITSUtil::auto_array_ptr<T> pArray(convert(inData));
362  T* array = pArray.get();
363 
364  m_isRead = false;
365  newNaxisN = 0;
366 
367  FITSUtil::MatchType<T> imageType;
368  long type(imageType());
369 
370  if (fits_write_imgnull(fPtr,type,first,nElements,array,
371  nullValue,&status)!= 0)
372  {
373  throw FitsError(status);
374  }
375  const size_t nDim=naxes.size();
376  long origTotSize=1;
377  for (size_t i=0; i<nDim; ++i)
378  origTotSize *= naxes[i];
379  const long highestOutputElem = first + nElements - 1;
380  if (highestOutputElem > origTotSize)
381  {
382  // NAXIS(nDIM) may have increased.
383  std::ostringstream oss;
384  oss <<"NAXIS" << nDim;
385  string keyname(oss.str());
386  long newVal = 1 + (highestOutputElem-1)/(origTotSize/naxes[nDim-1]);
387  if (newVal != naxes[nDim-1])
388  {
389  if (fits_update_key(fPtr,TLONG,(char *)keyname.c_str(),&newVal,0,&status) != 0)
390  {
391  throw FitsError(status);
392  }
393  newNaxisN = newVal;
394  }
395  }
396  if (fits_flush_file(fPtr,&status) != 0)
397  throw FitsError(status);
398 
399  }
400 
401  template <typename T>
402  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, long& newNaxisN)
403  {
404  // input vectors' size equality will be verified in prepareForSubset.
405  const size_t nDim = naxes.size();
406  FITSUtil::auto_array_ptr<long> pFPixel(new long[nDim]);
407  FITSUtil::auto_array_ptr<long> pLPixel(new long[nDim]);
408  std::valarray<T> subset;
409  m_isRead = false;
410  newNaxisN = 0;
411  prepareForSubset(naxes,firstVertex,lastVertex,stride,inData,subset);
412 
413  long* fPixel = pFPixel.get();
414  long* lPixel = pLPixel.get();
415  for (size_t i=0; i<nDim; ++i)
416  {
417  fPixel[i] = firstVertex[i];
418  lPixel[i] = lastVertex[i];
419  }
420 
421  FITSUtil::CAarray<T> convert;
422  FITSUtil::auto_array_ptr<T> pArray(convert(subset));
423  T* array = pArray.get();
424  FITSUtil::MatchType<T> imageType;
425  int status(0);
426 
427  if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status) )
428  throw FitsError(status);
429 
430  if (lPixel[nDim-1] > naxes[nDim-1])
431  {
432  std::ostringstream oss;
433  oss << "NAXIS" << nDim;
434  string keyname(oss.str());
435  long newVal = lPixel[nDim-1];
436  if (fits_update_key(fPtr,TLONG,(char *)keyname.c_str(),&newVal,0,&status) != 0)
437  {
438  throw FitsError(status);
439  }
440  newNaxisN = lPixel[nDim-1];
441  }
442  if (fits_flush_file(fPtr,&status) != 0)
443  throw FitsError(status);
444 
445  }
446 
447  template <typename T>
448  std::valarray<T>& Image<T>::image ()
449  {
450 
451  return m_fullImageCache;
452  }
453 
454  template <typename T>
455  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)
456  {
457 
458  // naxes, firstVertex, lastVertex, and stride must all be the same size.
459  const size_t N = naxes.size();
460  if (N != firstVertex.size() || N != lastVertex.size() || N != stride.size())
461  {
462  string errMsg("*** CCfits Error: Image write function requires that naxes, firstVertex,");
463  errMsg += " \nlastVertex, and stride vectors all be the same size.\n";
464  bool silent = false;
465  throw FitsException(errMsg, silent);
466  }
467  for (size_t i=0; i<N; ++i)
468  {
469  if (naxes[i] < 1)
470  {
471  bool silent = false;
472  throw FitsException("*** CCfits Error: Invalid naxes value sent to image write function.\n", silent);
473  }
474  string rangeErrMsg("*** CCfits Error: Out-of-range value sent to image write function in arg: ");
475  if (firstVertex[i] < 1 || (firstVertex[i] > naxes[i] && i != N-1))
476  {
477  bool silent = false;
478  rangeErrMsg += "firstVertex\n";
479  throw FitsException(rangeErrMsg, silent);
480  }
481  if (lastVertex[i] < firstVertex[i] || (lastVertex[i] > naxes[i] && i != N-1))
482  {
483  bool silent = false;
484  rangeErrMsg += "lastVertex\n";
485  throw FitsException(rangeErrMsg, silent);
486  }
487  if (stride[i] < 1)
488  {
489  bool silent = false;
490  rangeErrMsg += "stride\n";
491  throw FitsException(rangeErrMsg, silent);
492  }
493  }
494 
495  // nPoints refers to the subset of the image INCLUDING the zero'ed elements
496  // resulting from the stride parameter.
497  // subSizeWithStride refers to the same subset, not counting the zeros.
498  size_t subSizeWithStride = 1;
499  size_t nPoints = 1;
500  std::vector<size_t> subIncr(N);
501  for (size_t i=0; i<N; ++i)
502  {
503  subIncr[i] = nPoints;
504  nPoints *= static_cast<size_t>(1+lastVertex[i]-firstVertex[i]);
505  subSizeWithStride *= static_cast<size_t>(1+(lastVertex[i]-firstVertex[i])/stride[i]);
506  }
507  subset.resize(nPoints, 0);
508 
509  if (subSizeWithStride != inData.size())
510  {
511  bool silent = false;
512  string errMsg("*** CCfits Error: Data array size is not consistent with the values");
513  errMsg += "\n in range and stride vectors sent to the image write function.\n";
514  throw FitsException(errMsg, silent);
515  }
516 
517  size_t startPoint = 0;
518  size_t dimMult = 1;
519  std::vector<size_t> incr(N);
520  for (size_t j = 0; j < N; ++j)
521  {
522  startPoint += dimMult*(firstVertex[j]-1);
523  incr[j] = dimMult;
524  dimMult *= static_cast<size_t>(naxes[j]);
525  }
526 
527  size_t inDataPos = 0;
528  size_t iSub = 0;
529  loop(N-1, firstVertex, lastVertex, stride, startPoint, incr, inData, inDataPos, subIncr, subset, iSub);
530  }
531 
532  template <typename T>
533  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)
534  {
535  size_t start = static_cast<size_t>(firstVertex[iDim]);
536  size_t stop = static_cast<size_t>(lastVertex[iDim]);
537  size_t skip = static_cast<size_t>(stride[iDim]);
538  if (iDim == 0)
539  {
540  size_t length = stop - start + 1;
541  for (size_t i=0; i<length; i+=skip)
542  {
543  subset[i+iSub] = inData[iDat++];
544  }
545  }
546  else
547  {
548  size_t jump = incr[iDim]*skip;
549  size_t subJump = subIncr[iDim]*skip;
550  for (size_t i=start; i<=stop; i+=skip)
551  {
552  loop(iDim-1, firstVertex, lastVertex, stride, iPos, incr, inData, iDat, subIncr, subset, iSub);
553  iPos += jump;
554  iSub += subJump;
555  }
556  }
557  }
558 
559  template <typename T>
560  bool Image<T>::isNullValChanged(T* newNull) const
561  {
562  bool isChanged = false;
563  if (m_usingNullVal)
564  {
565  // If m_usingNullVal is true, we can assume m_lastNullVal != 0.
566  if (newNull)
567  {
568  T newVal = *newNull;
569  if (newVal != m_lastNullVal)
570  isChanged = true;
571  }
572  else
573  isChanged = true;
574  }
575  else
576  {
577  if (newNull && (*newNull != 0))
578  isChanged = true;
579  }
580 
581  return isChanged;
582  }
583 
584  template <typename T>
585  void Image<T>::setLastNullInfo(T* newNull)
586  {
587  if (!newNull || *newNull == 0)
588  {
589  m_usingNullVal = false;
590  m_lastNullVal = 0;
591  }
592  else
593  {
594  m_usingNullVal = true;
595  m_lastNullVal = *newNull;
596  }
597  }
598 
599  template <typename T>
600  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, long& newNaxisN)
601  {
602  std::vector<long> stride(firstVertex.size(), 1);
603  writeImage(fPtr, firstVertex, lastVertex, stride, inData, naxes, newNaxisN);
604  }
605 
606  template <typename T>
607  void Image<T>::scalingHasChanged()
608  {
609  m_isRead = false;
610  }
611 
612  template <typename T>
613  void Image<T>::resetRead()
614  {
615  m_isRead = false;
616  }
617 
618  // Additional Declarations
619 
620 } // namespace CCfits
621 
622 
623 #endif
Namespace enclosing all CCfits classes and globals definitions.
Definition: AsciiTable.cxx:26