1/* 2 * FITS image decoder 3 * Copyright (c) 2017 Paras Chadha 4 * 5 * This file is part of FFmpeg. 6 * 7 * FFmpeg is free software; you can redistribute it and/or 8 * modify it under the terms of the GNU Lesser General Public 9 * License as published by the Free Software Foundation; either 10 * version 2.1 of the License, or (at your option) any later version. 11 * 12 * FFmpeg is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 * Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public 18 * License along with FFmpeg; if not, write to the Free Software 19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 20 */ 21 22/** 23 * @file 24 * FITS image decoder 25 * 26 * Specification: https://fits.gsfc.nasa.gov/fits_standard.html Version 3.0 27 * 28 * Support all 2d images alongwith, bzero, bscale and blank keywords. 29 * RGBA images are supported as NAXIS3 = 3 or 4 i.e. Planes in RGBA order. Also CTYPE = 'RGB ' should be present. 30 * Also to interpret data, values are linearly scaled using min-max scaling but not RGB images. 31 */ 32 33#include "avcodec.h" 34#include "codec_internal.h" 35#include "internal.h" 36#include <float.h> 37#include "libavutil/intreadwrite.h" 38#include "libavutil/intfloat.h" 39#include "libavutil/dict.h" 40#include "libavutil/opt.h" 41#include "fits.h" 42 43typedef struct FITSContext { 44 const AVClass *class; 45 int blank_val; 46} FITSContext; 47 48/** 49 * Calculate the data_min and data_max values from the data. 50 * This is called if the values are not present in the header. 51 * @param ptr8 pointer to the data 52 * @param header pointer to the header 53 * @param end pointer to end of packet 54 * @return 0 if calculated successfully otherwise AVERROR_INVALIDDATA 55 */ 56static int fill_data_min_max(const uint8_t *ptr8, FITSHeader *header, const uint8_t *end) 57{ 58 uint8_t t8; 59 int16_t t16; 60 int32_t t32; 61 int64_t t64; 62 float tflt; 63 double tdbl; 64 int i, j; 65 66 header->data_min = DBL_MAX; 67 header->data_max = -DBL_MAX; 68 switch (header->bitpix) { 69#define CASE_N(a, t, rd) \ 70 case a: \ 71 for (i = 0; i < header->naxisn[1]; i++) { \ 72 for (j = 0; j < header->naxisn[0]; j++) { \ 73 t = rd; \ 74 if (!header->blank_found || t != header->blank) { \ 75 if (t > header->data_max) \ 76 header->data_max = t; \ 77 if (t < header->data_min) \ 78 header->data_min = t; \ 79 } \ 80 ptr8 += abs(a) >> 3; \ 81 } \ 82 } \ 83 break 84 85 CASE_N(-64, tdbl, av_int2double(AV_RB64(ptr8))); 86 CASE_N(-32, tflt, av_int2float(AV_RB32(ptr8))); 87 CASE_N(8, t8, ptr8[0]); 88 CASE_N(16, t16, AV_RB16(ptr8)); 89 CASE_N(32, t32, AV_RB32(ptr8)); 90 CASE_N(64, t64, AV_RB64(ptr8)); 91 default: 92 return AVERROR_INVALIDDATA; 93 } 94 return 0; 95} 96 97/** 98 * Read the fits header and store the values in FITSHeader pointed by header 99 * @param avctx AVCodec context 100 * @param ptr pointer to pointer to the data 101 * @param header pointer to the FITSHeader 102 * @param end pointer to end of packet 103 * @param metadata pointer to pointer to AVDictionary to store metadata 104 * @return 0 if calculated successfully otherwise AVERROR_INVALIDDATA 105 */ 106static int fits_read_header(AVCodecContext *avctx, const uint8_t **ptr, FITSHeader *header, 107 const uint8_t *end, AVDictionary **metadata) 108{ 109 const uint8_t *ptr8 = *ptr; 110 int lines_read, bytes_left, i, ret; 111 size_t size; 112 113 lines_read = 1; // to account for first header line, SIMPLE or XTENSION which is not included in packet... 114 avpriv_fits_header_init(header, STATE_BITPIX); 115 do { 116 if (end - ptr8 < 80) 117 return AVERROR_INVALIDDATA; 118 ret = avpriv_fits_header_parse_line(avctx, header, ptr8, &metadata); 119 ptr8 += 80; 120 lines_read++; 121 } while (!ret); 122 if (ret < 0) 123 return ret; 124 125 bytes_left = (((lines_read + 35) / 36) * 36 - lines_read) * 80; 126 if (end - ptr8 < bytes_left) 127 return AVERROR_INVALIDDATA; 128 ptr8 += bytes_left; 129 130 if (header->rgb && (header->naxis != 3 || (header->naxisn[2] != 3 && header->naxisn[2] != 4))) { 131 av_log(avctx, AV_LOG_ERROR, "File contains RGB image but NAXIS = %d and NAXIS3 = %d\n", header->naxis, header->naxisn[2]); 132 return AVERROR_INVALIDDATA; 133 } 134 135 if (!header->rgb && header->naxis != 2) { 136 av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions, NAXIS = %d\n", header->naxis); 137 return AVERROR_INVALIDDATA; 138 } 139 140 if (header->blank_found && (header->bitpix == -32 || header->bitpix == -64)) { 141 av_log(avctx, AV_LOG_WARNING, "BLANK keyword found but BITPIX = %d\n. Ignoring BLANK", header->bitpix); 142 header->blank_found = 0; 143 } 144 145 size = abs(header->bitpix) >> 3; 146 for (i = 0; i < header->naxis; i++) { 147 if (size == 0 || header->naxisn[i] > SIZE_MAX / size) { 148 av_log(avctx, AV_LOG_ERROR, "unsupported size of FITS image"); 149 return AVERROR_INVALIDDATA; 150 } 151 size *= header->naxisn[i]; 152 } 153 154 if (end - ptr8 < size) 155 return AVERROR_INVALIDDATA; 156 *ptr = ptr8; 157 158 if (!header->rgb && (!header->data_min_found || !header->data_max_found)) { 159 ret = fill_data_min_max(ptr8, header, end); 160 if (ret < 0) { 161 av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header->bitpix); 162 return ret; 163 } 164 } else { 165 /* 166 * instead of applying bscale and bzero to every element, 167 * we can do inverse transformation on data_min and data_max 168 */ 169 header->data_min = (header->data_min - header->bzero) / header->bscale; 170 header->data_max = (header->data_max - header->bzero) / header->bscale; 171 } 172 if (!header->rgb && header->data_min >= header->data_max) { 173 if (header->data_min > header->data_max) { 174 av_log(avctx, AV_LOG_ERROR, "data min/max (%g %g) is invalid\n", header->data_min, header->data_max); 175 return AVERROR_INVALIDDATA; 176 } 177 av_log(avctx, AV_LOG_WARNING, "data min/max indicates a blank image\n"); 178 header->data_max ++; 179 } 180 181 return 0; 182} 183 184static int fits_decode_frame(AVCodecContext *avctx, AVFrame *p, 185 int *got_frame, AVPacket *avpkt) 186{ 187 const uint8_t *ptr8 = avpkt->data, *end; 188 uint8_t t8; 189 int16_t t16; 190 int32_t t32; 191 int64_t t64; 192 float tflt; 193 double tdbl; 194 int ret, i, j, k; 195 const int map[] = {2, 0, 1, 3}; // mapping from GBRA -> RGBA as RGBA is to be stored in FITS file.. 196 uint8_t *dst8; 197 uint16_t *dst16; 198 uint64_t t; 199 FITSHeader header; 200 FITSContext * fitsctx = avctx->priv_data; 201 202 end = ptr8 + avpkt->size; 203 p->metadata = NULL; 204 ret = fits_read_header(avctx, &ptr8, &header, end, &p->metadata); 205 if (ret < 0) 206 return ret; 207 208 if (header.rgb) { 209 if (header.bitpix == 8) { 210 if (header.naxisn[2] == 3) { 211 avctx->pix_fmt = AV_PIX_FMT_GBRP; 212 } else { 213 avctx->pix_fmt = AV_PIX_FMT_GBRAP; 214 } 215 } else if (header.bitpix == 16) { 216 if (header.naxisn[2] == 3) { 217 avctx->pix_fmt = AV_PIX_FMT_GBRP16; 218 } else { 219 avctx->pix_fmt = AV_PIX_FMT_GBRAP16; 220 } 221 } else { 222 av_log(avctx, AV_LOG_ERROR, "unsupported BITPIX = %d\n", header.bitpix); 223 return AVERROR_INVALIDDATA; 224 } 225 } else { 226 if (header.bitpix == 8) { 227 avctx->pix_fmt = AV_PIX_FMT_GRAY8; 228 } else { 229 avctx->pix_fmt = AV_PIX_FMT_GRAY16; 230 } 231 } 232 233 if ((ret = ff_set_dimensions(avctx, header.naxisn[0], header.naxisn[1])) < 0) 234 return ret; 235 236 if ((ret = ff_get_buffer(avctx, p, 0)) < 0) 237 return ret; 238 239 /* 240 * FITS stores images with bottom row first. Therefore we have 241 * to fill the image from bottom to top. 242 */ 243 if (header.rgb) { 244 switch(header.bitpix) { 245#define CASE_RGB(cas, dst, type, dref) \ 246 case cas: \ 247 for (k = 0; k < header.naxisn[2]; k++) { \ 248 for (i = 0; i < avctx->height; i++) { \ 249 dst = (type *) (p->data[map[k]] + (avctx->height - i - 1) * p->linesize[map[k]]); \ 250 for (j = 0; j < avctx->width; j++) { \ 251 t32 = dref(ptr8); \ 252 if (!header.blank_found || t32 != header.blank) { \ 253 t = t32 * header.bscale + header.bzero; \ 254 } else { \ 255 t = fitsctx->blank_val; \ 256 } \ 257 *dst++ = (type) t; \ 258 ptr8 += cas >> 3; \ 259 } \ 260 } \ 261 } \ 262 break 263 264 CASE_RGB(8, dst8, uint8_t, *); 265 CASE_RGB(16, dst16, uint16_t, AV_RB16); 266 } 267 } else { 268 double scale = header.data_max - header.data_min; 269 270 if (scale <= 0 || !isfinite(scale)) { 271 scale = 1; 272 } 273 scale = 1/scale; 274 275 switch (header.bitpix) { 276#define CASE_GRAY(cas, dst, type, t, rd) \ 277 case cas: \ 278 for (i = 0; i < avctx->height; i++) { \ 279 dst = (type *) (p->data[0] + (avctx->height-i-1)* p->linesize[0]); \ 280 for (j = 0; j < avctx->width; j++) { \ 281 t = rd; \ 282 if (!header.blank_found || t != header.blank) { \ 283 *dst++ = lrint(((t - header.data_min) * ((1 << (sizeof(type) * 8)) - 1)) * scale); \ 284 } else { \ 285 *dst++ = fitsctx->blank_val; \ 286 } \ 287 ptr8 += abs(cas) >> 3; \ 288 } \ 289 } \ 290 break 291 292 CASE_GRAY(-64, dst16, uint16_t, tdbl, av_int2double(AV_RB64(ptr8))); 293 CASE_GRAY(-32, dst16, uint16_t, tflt, av_int2float(AV_RB32(ptr8))); 294 CASE_GRAY(8, dst8, uint8_t, t8, ptr8[0]); 295 CASE_GRAY(16, dst16, uint16_t, t16, AV_RB16(ptr8)); 296 CASE_GRAY(32, dst16, uint16_t, t32, AV_RB32(ptr8)); 297 CASE_GRAY(64, dst16, uint16_t, t64, AV_RB64(ptr8)); 298 default: 299 av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header.bitpix); 300 return AVERROR_INVALIDDATA; 301 } 302 } 303 304 p->key_frame = 1; 305 p->pict_type = AV_PICTURE_TYPE_I; 306 307 *got_frame = 1; 308 309 return avpkt->size; 310} 311 312static const AVOption fits_options[] = { 313 { "blank_value", "value that is used to replace BLANK pixels in data array", offsetof(FITSContext, blank_val), AV_OPT_TYPE_INT, { .i64 = 0 }, 0, 65535, AV_OPT_FLAG_DECODING_PARAM | AV_OPT_FLAG_VIDEO_PARAM}, 314 { NULL }, 315}; 316 317static const AVClass fits_decoder_class = { 318 .class_name = "FITS decoder", 319 .item_name = av_default_item_name, 320 .option = fits_options, 321 .version = LIBAVUTIL_VERSION_INT, 322}; 323 324const FFCodec ff_fits_decoder = { 325 .p.name = "fits", 326 .p.type = AVMEDIA_TYPE_VIDEO, 327 .p.id = AV_CODEC_ID_FITS, 328 .p.capabilities = AV_CODEC_CAP_DR1, 329 .p.long_name = NULL_IF_CONFIG_SMALL("Flexible Image Transport System"), 330 .p.priv_class = &fits_decoder_class, 331 .priv_data_size = sizeof(FITSContext), 332 FF_CODEC_DECODE_CB(fits_decode_frame), 333}; 334