LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
24
25namespace 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}
67namespace {
68
69void 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) {
89 return lsst::geom::Box2I(lsst::geom::Point2I(), _bbox.getDimensions());
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
143namespace {
144
145// Return a numpy-like description of a primitive type (e.g. uint8, float64).
146template <typename T>
147std::string makeDTypeString() {
148 using L = std::numeric_limits<T>;
149 static_assert(L::is_specialized, "makeDTypeString requires a primitive numeric template parameter.");
150 std::string result;
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
170template <typename T>
171ndarray::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
211INSTANTIATE(int);
212INSTANTIATE(float);
213INSTANTIATE(double);
215
216}}} // lsst::afw::image
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
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:308
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:782
Lifetime-management for memory that goes into FITS memory files.
Definition fits.h:125
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
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.
std::shared_ptr< daf::base::PropertyList > readMetadata(std::string const &fileName, int hdu=DEFAULT_HDU, bool strip=false)
Read FITS header.
Definition fits.cc:1696
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
Extent< int, 2 > Extent2I
Definition Extent.h:397
Point< int, 2 > Point2I
Definition Point.h:321
py::scoped_interpreter guard
Definition test_image.cc:19
T to_string(T... args)