LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
fitsCompression.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2
3#include "fitsio.h"
4extern "C" {
5#include "fitsio2.h"
6}
7
9
11
12extern float* fits_rand_value; // Random numbers, defined in cfitsio
13int const N_RESERVED_VALUES = 10; // Number of reserved values for float --> bitpix=32 conversions (cfitsio)
14
15namespace lsst {
16namespace afw {
17namespace fits {
18
20 if (name == "NONE") return ImageCompressionOptions::NONE;
21 if (name == "GZIP") return ImageCompressionOptions::GZIP;
22 if (name == "GZIP_SHUFFLE") return ImageCompressionOptions::GZIP_SHUFFLE;
23 if (name == "RICE") return ImageCompressionOptions::RICE;
24 if (name == "HCOMPRESS")
25 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "HCOMPRESS is unsupported");
26 if (name == "PLIO") return ImageCompressionOptions::PLIO;
27 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Unrecognised compression algorithm: " + name);
28}
29
31 switch (algorithm) {
33 return "NONE";
35 return "GZIP";
37 return "GZIP_SHUFFLE";
39 return "RICE";
41 return "PLIO";
42 default:
44 os << "Unrecognized compression algorithm: " << algorithm;
46 }
47}
48
50 switch (cfitsio) {
51 case 0:
53 case RICE_1:
55 case GZIP_1:
57 case GZIP_2:
59 case PLIO_1:
61 case HCOMPRESS_1:
63 "Unsupported compression algorithm: HCOMPRESS_1");
64 default:
66 os << "Unrecognized cfitsio compression: " << cfitsio;
68 }
69}
70
72 switch (algorithm) {
74 return 0;
76 return GZIP_1;
78 return GZIP_2;
80 return RICE_1;
82 return PLIO_1;
83 default:
85 os << "Unrecognized compression algorithm: " << algorithm;
87 }
88}
89
91 int rows, float quantizeLevel_)
92 : algorithm(algorithm_), tiles(ndarray::allocate(MAX_COMPRESS_DIM)), quantizeLevel(quantizeLevel_) {
93 tiles[0] = 0;
94 tiles[1] = rows;
95 for (int ii = 2; ii < MAX_COMPRESS_DIM; ++ii) tiles[ii] = 1;
96}
97
99 if (name == "NONE") return ImageScalingOptions::NONE;
100 if (name == "RANGE") return ImageScalingOptions::RANGE;
101 if (name == "STDEV_POSITIVE") return ImageScalingOptions::STDEV_POSITIVE;
102 if (name == "STDEV_NEGATIVE") return ImageScalingOptions::STDEV_NEGATIVE;
103 if (name == "STDEV_BOTH") return ImageScalingOptions::STDEV_BOTH;
104 if (name == "MANUAL") return ImageScalingOptions::MANUAL;
105 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "Unrecognized scaling algorithm: " + name);
106}
107
109 switch (algorithm) {
111 return "NONE";
113 return "RANGE";
115 return "STDEV_POSITIVE";
117 return "STDEV_NEGATIVE";
119 return "STDEV_BOTH";
121 return "MANUAL";
122 default:
124 os << "Unrecognized scaling algorithm: " << algorithm;
126 }
127}
128
130 std::vector<std::string> const& maskPlanes_, int seed_,
131 float quantizeLevel_, float quantizePad_, bool fuzz_, double bscale_,
132 double bzero_)
133 : algorithm(algorithm_),
134 bitpix(bitpix_),
135 fuzz(fuzz_),
136 seed(std::abs(seed_ - 1) % (N_RANDOM - 1) +
137 1), // zero is bad (cfitsio uses non-deterministic method)
138 maskPlanes(maskPlanes_),
139 quantizeLevel(quantizeLevel_),
140 quantizePad(quantizePad_),
141 bscale(bscale_),
142 bzero(bzero_) {}
143
144namespace {
145
147template <typename T, int N>
148std::pair<T, T> calculateMedianStdev(ndarray::Array<T const, N, N> const& image,
149 ndarray::Array<bool, N, N> const& mask) {
150 std::size_t num = 0;
151 auto const& flatMask = ndarray::flatten<1>(mask);
152 for (auto mm = flatMask.begin(); mm != flatMask.end(); ++mm) {
153 if (!*mm) ++num;
154 }
155 ndarray::Array<T, 1, 1> array = ndarray::allocate(num);
156 auto const& flatImage = ndarray::flatten<1>(image);
157 auto mm = ndarray::flatten<1>(mask).begin();
158 auto aa = array.begin();
159 for (auto ii = flatImage.begin(); ii != flatImage.end(); ++ii, ++mm) {
160 if (*mm) continue;
161 *aa = *ii;
162 ++aa;
163 }
164
165 // Quartiles; from https://stackoverflow.com/a/11965377/834250
166 auto const q1 = num / 4;
167 auto const q2 = num / 2;
168 auto const q3 = q1 + q2;
169 std::nth_element(array.begin(), array.begin() + q1, array.end());
170 std::nth_element(array.begin() + q1 + 1, array.begin() + q2, array.end());
171 std::nth_element(array.begin() + q2 + 1, array.begin() + q3, array.end());
172
173 T const median = num % 2 ? array[num / 2] : 0.5 * (array[num / 2] + array[num / 2 - 1]);
174 // No, we're not doing any interpolation for the lower and upper quartiles.
175 // We're estimating the noise, so it doesn't need to be super precise.
176 T const lq = array[q1];
177 T const uq = array[q3];
178 return std::make_pair(median, 0.741 * (uq - lq));
179}
180
182template <typename T, int N>
183std::pair<T, T> calculateMinMax(ndarray::Array<T const, N, N> const& image,
184 ndarray::Array<bool, N, N> const& mask) {
186 auto mm = ndarray::flatten<1>(mask).begin();
187 auto const& flatImage = ndarray::flatten<1>(image);
188 for (auto ii = flatImage.begin(); ii != flatImage.end(); ++ii, ++mm) {
189 if (*mm) continue;
190 if (!std::isfinite(*ii)) continue;
191 if (*ii > max) max = *ii;
192 if (*ii < min) min = *ii;
193 }
194 return std::make_pair(min, max);
195}
196
197// Return range of values for target BITPIX
198template <typename T>
199double rangeForBitpix(int bitpix, bool cfitsioPadding) {
200 if (bitpix == 0) {
201 bitpix = detail::Bitpix<T>::value;
202 }
203 double range = std::pow(2.0, bitpix) - 1; // Range of values for target BITPIX
204 if (cfitsioPadding) {
205 range -= N_RESERVED_VALUES;
206 }
207 return range;
208}
209
210} // anonymous namespace
211
212template <typename T, int N>
213ImageScale ImageScalingOptions::determineFromRange(ndarray::Array<T const, N, N> const& image,
214 ndarray::Array<bool, N, N> const& mask, bool isUnsigned,
215 bool cfitsioPadding) const {
216 auto minMax = calculateMinMax(image, mask);
217 T const min = minMax.first;
218 T const max = minMax.second;
219 if (min == max) return ImageScale(bitpix, 1.0, min);
220 double range = rangeForBitpix<T>(bitpix, cfitsioPadding);
221 range -= 2; // To allow for rounding and fuzz at either end
222 double const bscale = static_cast<T>((max - min) / range);
223 double bzero = static_cast<T>(isUnsigned ? min : min + 0.5 * range * bscale);
224 if (cfitsioPadding) {
226 }
227 bzero -= bscale; // Allow for rounding and fuzz on the low end
228 return ImageScale(bitpix, bscale, bzero);
229}
230
231template <typename T, int N>
232ImageScale ImageScalingOptions::determineFromStdev(ndarray::Array<T const, N, N> const& image,
233 ndarray::Array<bool, N, N> const& mask, bool isUnsigned,
234 bool cfitsioPadding) const {
235 auto stats = calculateMedianStdev(image, mask);
236 auto const median = stats.first, stdev = stats.second;
237 double const bscale = static_cast<T>(stdev / quantizeLevel);
238
240 auto minMax = calculateMinMax(image, mask);
241 T const min = minMax.first;
242 T const max = minMax.second;
243 double range = rangeForBitpix<T>(bitpix, cfitsioPadding); // Range of values for target BITPIX
244 double const numUnique = (max - min) / bscale; // Number of unique values
245
246 double imageVal; // Value on image
247 long diskVal; // Corresponding quantized value
248 if (numUnique < range) {
249 imageVal = median;
250 diskVal = cfitsioPadding ? 0.5 * N_RESERVED_VALUES : 0;
251 } else {
252 switch (algorithm) {
254 // Put (median - N sigma) at the lowest possible value: predominantly positive images
255 imageVal = median - quantizePad * stdev;
256 diskVal = (bitpix == 8) ? 0 : -(1L << (bitpix - 1)); // Lowest value: -2^(bitpix-1)
257 if (cfitsioPadding) diskVal -= N_RESERVED_VALUES;
258 break;
260 // Put (median + N sigma) at the highest possible value: predominantly negative images
261 imageVal = median + quantizePad * stdev;
262 diskVal = (bitpix == 8) ? 255 : (1L << (bitpix - 1)) - 1; // Lowest value: 2^(bitpix-1)-1
263 break;
265 // Put median right in the middle: images with an equal abundance of positive and negative
266 // values
267 imageVal = median;
268 diskVal = cfitsioPadding ? 0.5 * N_RESERVED_VALUES : 0;
269 break;
270 default:
271 std::abort(); // Programming error: should never get here
272 }
273 }
274
275 double bzero = static_cast<T>(imageVal - bscale * diskVal);
276 return ImageScale(bitpix, bscale, bzero);
277}
278
280template <typename T, class Enable = void>
281struct Bzero {
282 static double constexpr value = 0.0;
283};
284
285// uint64 version
286// 'double' doesn't have sufficient bits to represent the appropriate BZERO,
287// so let cfitsio handle it.
288template <>
289struct Bzero<std::uint64_t> {
290 static double constexpr value = 0.0;
291};
292
293// Unsigned integer version
294template <typename T>
295struct Bzero<T, typename std::enable_if<std::numeric_limits<T>::is_integer &&
296 !std::numeric_limits<T>::is_signed>::type> {
297 static double constexpr value = (std::numeric_limits<T>::max() >> 1) + 1;
298};
299
300#ifndef DOXYGEN // suppress a bogus Doxygen complaint about an documented symbol
301template <typename T, int N>
302ImageScale ImageScalingOptions::determine(ndarray::Array<T const, N, N> const& image,
303 ndarray::Array<bool, N, N> const& mask) const {
305 algorithm != NONE) {
307 "Image scaling not supported for integral types");
308 }
309 bool const isUnsigned = bitpix == 8 || (bitpix == 0 && detail::Bitpix<T>::value == 8);
310 bool const cfitsioPadding = !std::numeric_limits<T>::is_integer && bitpix == 32;
311 switch (algorithm) {
312 case NONE:
313 return ImageScale(detail::Bitpix<T>::value, 1.0, Bzero<T>::value);
314 case RANGE:
315 return determineFromRange(image, mask, isUnsigned, cfitsioPadding);
316 case MANUAL:
317 return ImageScale(bitpix, bscale, bzero);
321 return determineFromStdev(image, mask, isUnsigned, cfitsioPadding);
322 default:
323 std::abort(); // should never get here
324 }
325}
326#endif
327
328namespace {
329
334class CfitsioRandom {
335public:
337 CfitsioRandom(int seed) : _seed(seed) {
338 assert(seed != 0 && seed < N_RANDOM);
339 // Ensure that cfitsio has the correct locks initialized
340 // prior to initializing the random number table.
341 if (fitsio_init_lock()) {
342 throw LSST_EXCEPT(pex::exceptions::RuntimeError,
343 "Failed to initialize cfitsio locks for random table");
344 }
345 if (fits_init_randoms()) {
346 throw LSST_EXCEPT(pex::exceptions::RuntimeError,
347 "Failed to initialize cfitsio random table");
348 }
349 resetForTile(0);
350 }
351
353 void resetForTile(int iTile) {
354 _start = (iTile + _seed - 1) % N_RANDOM;
355 reseed();
356 }
357
359 float getNext() {
360 float const value = fits_rand_value[_index];
361 increment();
362 return value;
363 }
364
366 void increment() {
367 ++_index;
368 if (_index == N_RANDOM) {
369 ++_start;
370 if (_start == N_RANDOM) {
371 _start = 0;
372 }
373 reseed();
374 }
375 }
376
378 template <typename T>
379 ndarray::Array<T, 1, 1> forImage(typename ndarray::Array<T const, 2, 2>::Index const& shape,
380 ndarray::Array<long, 1> const& tiles) {
381 std::size_t const xSize = shape[1], ySize = shape[0];
382 ndarray::Array<T, 1, 1> out = ndarray::allocate(xSize * ySize);
383 std::size_t const xTileSize = tiles[0] <= 0 ? xSize : tiles[0];
384 std::size_t const yTileSize = tiles[1] < 0 ? ySize : (tiles[1] == 0 ? 1 : tiles[1]);
385 int const xNumTiles = std::ceil(xSize / static_cast<float>(xTileSize));
386 int const yNumTiles = std::ceil(ySize / static_cast<float>(yTileSize));
387 for (int iTile = 0, yTile = 0; yTile < yNumTiles; ++yTile) {
388 int const yStart = yTile * yTileSize;
389 int const yStop = std::min(yStart + yTileSize, ySize);
390 for (int xTile = 0; xTile < xNumTiles; ++xTile, ++iTile) {
391 int const xStart = xTile * xTileSize;
392 int const xStop = std::min(xStart + xTileSize, xSize);
393 resetForTile(iTile);
394 for (int y = yStart; y < yStop; ++y) {
395 auto iter = out.begin() + y * xSize + xStart;
396 for (int x = xStart; x < xStop; ++x, ++iter) {
397 *iter = static_cast<T>(getNext());
398 }
399 }
400 }
401 }
402 return out;
403 }
404
405private:
407 void reseed() { _index = static_cast<int>(fits_rand_value[_start] * 500); }
408
409 int _seed; // Initial seed
410 int _start; // Starting index for tile; "iseed" in cfitsio
411 int _index; // Index of next value; "nextrand" in cfitsio
412};
413
414} // anonymous namespace
415
416template <typename T>
418 bool forceNonfiniteRemoval, bool fuzz,
419 ndarray::Array<long, 1> const& tiles,
420 int seed) const {
424 "Floating-point images may not be converted to different floating-point types");
425 }
426 if (bscale != 1.0 || bzero != 0.0) {
428 "Scaling may not be applied to floating-point images");
429 }
430 }
431
432 if (bitpix < 0 || (bitpix == 0 && !std::numeric_limits<T>::is_integer) ||
433 (bscale == 1.0 && bzero == 0.0 && !fuzz)) {
434 if (!forceNonfiniteRemoval) {
435 // Type conversion only
436 return detail::makePixelArray(bitpix, ndarray::Array<T const, 1, 1>(ndarray::flatten<1>(image)));
437 }
439 ndarray::Array<T, 1, 1> out = ndarray::allocate(image.getNumElements());
440 auto outIter = out.begin();
441 auto const& flatImage = ndarray::flatten<1>(image);
442 for (auto inIter = flatImage.begin(); inIter != flatImage.end(); ++inIter, ++outIter) {
443 *outIter = std::isfinite(*inIter) ? *inIter : std::numeric_limits<T>::max();
444 }
445 return detail::makePixelArray(bitpix, out);
446 }
447 // Fall through for explicit scaling
448 }
449
450 // Note: BITPIX=8 treated differently, since it uses unsigned values; the rest use signed */
451 double min = bitpix == 8 ? 0 : -std::pow(2.0, bitpix - 1);
452 double max = bitpix == 8 ? 255 : (std::pow(2.0, bitpix - 1) - 1.0);
453
455 // cfitsio saves space for N_RESERVED_VALUES=10 values at the low end
457 }
458
459 double const scale = 1.0 / bscale;
460 std::size_t const num = image.getNumElements();
461 bool const applyFuzz = fuzz && !std::numeric_limits<T>::is_integer && bitpix > 0;
462 ndarray::Array<double, 1, 1> out;
463 if (applyFuzz) {
464 if (tiles.isEmpty()) {
466 "Tile sizes must be provided if fuzzing is desired");
467 }
468 out = CfitsioRandom(seed).forImage<double>(image.getShape(), tiles);
469 } else {
470 out = ndarray::allocate(num);
471 out.deep() = 0;
472 }
473 auto outIter = out.begin();
474 auto const& flatImage = ndarray::flatten<1>(image);
475 for (auto inIter = flatImage.begin(); inIter != flatImage.end(); ++inIter, ++outIter) {
476 double value = (*inIter - bzero) * scale;
477 if (!std::isfinite(value)) {
478 // This choice of "max" for non-finite and overflow pixels is mainly cosmetic --- it has to be
479 // something, and "min" would produce holes in the cores of bright stars.
480 *outIter = blank;
481 continue;
482 }
483 if (applyFuzz) {
484 // Add random factor [0.0,1.0): adds a variance of 1/12,
485 // but preserves the expectation value given the floor()
486 value += *outIter;
487 }
488 *outIter = (value < min ? blank : (value > max ? blank : std::floor(value)));
489 }
490 return detail::makePixelArray(bitpix, out);
491}
492
493template <typename T>
494ndarray::Array<T, 2, 2> ImageScale::fromFits(ndarray::Array<T, 2, 2> const& image) const {
495 ndarray::Array<T, 2, 2> memory = ndarray::allocate(image.getShape());
496 memory.deep() = bscale * image + bzero;
497 return memory;
498}
499
500// Explicit instantiation
501#define INSTANTIATE(TYPE) \
502 template ImageScale ImageScalingOptions::determine<TYPE, 2>( \
503 ndarray::Array<TYPE const, 2, 2> const& image, ndarray::Array<bool, 2, 2> const& mask) const; \
504 template std::shared_ptr<detail::PixelArrayBase> ImageScale::toFits<TYPE>( \
505 ndarray::Array<TYPE const, 2, 2> const&, bool, bool, ndarray::Array<long, 1> const&, int) const; \
506 template ndarray::Array<TYPE, 2, 2> ImageScale::fromFits<TYPE>(ndarray::Array<TYPE, 2, 2> const&) const;
507
515INSTANTIATE(boost::float32_t);
516INSTANTIATE(boost::float64_t);
517
518} // namespace fits
519} // namespace afw
520} // namespace lsst
int min
int max
#define INSTANTIATE(FROMSYS, TOSYS)
Definition Detector.cc:509
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< MaskPixelT > > mask
std::ostream * os
Definition Schema.cc:557
int y
Definition SpanSet.cc:48
T abort(T... args)
T ceil(T... args)
ScalingAlgorithm algorithm
Scaling algorithm to use.
float quantizeLevel
Divisor of the standard deviation for STDEV_* scaling.
@ STDEV_NEGATIVE
Scale based on the standard deviation, dynamic range negative.
@ STDEV_POSITIVE
Scale based on the standard deviation. dynamic range positive.
@ STDEV_BOTH
Scale based on the standard deviation, dynamic range positive+negative.
@ RANGE
Scale to preserve dynamic range.
ImageScale determine(image::ImageBase< T > const &image, image::Mask< image::MaskPixel > const *mask=nullptr) const
Determine the scaling for a particular image.
float quantizePad
Number of stdev to allow on the low/high side (for STDEV_POSITIVE/NEGATIVE)
int bitpix
Bits per pixel (0, 8,16,32,64,-32,-64)
double bscale
Manually specified BSCALE (for MANUAL scaling)
double bzero
Manually specified BZERO (for MANUAL scaling)
Reports invalid arguments.
Definition Runtime.h:66
float * fits_rand_value
int const N_RESERVED_VALUES
T floor(T... args)
T isfinite(T... args)
T make_pair(T... args)
T max(T... args)
T min(T... args)
std::shared_ptr< PixelArrayBase > makePixelArray(int bitpix, ndarray::Array< T, 1, 1 > const &array)
Create a PixelArray suitable for an image with the nominated BITPIX.
std::string compressionAlgorithmToString(ImageCompressionOptions::CompressionAlgorithm algorithm)
Provide string version of compression algorithm.
ImageScalingOptions::ScalingAlgorithm scalingAlgorithmFromString(std::string const &name)
Interpret scaling algorithm expressed in string.
int compressionAlgorithmToCfitsio(ImageCompressionOptions::CompressionAlgorithm algorithm)
Convert ImageCompressionOptions::CompressionAlgorithm to cfitsio.
std::string scalingAlgorithmToString(ImageScalingOptions::ScalingAlgorithm algorithm)
Provide string version of compression algorithm.
ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromCfitsio(int cfitsio)
Convert compression algorithm from cfitsio to ImageCompressionOptions::CompressionAlgorithm.
ImageCompressionOptions::CompressionAlgorithm compressionAlgorithmFromString(std::string const &name)
Interpret compression algorithm expressed in string.
STL namespace.
T nth_element(T... args)
T pow(T... args)
Scaling zero-point, set according to pixel type.
static double constexpr value
ImageCompressionOptions(CompressionAlgorithm algorithm_, Tiles tiles_, float quantizeLevel_=0.0)
Custom compression.
Tiles tiles
Tile size; a dimension with 0 means infinite (e.g., to specify one row: 0,1)
@ GZIP_SHUFFLE
GZIP compression with shuffle (most-significant byte first)
Scale to apply to image.
double bscale
Scale to apply when reading from FITS.
long blank
Value for integer images indicating non-finite values.
int bitpix
Bits per pixel; negative means floating-point: 8,16,32,64,-32,-64.
double bzero
Zero-point to apply when reading from FITS.
ndarray::Array< T, 2, 2 > fromFits(ndarray::Array< T, 2, 2 > const &image) const
Convert to an array.
std::shared_ptr< detail::PixelArrayBase > toFits(ndarray::Array< T const, 2, 2 > const &image, bool forceNonfiniteRemoval, bool fuzz=true, ndarray::Array< long, 1 > const &tiles=ndarray::Array< long, 1, 1 >(), int seed=1) const
Convert to an array of pixel values to write to FITS.
FITS BITPIX header value by C++ type.