LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
ImageBaseFitsReader.cc
Go to the documentation of this file.
1 /*
2  * Developed for the LSST Data Management System.
3  * This product includes software developed by the LSST Project
4  * (https://www.lsst.org).
5  * See the COPYRIGHT file at the top-level directory of this distribution
6  * for details of code ownership.
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program. If not, see <https://www.gnu.org/licenses/>.
20  */
21 
23 #include "lsst/afw/geom/wcsUtils.h"
24 
25 namespace lsst { namespace afw { namespace image {
26 
28  _ownsFitsFile(true),
29  _hdu(0),
30  _fitsFile(new fits::Fits(fileName, "r", fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK))
31 {
32  _fitsFile->setHdu(hdu);
33  _fitsFile->checkCompressedImagePhu();
34  _hdu = _fitsFile->getHdu();
35 }
36 
38  _ownsFitsFile(true),
39  _hdu(0),
40  _fitsFile(new fits::Fits(manager, "r", fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK))
41 {
42  _fitsFile->setHdu(hdu);
43  _fitsFile->checkCompressedImagePhu();
44  _hdu = _fitsFile->getHdu();
45 }
46 
48  _ownsFitsFile(false),
49  _hdu(0),
50  _fitsFile(fitsFile)
51 {
52  if (_fitsFile) {
53  if (_fitsFile->getHdu() == 0 && _fitsFile->getImageDim() == 0) {
54  _fitsFile->setHdu(fits::DEFAULT_HDU);
55  }
56  _fitsFile->checkCompressedImagePhu();
57  _hdu = _fitsFile->getHdu();
58  }
59 }
60 
62  if (_ownsFitsFile) {
63  delete _fitsFile;
64  }
65 }
66 
67 namespace {
68 
69 void checkFitsFile(fits::Fits * f) {
70  if (f == nullptr) {
71  throw LSST_EXCEPT(
73  "FitsReader not initialized; desired HDU is probably missing."
74  );
75  }
76 }
77 
78 } // anonymous
79 
81  checkFitsFile(_fitsFile);
82  fits::HduMoveGuard guard(*_fitsFile, _hdu);
83  return _fitsFile->getImageDType();
84 }
85 
87  readMetadata(); // guarantees _bbox is initialized
88  if (origin == LOCAL) {
90  }
91  return _bbox;
92 }
93 
95  if (bbox.isEmpty()) {
96  return readBBox().getMin();
97  } else if (origin == PARENT) {
98  return bbox.getMin();
99  } else {
100  auto full = readBBox();
101  return full.getMin() + lsst::geom::Extent2I(bbox.getMin());
102  }
103 }
104 
106  checkFitsFile(_fitsFile);
107  if (_metadata == nullptr) {
108  fits::HduMoveGuard guard(*_fitsFile, _hdu);
109  auto metadata = fits::readMetadata(*_fitsFile, /*strip=*/true);
110  // One-use lambda avoids declare-then-initialize; see C++ Core Guidelines, ES28.
111  auto computeShape = [this](){
112  int nAxis = _fitsFile->getImageDim();
113  if (nAxis == 2) {
114  return _fitsFile->getImageShape<2>();
115  }
116  if (nAxis == 3) {
117  ndarray::Vector<ndarray::Size, 3> shape3 = _fitsFile->getImageShape<3>();
118  if (shape3[0] != 1) {
119  throw LSST_EXCEPT(
121  str(boost::format("Error reading %s: HDU %d has 3rd dimension %d != 1") %
122  getFileName() % _hdu % shape3[0]));
123  }
124  return shape3.last<2>();
125  }
126  throw LSST_EXCEPT(
128  str(boost::format("Error reading %s: HDU %d has %d dimensions") %
129  getFileName() % _hdu % nAxis)
130  );
131  };
132  // Construct _bbox and check dimensions at the same time both to make
133  // one less valid state for this class to be in and to make sure
134  // the appropriate metadata is stripped before we ever return it.
135  ndarray::Vector<ndarray::Size, 2> shape = computeShape();
136  auto xy0 = afw::geom::getImageXY0FromMetadata(*metadata, detail::wcsNameForXY0, /*strip=*/true);
137  _bbox = lsst::geom::Box2I(xy0, lsst::geom::Extent2I(shape[1], shape[0]));
138  _metadata = std::move(metadata);
139  }
140  return _metadata;
141 }
142 
143 namespace {
144 
145 // Return a numpy-like description of a primitive type (e.g. uint8, float64).
146 template <typename T>
147 std::string makeDTypeString() {
148  using L = std::numeric_limits<T>;
149  static_assert(L::is_specialized, "makeDTypeString requires a primitive numeric template parameter.");
152  result += "bool";
153  } else {
154  if (L::is_integer) {
155  if (L::is_signed) {
156  result += "int";
157  } else {
158  result += "uint";
159  }
160  } else {
161  result += "float";
162  }
163  result += std::to_string(sizeof(T)*8);
164  }
165  return result;
166 }
167 
168 } // anonymous
169 
170 template <typename T>
171 ndarray::Array<T, 2, 2> ImageBaseFitsReader::readArray(lsst::geom::Box2I const & bbox, ImageOrigin origin,
172  bool allowUnsafe) {
173  checkFitsFile(_fitsFile);
174  auto fullBBox = readBBox(origin);
175  auto subBBox = bbox;
176  if (subBBox.isEmpty()) {
177  subBBox = fullBBox;
178  } else if (!fullBBox.contains(subBBox)) {
179  throw LSST_EXCEPT(
181  str(boost::format("Subimage box (%d,%d) %dx%d doesn't fit in image (%d,%d) %dx%d in HDU %d") %
182  subBBox.getMinX() % subBBox.getMinY() % subBBox.getWidth() % subBBox.getHeight() %
183  fullBBox.getMinX() % fullBBox.getMinY() % fullBBox.getWidth() % fullBBox.getHeight() % _hdu)
184  );
185  }
186  fits::HduMoveGuard guard(*_fitsFile, _hdu);
187  if (!allowUnsafe & !_fitsFile->checkImageType<T>()) {
188  throw LSST_FITS_EXCEPT(
189  fits::FitsTypeError, *_fitsFile,
190  str(boost::format("Incompatible type for FITS image: on disk is %s (HDU %d), in-memory is %s. "
191  "Read with allowUnsafe=True to permit conversions that may overflow.") %
192  _fitsFile->getImageDType() % _hdu % makeDTypeString<T>())
193  );
194  }
195  ndarray::Array<T, 2, 2> result = ndarray::allocate(subBBox.getHeight(), subBBox.getWidth());
196  ndarray::Vector<int, 2> offset = ndarray::makeVector(subBBox.getMinY() - fullBBox.getMinY(),
197  subBBox.getMinX() - fullBBox.getMinX());
198  _fitsFile->readImage(result, offset);
199  return result;
200 }
201 
202 
203 #define INSTANTIATE(T) \
204  template ndarray::Array<T, 2, 2> ImageBaseFitsReader::readArray( \
205  lsst::geom::Box2I const & bbox, \
206  ImageOrigin origin, \
207  bool \
208  )
209 
211 INSTANTIATE(int);
212 INSTANTIATE(float);
213 INSTANTIATE(double);
215 
216 }}} // lsst::afw::image
py::object result
Definition: _schema.cc:429
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Fits * fits
Definition: FitsWriter.cc:90
#define INSTANTIATE(T)
An exception thrown when problems are found when reading or writing FITS files.
Definition: fits.h:36
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Definition: fits.h:297
void readImage(ndarray::Array< T, N, N > const &array, ndarray::Vector< int, N > const &offset)
Read an array from a FITS image.
Definition: fits.h:541
int getImageDim()
Return the number of dimensions in the current HDU.
Definition: fits.cc:1422
void setHdu(int hdu, bool relative=false)
Set the current HDU.
Definition: fits.cc:513
bool checkCompressedImagePhu()
Go to the first image header in the FITS file.
Definition: fits.cc:1726
ndarray::Vector< ndarray::Size, N > getImageShape()
Return the shape of the current (image) HDU.
Definition: fits.h:512
std::string getImageDType()
Return the numpy dtype equivalent of the image pixel type (e.g.
Definition: fits.cc:1466
int getHdu()
Return the current HDU (0-indexed; 0 is the Primary HDU).
Definition: fits.cc:507
bool checkImageType()
Return true if the current HDU is compatible with the given pixel type.
Definition: fits.cc:1435
An exception thrown when a FITS file has the wrong type.
Definition: fits.h:41
RAII scoped guard for moving the HDU in a Fits object.
Definition: fits.h:724
Lifetime-management for memory that goes into FITS memory files.
Definition: fits.h:121
lsst::geom::Point2I readXY0(lsst::geom::Box2I const &bbox=lsst::geom::Box2I(), ImageOrigin origin=PARENT)
Read the image origin from the on-disk image or a subimage thereof.
std::shared_ptr< daf::base::PropertyList > readMetadata()
Read the image's FITS header.
ndarray::Array< T, 2, 2 > readArray(lsst::geom::Box2I const &bbox, ImageOrigin origin=PARENT, bool allowUnsafe=false)
Read the image's data array.
std::string getFileName() const
Return the name of the file this reader targets.
lsst::geom::Box2I readBBox(ImageOrigin origin=PARENT)
Read the bounding box of the on-disk image.
std::string readDType() const
Read a string describing the pixel type of the on-disk image.
ImageBaseFitsReader(std::string const &fileName, int hdu=fits::DEFAULT_HDU)
Construct a FITS reader object.
An integer coordinate rectangle.
Definition: Box.h:55
Point2I const getMin() const noexcept
Definition: Box.h:156
Extent2I const getDimensions() const noexcept
Definition: Box.h:186
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
#define LSST_FITS_EXCEPT(type, fitsObj,...)
A FITS-related replacement for LSST_EXCEPT that takes an additional Fits object and uses makeErrorMes...
Definition: fits.h:105
T move(T... args)
const int DEFAULT_HDU
Specify that the default HDU should be read.
Definition: fitsDefaults.h:18
std::shared_ptr< daf::base::PropertyList > readMetadata(std::string const &fileName, int hdu=DEFAULT_HDU, bool strip=false)
Read FITS header.
Definition: fits.cc:1657
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:95
std::string const wcsNameForXY0
Definition: ImageBase.h:70
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Extent< int, 2 > Extent2I
Definition: Extent.h:397
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T to_string(T... args)