相关文章推荐
小百科  ›  CCfits: Image.h Source File
const string vector
小眼睛的脆皮肠
1 年前
1 // Astrophysics Science Division,
2 // NASA/ Goddard Space Flight Center
3 // HEASARC
4 // http://heasarc.gsfc.nasa.gov
5 // e-mail: [email protected]
6 //
7 // Original author: Ben Dorman
9 #ifndef IMAGE_H
10 #define IMAGE_H 1
12 // functional
13 #include <functional>
14 // valarray
15 #include <valarray>
16 // vector
17 #include <vector>
18 // numeric
19 #include <numeric>
20 #include <sstream>
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"
30 namespace CCfits {
34 template < typename T>
35 class Image
36 {
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);
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 ;
60 // Additional Public Declarations
62 protected :
63 // Additional Protected Declarations
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);
72 // Additional Private Declarations
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;
79 // Information regarding the usage of null values for the
80 // most recent read operation.
81 bool m_usingNullVal;
82 T m_lastNullVal;
84 // Data Members for Associations
85 std::valarray< T > m_fullImageCache;
86 std::valarray<T> m_currentRead;
88 // Additional Implementation Declarations
90 };
92 // Parameterized Class CCfits::Image
94 template < typename T>
95 inline bool Image<T>::isRead () const
96 {
97 return m_isRead;
98 }
100 template < typename T>
101 inline const std::valarray< T >& Image<T>::image () const
102 {
103 return m_fullImageCache;
104 }
107 // Parameterized Class CCfits::Image
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 }
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 }
130 template < typename T>
131 Image<T>::~Image()
132 {
133 }
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 }
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>()));
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;
185 int status(0);
186 int any (0);
187 FITSUtil::MatchType<T> imageType;
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;
236 return m_currentRead;
238 }
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 }
257 FITSUtil::CVarray<long> carray;
258 int any(0);
259 int status(0);
260 long requestedSize=1;
261 long nTotalSize=1;
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));
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);
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 }
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 }
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();
364 m_isRead = false ;
365 newNaxisN = 0;
367 FITSUtil::MatchType<T> imageType;
368 long type(imageType());
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);
399 }
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);
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 }
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);
427 if ( fits_write_subset(fPtr,imageType(),fPixel,lPixel,array,&status) )
428 throw FitsError(status);
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);
445 }
447 template < typename T>
448 std::valarray<T>& Image<T>::image ()
449 {
451 return m_fullImageCache;
452 }
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 {
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 }
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);
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 }
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 }
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 }
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 }
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 }
581 return isChanged;
582 }
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 }
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 }
606 template < typename T>
607 void Image<T>::scalingHasChanged()
608 {
609 m_isRead = false ;
610 }
612 template < typename T>
613 void Image<T>::resetRead()
614 {
615 m_isRead = false ;
616 }
618 // Additional Declarations
620 } // namespace CCfits
623 #endif
 
推荐文章
Link管理   ·   51好读   ·   Sov5搜索   ·   小百科
小百科 - 百科知识指南