diff mbox

[FFmpeg-devel] Add FITS Decoder

Message ID 1498802412-24766-1-git-send-email-paraschadha18@gmail.com
State Superseded
Headers show

Commit Message

Paras June 30, 2017, 6 a.m. UTC
Made the changes suggested above

Signed-off-by: Paras Chadha <paraschadha18@gmail.com>
---
 Changelog               |   1 +
 doc/general.texi        |   2 +
 libavcodec/Makefile     |   1 +
 libavcodec/allcodecs.c  |   1 +
 libavcodec/avcodec.h    |   1 +
 libavcodec/codec_desc.c |   8 +
 libavcodec/fitsdec.c    | 580 ++++++++++++++++++++++++++++++++++++++++++++++++
 libavcodec/version.h    |   4 +-
 libavformat/img2.c      |   1 +
 9 files changed, 597 insertions(+), 2 deletions(-)
 create mode 100644 libavcodec/fitsdec.c

--
2.4.11

Comments

Nicolas George June 30, 2017, 1:08 p.m. UTC | #1
Hi. A few technical / cosmetic remarks below. I do not know the FITS
format itself more than in passing.

Le duodi 12 messidor, an CCXXV, Paras Chadha a écrit :
> Made the changes suggested above
> 
> Signed-off-by: Paras Chadha <paraschadha18@gmail.com>
> ---
>  Changelog               |   1 +
>  doc/general.texi        |   2 +
>  libavcodec/Makefile     |   1 +
>  libavcodec/allcodecs.c  |   1 +
>  libavcodec/avcodec.h    |   1 +
>  libavcodec/codec_desc.c |   8 +
>  libavcodec/fitsdec.c    | 580 ++++++++++++++++++++++++++++++++++++++++++++++++
>  libavcodec/version.h    |   4 +-
>  libavformat/img2.c      |   1 +
>  9 files changed, 597 insertions(+), 2 deletions(-)
>  create mode 100644 libavcodec/fitsdec.c
> 
> diff --git a/Changelog b/Changelog
> index a8726c6..2c2bdec 100644
> --- a/Changelog
> +++ b/Changelog
> @@ -26,6 +26,7 @@ version <next>:
>    --x86asmexe=yasm to configure to restore the old behavior.
>  - additional frame format support for Interplay MVE movies
>  - support for decoding through D3D11VA in ffmpeg
> +- FITS demuxer and decoder
> 
>  version 3.3:
>  - CrystalHD decoder moved to new decode API
> diff --git a/doc/general.texi b/doc/general.texi
> index 8f582d5..c00ce32 100644
> --- a/doc/general.texi
> +++ b/doc/general.texi
> @@ -591,6 +591,8 @@ following image formats are supported:
>      @tab Digital Picture Exchange
>  @item EXR          @tab   @tab X
>      @tab OpenEXR
> +@item FITS         @tab   @tab X
> +    @tab Flexible Image Transport System
>  @item JPEG         @tab X @tab X
>      @tab Progressive JPEG is not supported.
>  @item JPEG 2000    @tab X @tab X
> diff --git a/libavcodec/Makefile b/libavcodec/Makefile
> index b440a00..729e95e 100644
> --- a/libavcodec/Makefile
> +++ b/libavcodec/Makefile
> @@ -291,6 +291,7 @@ OBJS-$(CONFIG_FFV1_DECODER)            += ffv1dec.o ffv1.o
>  OBJS-$(CONFIG_FFV1_ENCODER)            += ffv1enc.o ffv1.o
>  OBJS-$(CONFIG_FFWAVESYNTH_DECODER)     += ffwavesynth.o
>  OBJS-$(CONFIG_FIC_DECODER)             += fic.o
> +OBJS-$(CONFIG_FITS_DECODER)            += fitsdec.o
>  OBJS-$(CONFIG_FLAC_DECODER)            += flacdec.o flacdata.o flac.o
>  OBJS-$(CONFIG_FLAC_ENCODER)            += flacenc.o flacdata.o flac.o vorbis_data.o
>  OBJS-$(CONFIG_FLASHSV_DECODER)         += flashsv.o
> diff --git a/libavcodec/allcodecs.c b/libavcodec/allcodecs.c
> index 0243f47..a4cfd80 100644
> --- a/libavcodec/allcodecs.c
> +++ b/libavcodec/allcodecs.c
> @@ -192,6 +192,7 @@ static void register_all(void)
>      REGISTER_ENCDEC (FFV1,              ffv1);
>      REGISTER_ENCDEC (FFVHUFF,           ffvhuff);
>      REGISTER_DECODER(FIC,               fic);
> +    REGISTER_DECODER(FITS,              fits);
>      REGISTER_ENCDEC (FLASHSV,           flashsv);
>      REGISTER_ENCDEC (FLASHSV2,          flashsv2);
>      REGISTER_DECODER(FLIC,              flic);
> diff --git a/libavcodec/avcodec.h b/libavcodec/avcodec.h
> index b697afa..8eba460 100644
> --- a/libavcodec/avcodec.h
> +++ b/libavcodec/avcodec.h
> @@ -447,6 +447,7 @@ enum AVCodecID {
>      AV_CODEC_ID_SRGC,
>      AV_CODEC_ID_SVG,
>      AV_CODEC_ID_GDV,
> +    AV_CODEC_ID_FITS,
> 
>      /* various PCM "codecs" */
>      AV_CODEC_ID_FIRST_AUDIO = 0x10000,     ///< A dummy id pointing at the start of audio codecs
> diff --git a/libavcodec/codec_desc.c b/libavcodec/codec_desc.c
> index cf1246e..0112517 100644
> --- a/libavcodec/codec_desc.c
> +++ b/libavcodec/codec_desc.c
> @@ -1464,6 +1464,14 @@ static const AVCodecDescriptor codec_descriptors[] = {
>                       AV_CODEC_PROP_LOSSLESS,
>      },
>      {
> +        .id        = AV_CODEC_ID_FITS,
> +        .type      = AVMEDIA_TYPE_VIDEO,
> +        .name      = "fits",
> +        .long_name = NULL_IF_CONFIG_SMALL("Flexible Image Transport System"),
> +        .props     = AV_CODEC_PROP_INTRA_ONLY | AV_CODEC_PROP_LOSSY |
> +                     AV_CODEC_PROP_LOSSLESS,
> +    },
> +    {
>          .id        = AV_CODEC_ID_GIF,
>          .type      = AVMEDIA_TYPE_VIDEO,
>          .name      = "gif",
> diff --git a/libavcodec/fitsdec.c b/libavcodec/fitsdec.c
> new file mode 100644
> index 0000000..4eaf3c8
> --- /dev/null
> +++ b/libavcodec/fitsdec.c
> @@ -0,0 +1,580 @@
> +/*
> + * FITS image decoder
> + * Copyright (c) 2017 Paras Chadha
> + *
> + * This file is part of FFmpeg.
> + *
> + * FFmpeg is free software; you can redistribute it and/or
> + * modify it under the terms of the GNU Lesser General Public
> + * License as published by the Free Software Foundation; either
> + * version 2.1 of the License, or (at your option) any later version.
> + *
> + * FFmpeg is distributed in the hope that it will be useful,
> + * but WITHOUT ANY WARRANTY; without even the implied warranty of
> + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> + * Lesser General Public License for more details.
> + *
> + * You should have received a copy of the GNU Lesser General Public
> + * License along with FFmpeg; if not, write to the Free Software
> + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
> + */
> +
> +/**

> + * @file
> + * FITS image decoder
> + * It supports all 2-d images alongwith, bzero, bscale and blank keywords.
> + * RGBA images are supported as NAXIS3 = 3 or 4 i.e. Planes in RGBA order. Also CTYPE = 'RGB ' should be present.
> + * Also to interpret data, values are linearly scaled using min-max scaling but not RGB images.

Please add a link to the specification you followed.

> + */
> +
> +#include "avcodec.h"
> +#include "internal.h"
> +#include <float.h>
> +#include "libavutil/intreadwrite.h"
> +#include "libavutil/intfloat.h"
> +#include "libavutil/dict.h"
> +
> +/**
> + * Structure to store the header keywords in FITS file
> + */
> +typedef struct FITSContext {

> +    char simple;

This field is only used in fits_read_header(), it should be a local
variable.

> +    int bitpix;
> +    int64_t blank;

> +    int naxis;

Ditto.

> +    int naxisn[3];
> +    int rgb; /**< 1 if file contains RGB image, 0 otherwise */

> +    int xtension;

Ditto.

> +    double bscale;
> +    double bzero;
> +    double data_min;
> +    double data_max;
> +} FITSDecContext;

In fact, I believe not a single of these fields has a value that
survives the decode_frame() call: it the whole header structure should
be a local variable in decode_frame().

> +
> +/**

> + * function calculates the data_min and data_max values from the data.

The style for doxygen comments is impersonal verbs: "Calculate the ...".

> + * This is called if the values are not present in the header.

> + * @param ptr8 - pointer to the data
> + * @param header - pointer to the header
> + * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA

Please drop the dashes and comma. The last two comments apply for all
doxygen comments.

And 0 is usually used for success.

> + */
> +static int fill_data_min_max(const uint8_t * ptr8, FITSDecContext * header, const uint8_t * end)
> +{
> +    int16_t t16;
> +    int32_t t32;
> +    int64_t t64;
> +    float tflt;
> +    double tdbl;
> +    int i, j;
> +
> +    header->data_min = DBL_MAX;
> +    header->data_max = DBL_MIN;
> +    switch (header->bitpix) {
> +        case -64:

> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    tdbl = av_int2double(AV_RB64(ptr8));
> +                    if (tdbl > header->data_max)
> +                        header->data_max = tdbl;
> +                    if (tdbl < header->data_min)
> +                        header->data_min = tdbl;
> +                    ptr8 += 8;
> +                }
> +            }
> +            break;
> +        case -32:
> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    tflt = av_int2float(AV_RB32(ptr8));
> +                    if (tflt > header->data_max)
> +                        header->data_max = tflt;
> +                    if (tflt < header->data_min)
> +                        header->data_min = tflt;
> +                    ptr8 += 4;
> +                }
> +            }
> +            break;
> +        case 8:
> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    if (ptr8[0] != header->blank) {
> +                        if (ptr8[0] > header->data_max)
> +                            header->data_max = ptr8[0];
> +                        if (ptr8[0] < header->data_min)
> +                            header->data_min = ptr8[0];
> +                    }
> +                    ptr8++;
> +                }
> +            }
> +            break;
> +        case 16:
> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    t16 = AV_RB16(ptr8);
> +                    if (t16 != header->blank) {
> +                        if (t16 > header->data_max)
> +                            header->data_max = t16;
> +                        if (t16 < header->data_min)
> +                            header->data_min = t16;
> +                    }
> +                    ptr8 += 2;
> +                }
> +            }
> +            break;
> +        case 32:
> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    t32 = AV_RB32(ptr8);
> +                    if (t32 != header->blank) {
> +                        if (t32 > header->data_max)
> +                            header->data_max = t32;
> +                        if (t32 < header->data_min)
> +                            header->data_min = t32;
> +                    }
> +                    ptr8 += 4;
> +                }
> +            }
> +            break;
> +        case 64:
> +            for (i = 0; i < header->naxisn[1]; i++) {
> +                for (j = 0; j < header->naxisn[0]; j++) {
> +                    t64 = AV_RB64(ptr8);
> +                    if (t64 != header->blank) {
> +                        if (t64 > header->data_max)
> +                            header->data_max = t64;
> +                        if (t64 < header->data_min)
> +                            header->data_min = t64;
> +                    }
> +                    ptr8 += 8;
> +                }
> +            }
> +            break;

Please find a way of factoring these almost-identical loops. A macro
would be the obvious choice.

> +        default:
> +            return AVERROR_INVALIDDATA;
> +    }
> +    return 1;
> +}
> +
> +/**
> + * function reads the fits header and stores the values in FITSDecContext pointed by header
> + * @param avctx - AVCodec context
> + * @param ptr - pointer to pointer to the data
> + * @param header - pointer to the FITSDecContext
> + * @param end - pointer to end of packet
> + * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA
> + */
> +static int fits_read_header(AVCodecContext *avctx, const uint8_t **ptr, FITSDecContext * header,
> +                            const uint8_t * end, AVDictionary **meta)
> +{
> +    const uint8_t *ptr8 = *ptr;
> +    int lines_read = 0, i, dim_no, t, data_min_found = 0, data_max_found = 0, ret;
> +    uint64_t size=1;
> +    double d;
> +    AVDictionary *metadata = NULL;
> +    char keyword[10], value[72];
> +

> +    header->blank = LLONG_MIN;

This is almost certainly not the value you want.

> +    header->bscale = 1.0;
> +    header->bzero = 0;
> +    header->rgb = 0;
> +

> +    if (end - ptr8 < 80)
> +        return AVERROR_INVALIDDATA;
> +
> +    if (sscanf(ptr8, "SIMPLE = %c", &header->simple) == 1) {

You did not check that there is a NUL character. A crafted input could
cause an over-read, and relying on the padding is bad style. I suggest
you check that there is at least one padding space at the end (and make
sure the sscanf patterns do not). As a macro or a small helper function,
for example.

Another possibility, more strictly conforming to the standard, would be
to copy the line into a 81-sized buffer.

> +        if (header->simple == 'F') {
> +            av_log(avctx, AV_LOG_WARNING, "not a standard FITS file\n");
> +            av_dict_set(&metadata, "SIMPLE", "F", 0);

> +        } else if (header->simple != 'T') {
> +            av_log(avctx, AV_LOG_ERROR, "invalid SIMPLE value, SIMPLE = %c\n", header->simple);
> +            return AVERROR_INVALIDDATA;
> +        } else {
> +            av_dict_set(&metadata, "SIMPLE", "T", 0);
> +        }

I strongly suggest to swap the if and else clause.

> +        header->xtension = 0;
> +    } else if (!strncmp(ptr8, "XTENSION= 'IMAGE", 16)) {
> +        header->xtension = 1;
> +        av_dict_set(&metadata, "XTENSION", "'IMAGE   '", 0);
> +    } else {
> +        av_log(avctx, AV_LOG_ERROR, "missing SIMPLE keyword or invalid XTENSION\n");
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    ptr8 += 80;
> +    lines_read++;
> +
> +    if (end - ptr8 < 80)
> +        return AVERROR_INVALIDDATA;
> +
> +    if (sscanf(ptr8, "BITPIX = %d", &header->bitpix) != 1) {
> +        av_log(avctx, AV_LOG_ERROR, "missing BITPIX keyword\n");
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    av_dict_set_int(&metadata, "BITPIX", header->bitpix, 0);
> +    size = abs(header->bitpix) >> 3;
> +    ptr8 += 80;
> +    lines_read++;
> +
> +    if (end - ptr8 < 80)
> +        return AVERROR_INVALIDDATA;
> +
> +    if (sscanf(ptr8, "NAXIS = %d", &header->naxis) != 1) {
> +        av_log(avctx, AV_LOG_ERROR, "missing NAXIS keyword\n");
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    if (!header->naxis) {
> +        av_log(avctx, AV_LOG_ERROR, "No image data found, NAXIS = %d\n", header->naxis);
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    if (header->naxis != 2 && header->naxis != 3) {
> +        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions, NAXIS = %d\n", header->naxis);
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    av_dict_set_int(&metadata, "NAXIS", header->naxis, 0);
> +    ptr8 += 80;
> +    lines_read++;
> +
> +    for (i = 0; i < header->naxis; i++) {
> +        if (end - ptr8 < 80)
> +            return AVERROR_INVALIDDATA;
> +
> +        if (sscanf(ptr8, "NAXIS%d = %d", &dim_no, &header->naxisn[i]) != 2 || dim_no != i+1) {
> +            av_log(avctx, AV_LOG_ERROR, "missing NAXIS%d keyword\n", i+1);
> +            return AVERROR_INVALIDDATA;
> +        }
> +

> +        ret = snprintf(keyword, 10, "NAXIS%d", dim_no);
> +        if (ret < 0 || ret >= 10) {

Use sizeof(keyword).

> +            return AVERROR_INVALIDDATA;
> +        }
> +
> +        av_dict_set_int(&metadata, keyword, header->naxisn[i], 0);

> +        size *= header->naxisn[i];

This can overflow.

> +        ptr8 += 80;
> +        lines_read++;
> +    }
> +

> +    if (end - ptr8 < 80)
> +            return AVERROR_INVALIDDATA;

Strange indentation.

> +
> +    while (strncmp(ptr8, "END", 3)) {
> +        if (sscanf(ptr8, "BLANK = %d", &t) == 1) {
> +            header->blank = t;
> +        } else if (sscanf(ptr8, "BSCALE = %lf", &d) == 1) {
> +            header->bscale = d;
> +        } else if (sscanf(ptr8, "BZERO = %lf", &d) == 1) {
> +            header->bzero = d;
> +        } else if (sscanf(ptr8, "DATAMAX = %lf", &d) == 1) {
> +            data_max_found = 1;
> +            header->data_max = d;
> +        } else if (sscanf(ptr8, "DATAMIN = %lf", &d) == 1) {
> +            data_min_found = 1;
> +            header->data_min = d;
> +        } else if (!strncmp(ptr8, "CTYPE3  = 'RGB", 14)) {
> +            header->rgb = 1;
> +            if (header->naxis != 3 || (header->naxisn[2] != 3 && header->naxisn[2] != 4)) {
> +                av_log(avctx, AV_LOG_ERROR, "File contains RGB image but NAXIS = %d and NAXIS3 = %d\n", header->naxis, header->naxisn[2]);
> +                return AVERROR_INVALIDDATA;
> +            }
> +        }
> +
> +        if (ptr8[8] == '=') {
> +            for (i = 0; i < 8 && ptr8[i] != ' '; i++) {
> +                keyword[i] = ptr8[i];
> +            }
> +            keyword[i] = '\0';
> +
> +            t = 0;
> +            i = 10;
> +            while (i < 80 && ptr8[i] == ' ')
> +                i++;
> +
> +            if (i < 80) {
> +                value[t] = ptr8[i];
> +                i++;
> +                t++;
> +                if (ptr8[i-1] == '\'') {
> +                    while (i < 80 && ptr8[i] != '\'') {
> +                        value[t] = ptr8[i];
> +                        i++;
> +                        t++;
> +                    }
> +                    value[t] = '\'';
> +                    t++;
> +                } else if (ptr8[i-1] == '(') {
> +                    while (i < 80 && ptr8[i] != ')') {
> +                        value[t] = ptr8[i];
> +                        i++;
> +                        t++;
> +                    }
> +                    value[t] = ')';
> +                    t++;
> +                } else {
> +                    while (i < 80 && ptr8[i] != ' ' && ptr8[i] != '/') {
> +                        value[t] = ptr8[i];
> +                        i++;
> +                        t++;
> +                    }
> +                }
> +            }
> +
> +            value[t] = '\0';
> +            av_dict_set(&metadata, keyword, value, 0);
> +        }
> +
> +        ptr8 += 80;
> +        lines_read++;
> +
> +        if (end - ptr8 < 80)
> +            return AVERROR_INVALIDDATA;
> +    }
> +
> +    if (!header->rgb && header->naxis != 2){
> +        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions, NAXIS = %d\n", header->naxis);
> +        return AVERROR_INVALIDDATA;
> +    }
> +
> +    ptr8 += 80;
> +    lines_read++;
> +    lines_read %= 36;
> +
> +    t = ((36 - lines_read) % 36) * 80;
> +    if (end - ptr8 < t)
> +        return AVERROR_INVALIDDATA;
> +    ptr8 += t;
> +    *ptr = ptr8;
> +
> +    if (end - ptr8 < size)
> +        return AVERROR_INVALIDDATA;
> +

> +    if (!header->rgb && (!data_min_found || !data_max_found)) {
> +        if ((ret = fill_data_min_max(ptr8, header, end)) < 0) {

If only one of the min/max is provided in the file, fill_data_min_max()
will be called and override both. It is on purpose?

> +            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header->bitpix);

> +            return AVERROR_INVALIDDATA;

return ret;

> +        }
> +    } else {
> +        /*
> +         * instead of applying bscale and bzero to every element, we can do inverse transformation on data_min and
> +         * data_max
> +         */
> +        header->data_min = (header->data_min - header->bzero) / header->bscale;
> +        header->data_max = (header->data_max - header->bzero) / header->bscale;
> +    }
> +
> +    *meta = metadata;
> +    return 1;
> +}
> +
> +static int fits_decode_frame(AVCodecContext *avctx, void *data, int *got_frame, AVPacket *avpkt)
> +{
> +    AVFrame *p=data;
> +    const uint8_t *ptr8 = avpkt->data, *end;
> +    int16_t t16;
> +    int32_t t32;
> +    int64_t t64;
> +    float   tflt;
> +    double  tdbl;
> +    int ret, i, j;
> +    uint8_t *dst8;
> +    uint16_t *dst16;
> +    uint32_t *dst32;
> +    uint64_t *dst64, size, r, g, b, a, t;
> +    FITSDecContext * header = avctx->priv_data;
> +
> +    end = ptr8 + avpkt->size;
> +    if ((ret = fits_read_header(avctx, &ptr8, header, end, &p->metadata)) < 0)
> +        return ret;
> +
> +    size = (header->naxisn[0]) * (header->naxisn[1]);
> +
> +    if (header->rgb) {
> +        if (header->bitpix == 8) {
> +            avctx->pix_fmt = AV_PIX_FMT_RGB32;
> +        } else if (header->bitpix == 16) {
> +            avctx->pix_fmt = AV_PIX_FMT_RGBA64;
> +        } else {
> +            av_log(avctx, AV_LOG_ERROR, "unsupported BITPIX = %d\n", header->bitpix);
> +            return AVERROR_INVALIDDATA;
> +        }
> +    } else {
> +        if (header->bitpix == 8) {
> +            avctx->pix_fmt = AV_PIX_FMT_GRAY8;
> +        } else {
> +            avctx->pix_fmt = AV_PIX_FMT_GRAY16;
> +        }
> +    }
> +
> +    if ((ret = ff_set_dimensions(avctx, header->naxisn[0], header->naxisn[1])) < 0)
> +        return ret;
> +
> +    if ((ret = ff_get_buffer(avctx, p, 0)) < 0)
> +        return ret;
> +
> +    if (header->rgb) {

> +        if (header->bitpix == 8) {
> +            for (i = 0; i < avctx->height; i++) {
> +                /*
> +                 * FITS stores images with bottom row first. Therefore we have
> +                 * to fill the image from bottom to top.
> +                 */
> +                dst32 = (uint32_t *)(p->data[0] + (avctx->height-i-1)* p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    if (header->naxisn[2] == 4) {
> +                        if (ptr8[size * 3] != header->blank)
> +                            t = ptr8[size * 3] * header->bscale + header->bzero;
> +                        a = t << 24;
> +                    } else {
> +                        a = (255 << 24);
> +                    }
> +
> +                    if (ptr8[0] != header->blank)
> +                        t = ptr8[0] * header->bscale + header->bzero;
> +                    r = t << 16;
> +
> +                    if (ptr8[size] != header->blank)
> +                        t = ptr8[size] * header->bscale + header->bzero;
> +                    g = t << 8;
> +
> +                    if (ptr8[size * 2] != header->blank)
> +                        t = ptr8[size * 2] * header->bscale + header->bzero;
> +                    b = t;
> +
> +                    *dst32++ = ((uint32_t)a) | ((uint32_t)r) | ((uint32_t)g) | ((uint32_t)b);
> +                    ptr8++;
> +                }
> +            }
> +        } else if (header->bitpix == 16) {
> +            // not tested ....
> +            for (i = 0; i < avctx->height; i++) {
> +                dst64 = (uint64_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +
> +                    if (header->naxisn[2] == 4) {
> +                        t = ((ptr8[size * 3] << 8) | ptr8[size * 3 + 1]);
> +                        if (t != header->blank)
> +                            t = t*header->bscale + header->bzero;
> +                        a = t << 48;
> +                    } else {
> +                        a = 65535ULL << 48;
> +                    }
> +
> +                    t = ptr8[0] << 8 | ptr8[1];
> +                    if (t != header->blank)
> +                        t = t*header->bscale + header->bzero;
> +                    r = t << 32;
> +
> +                    t = ptr8[size] << 8 | ptr8[size + 1];
> +                    if (t != header->blank)
> +                        t = t*header->bscale + header->bzero;
> +                    g = t << 16;
> +
> +                    t = ptr8[size * 2] << 8 | ptr8[size * 2 + 1];
> +                    if (t != header->blank)
> +                        t = t*header->bscale + header->bzero;
> +                    b = t;
> +
> +                    *dst64++ = a | r | g | b;
> +                    ptr8 += 2;
> +                }
> +            }
> +        }

These two almost identical blocks need to be factored.

> +    } else {
> +        if (header->bitpix == 8) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst8 = (uint8_t *) (p->data[0] + (avctx->height-i-1)* p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    if (ptr8[0] != header->blank) {
> +                        *dst8++ = ((ptr8[0] - header->data_min) * 255) / (header->data_max - header->data_min);
> +                    } else {
> +                        *dst8++ = 0;
> +                    }
> +                    ptr8++;
> +                }
> +            }
> +        } else if (header->bitpix == 16) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    t16 = AV_RB16(ptr8);
> +                    if (t16 != header->blank) {
> +                        t16 = ((t16 - header->data_min) * 65535) / (header->data_max - header->data_min);
> +                    } else {
> +                        t16 = 0;
> +                    }
> +                    *dst16++ = t16;
> +                    ptr8 += 2;
> +                }
> +            }
> +        } else if (header->bitpix == 32) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    t32 = AV_RB32(ptr8);
> +                    if (t32 != header->blank) {
> +                        t16 = ((t32 - header->data_min) * 65535) / (header->data_max - header->data_min);
> +                    } else {
> +                        t16 = 0;
> +                    }
> +                    *dst16++ = t16;
> +                    ptr8 += 4;
> +                }
> +            }
> +        } else if (header->bitpix == 64) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    t64 = AV_RB64(ptr8);
> +                    if (t64 != header->blank) {
> +                        t16 = ((t64 - header->data_min) * 65535) / (header->data_max - header->data_min);
> +                    } else {
> +                        t16 = 0;
> +                    }
> +                    *dst16++ = t16;
> +                    ptr8 += 8;
> +                }
> +            }

Ditto for these four blocks.

> +        } else if (header->bitpix == -32) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    tflt = av_int2float(AV_RB32(ptr8));
> +                    *dst16++ = ((tflt - header->data_min) * 65535) / (header->data_max - header->data_min);
> +                    ptr8 += 4;
> +                }
> +            }
> +        } else if (header->bitpix == -64) {
> +            for (i = 0; i < avctx->height; i++) {
> +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
> +                for (j = 0; j < avctx->width; j++) {
> +                    tdbl = av_int2double(AV_RB64(ptr8));
> +                    *dst16++ = ((tdbl - header->data_min) * 65535) / (header->data_max - header->data_min);
> +                    ptr8 += 8;
> +                }
> +            }

And ditto for these two, or even possibly these six.

> +        } else {
> +            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header->bitpix);
> +            return AVERROR_INVALIDDATA;
> +        }
> +    }
> +
> +    p->key_frame = 1;
> +    p->pict_type = AV_PICTURE_TYPE_I;
> +
> +    *got_frame = 1;
> +
> +    return avpkt->size;
> +}
> +
> +AVCodec ff_fits_decoder = {
> +    .name           = "fits",
> +    .type           = AVMEDIA_TYPE_VIDEO,
> +    .id             = AV_CODEC_ID_FITS,
> +    .priv_data_size = sizeof(FITSDecContext),
> +    .decode         = fits_decode_frame,
> +    .capabilities   = AV_CODEC_CAP_DR1,
> +    .long_name      = NULL_IF_CONFIG_SMALL("Flexible Image Transport System")
> +};
> diff --git a/libavcodec/version.h b/libavcodec/version.h
> index 0661526..5b99785 100644
> --- a/libavcodec/version.h
> +++ b/libavcodec/version.h
> @@ -28,8 +28,8 @@
>  #include "libavutil/version.h"
> 
>  #define LIBAVCODEC_VERSION_MAJOR  57
> -#define LIBAVCODEC_VERSION_MINOR 100
> -#define LIBAVCODEC_VERSION_MICRO 102
> +#define LIBAVCODEC_VERSION_MINOR 101
> +#define LIBAVCODEC_VERSION_MICRO 100
> 
>  #define LIBAVCODEC_VERSION_INT  AV_VERSION_INT(LIBAVCODEC_VERSION_MAJOR, \
>                                                 LIBAVCODEC_VERSION_MINOR, \
> diff --git a/libavformat/img2.c b/libavformat/img2.c
> index 8432cc0..e405df8 100644
> --- a/libavformat/img2.c
> +++ b/libavformat/img2.c
> @@ -80,6 +80,7 @@ const IdStrMap ff_img_tags[] = {
>      { AV_CODEC_ID_XPM,        "xpm"      },
>      { AV_CODEC_ID_XFACE,      "xface"    },
>      { AV_CODEC_ID_XWD,        "xwd"      },
> +    { AV_CODEC_ID_FITS,       "fits"     },
>      { AV_CODEC_ID_NONE,       NULL       }
>  };
> 

Regards,
Moritz Barsnick June 30, 2017, 7:59 p.m. UTC | #2
On Fri, Jun 30, 2017 at 15:08:17 +0200, Nicolas George wrote:
> > + * function calculates the data_min and data_max values from the data.
> 
> The style for doxygen comments is impersonal verbs: "Calculate the ...".

I don't want to nitpick, and I *may* be wrong, but the grammar freak in
me thinks the suggested form is actually imperative.

Moritz
Paras July 2, 2017, 8:16 a.m. UTC | #3
On Fri, Jun 30, 2017 at 6:38 PM, Nicolas George <george@nsup.org> wrote:

> Hi. A few technical / cosmetic remarks below. I do not know the FITS
> format itself more than in passing.
>
> Le duodi 12 messidor, an CCXXV, Paras Chadha a écrit :
> > Made the changes suggested above
> >
> > Signed-off-by: Paras Chadha <paraschadha18@gmail.com>
> > ---
> >  Changelog               |   1 +
> >  doc/general.texi        |   2 +
> >  libavcodec/Makefile     |   1 +
> >  libavcodec/allcodecs.c  |   1 +
> >  libavcodec/avcodec.h    |   1 +
> >  libavcodec/codec_desc.c |   8 +
> >  libavcodec/fitsdec.c    | 580 ++++++++++++++++++++++++++++++
> ++++++++++++++++++
> >  libavcodec/version.h    |   4 +-
> >  libavformat/img2.c      |   1 +
> >  9 files changed, 597 insertions(+), 2 deletions(-)
> >  create mode 100644 libavcodec/fitsdec.c
> >
> > diff --git a/Changelog b/Changelog
> > index a8726c6..2c2bdec 100644
> > --- a/Changelog
> > +++ b/Changelog
> > @@ -26,6 +26,7 @@ version <next>:
> >    --x86asmexe=yasm to configure to restore the old behavior.
> >  - additional frame format support for Interplay MVE movies
> >  - support for decoding through D3D11VA in ffmpeg
> > +- FITS demuxer and decoder
> >
> >  version 3.3:
> >  - CrystalHD decoder moved to new decode API
> > diff --git a/doc/general.texi b/doc/general.texi
> > index 8f582d5..c00ce32 100644
> > --- a/doc/general.texi
> > +++ b/doc/general.texi
> > @@ -591,6 +591,8 @@ following image formats are supported:
> >      @tab Digital Picture Exchange
> >  @item EXR          @tab   @tab X
> >      @tab OpenEXR
> > +@item FITS         @tab   @tab X
> > +    @tab Flexible Image Transport System
> >  @item JPEG         @tab X @tab X
> >      @tab Progressive JPEG is not supported.
> >  @item JPEG 2000    @tab X @tab X
> > diff --git a/libavcodec/Makefile b/libavcodec/Makefile
> > index b440a00..729e95e 100644
> > --- a/libavcodec/Makefile
> > +++ b/libavcodec/Makefile
> > @@ -291,6 +291,7 @@ OBJS-$(CONFIG_FFV1_DECODER)            += ffv1dec.o
> ffv1.o
> >  OBJS-$(CONFIG_FFV1_ENCODER)            += ffv1enc.o ffv1.o
> >  OBJS-$(CONFIG_FFWAVESYNTH_DECODER)     += ffwavesynth.o
> >  OBJS-$(CONFIG_FIC_DECODER)             += fic.o
> > +OBJS-$(CONFIG_FITS_DECODER)            += fitsdec.o
> >  OBJS-$(CONFIG_FLAC_DECODER)            += flacdec.o flacdata.o flac.o
> >  OBJS-$(CONFIG_FLAC_ENCODER)            += flacenc.o flacdata.o flac.o
> vorbis_data.o
> >  OBJS-$(CONFIG_FLASHSV_DECODER)         += flashsv.o
> > diff --git a/libavcodec/allcodecs.c b/libavcodec/allcodecs.c
> > index 0243f47..a4cfd80 100644
> > --- a/libavcodec/allcodecs.c
> > +++ b/libavcodec/allcodecs.c
> > @@ -192,6 +192,7 @@ static void register_all(void)
> >      REGISTER_ENCDEC (FFV1,              ffv1);
> >      REGISTER_ENCDEC (FFVHUFF,           ffvhuff);
> >      REGISTER_DECODER(FIC,               fic);
> > +    REGISTER_DECODER(FITS,              fits);
> >      REGISTER_ENCDEC (FLASHSV,           flashsv);
> >      REGISTER_ENCDEC (FLASHSV2,          flashsv2);
> >      REGISTER_DECODER(FLIC,              flic);
> > diff --git a/libavcodec/avcodec.h b/libavcodec/avcodec.h
> > index b697afa..8eba460 100644
> > --- a/libavcodec/avcodec.h
> > +++ b/libavcodec/avcodec.h
> > @@ -447,6 +447,7 @@ enum AVCodecID {
> >      AV_CODEC_ID_SRGC,
> >      AV_CODEC_ID_SVG,
> >      AV_CODEC_ID_GDV,
> > +    AV_CODEC_ID_FITS,
> >
> >      /* various PCM "codecs" */
> >      AV_CODEC_ID_FIRST_AUDIO = 0x10000,     ///< A dummy id pointing at
> the start of audio codecs
> > diff --git a/libavcodec/codec_desc.c b/libavcodec/codec_desc.c
> > index cf1246e..0112517 100644
> > --- a/libavcodec/codec_desc.c
> > +++ b/libavcodec/codec_desc.c
> > @@ -1464,6 +1464,14 @@ static const AVCodecDescriptor
> codec_descriptors[] = {
> >                       AV_CODEC_PROP_LOSSLESS,
> >      },
> >      {
> > +        .id        = AV_CODEC_ID_FITS,
> > +        .type      = AVMEDIA_TYPE_VIDEO,
> > +        .name      = "fits",
> > +        .long_name = NULL_IF_CONFIG_SMALL("Flexible Image Transport
> System"),
> > +        .props     = AV_CODEC_PROP_INTRA_ONLY | AV_CODEC_PROP_LOSSY |
> > +                     AV_CODEC_PROP_LOSSLESS,
> > +    },
> > +    {
> >          .id        = AV_CODEC_ID_GIF,
> >          .type      = AVMEDIA_TYPE_VIDEO,
> >          .name      = "gif",
> > diff --git a/libavcodec/fitsdec.c b/libavcodec/fitsdec.c
> > new file mode 100644
> > index 0000000..4eaf3c8
> > --- /dev/null
> > +++ b/libavcodec/fitsdec.c
> > @@ -0,0 +1,580 @@
> > +/*
> > + * FITS image decoder
> > + * Copyright (c) 2017 Paras Chadha
> > + *
> > + * This file is part of FFmpeg.
> > + *
> > + * FFmpeg is free software; you can redistribute it and/or
> > + * modify it under the terms of the GNU Lesser General Public
> > + * License as published by the Free Software Foundation; either
> > + * version 2.1 of the License, or (at your option) any later version.
> > + *
> > + * FFmpeg is distributed in the hope that it will be useful,
> > + * but WITHOUT ANY WARRANTY; without even the implied warranty of
> > + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
> > + * Lesser General Public License for more details.
> > + *
> > + * You should have received a copy of the GNU Lesser General Public
> > + * License along with FFmpeg; if not, write to the Free Software
> > + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
> 02110-1301 USA
> > + */
> > +
> > +/**
>
> > + * @file
> > + * FITS image decoder
> > + * It supports all 2-d images alongwith, bzero, bscale and blank
> keywords.
> > + * RGBA images are supported as NAXIS3 = 3 or 4 i.e. Planes in RGBA
> order. Also CTYPE = 'RGB ' should be present.
> > + * Also to interpret data, values are linearly scaled using min-max
> scaling but not RGB images.
>
> Please add a link to the specification you followed.
>

Done


>
> > + */
> > +
> > +#include "avcodec.h"
> > +#include "internal.h"
> > +#include <float.h>
> > +#include "libavutil/intreadwrite.h"
> > +#include "libavutil/intfloat.h"
> > +#include "libavutil/dict.h"
> > +
> > +/**
> > + * Structure to store the header keywords in FITS file
> > + */
> > +typedef struct FITSContext {
>
> > +    char simple;
>
> This field is only used in fits_read_header(), it should be a local
> variable.
>
> > +    int bitpix;
> > +    int64_t blank;
>
> > +    int naxis;
>
> Ditto.
>
> > +    int naxisn[3];
> > +    int rgb; /**< 1 if file contains RGB image, 0 otherwise */
>
> > +    int xtension;
>
> Ditto.
>
> > +    double bscale;
> > +    double bzero;
> > +    double data_min;
> > +    double data_max;
> > +} FITSDecContext;
>
> In fact, I believe not a single of these fields has a value that
> survives the decode_frame() call: it the whole header structure should
> be a local variable in decode_frame().
>

Yes, i made it a local variable in decode_frame.


>
> > +
> > +/**
>
> > + * function calculates the data_min and data_max values from the data.
>
> The style for doxygen comments is impersonal verbs: "Calculate the ...".
>

ok


>
> > + * This is called if the values are not present in the header.
>
> > + * @param ptr8 - pointer to the data
> > + * @param header - pointer to the header
> > + * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA
>
> Please drop the dashes and comma. The last two comments apply for all
> doxygen comments.
>
> And 0 is usually used for success.
>

ok


>
> > + */
> > +static int fill_data_min_max(const uint8_t * ptr8, FITSDecContext *
> header, const uint8_t * end)
> > +{
> > +    int16_t t16;
> > +    int32_t t32;
> > +    int64_t t64;
> > +    float tflt;
> > +    double tdbl;
> > +    int i, j;
> > +
> > +    header->data_min = DBL_MAX;
> > +    header->data_max = DBL_MIN;
> > +    switch (header->bitpix) {
> > +        case -64:
>
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    tdbl = av_int2double(AV_RB64(ptr8));
> > +                    if (tdbl > header->data_max)
> > +                        header->data_max = tdbl;
> > +                    if (tdbl < header->data_min)
> > +                        header->data_min = tdbl;
> > +                    ptr8 += 8;
> > +                }
> > +            }
> > +            break;
> > +        case -32:
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    tflt = av_int2float(AV_RB32(ptr8));
> > +                    if (tflt > header->data_max)
> > +                        header->data_max = tflt;
> > +                    if (tflt < header->data_min)
> > +                        header->data_min = tflt;
> > +                    ptr8 += 4;
> > +                }
> > +            }
> > +            break;
> > +        case 8:
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    if (ptr8[0] != header->blank) {
> > +                        if (ptr8[0] > header->data_max)
> > +                            header->data_max = ptr8[0];
> > +                        if (ptr8[0] < header->data_min)
> > +                            header->data_min = ptr8[0];
> > +                    }
> > +                    ptr8++;
> > +                }
> > +            }
> > +            break;
> > +        case 16:
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    t16 = AV_RB16(ptr8);
> > +                    if (t16 != header->blank) {
> > +                        if (t16 > header->data_max)
> > +                            header->data_max = t16;
> > +                        if (t16 < header->data_min)
> > +                            header->data_min = t16;
> > +                    }
> > +                    ptr8 += 2;
> > +                }
> > +            }
> > +            break;
> > +        case 32:
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    t32 = AV_RB32(ptr8);
> > +                    if (t32 != header->blank) {
> > +                        if (t32 > header->data_max)
> > +                            header->data_max = t32;
> > +                        if (t32 < header->data_min)
> > +                            header->data_min = t32;
> > +                    }
> > +                    ptr8 += 4;
> > +                }
> > +            }
> > +            break;
> > +        case 64:
> > +            for (i = 0; i < header->naxisn[1]; i++) {
> > +                for (j = 0; j < header->naxisn[0]; j++) {
> > +                    t64 = AV_RB64(ptr8);
> > +                    if (t64 != header->blank) {
> > +                        if (t64 > header->data_max)
> > +                            header->data_max = t64;
> > +                        if (t64 < header->data_min)
> > +                            header->data_min = t64;
> > +                    }
> > +                    ptr8 += 8;
> > +                }
> > +            }
> > +            break;
>
> Please find a way of factoring these almost-identical loops. A macro
> would be the obvious choice.
>

Done


>
> > +        default:
> > +            return AVERROR_INVALIDDATA;
> > +    }
> > +    return 1;
> > +}
> > +
> > +/**
> > + * function reads the fits header and stores the values in
> FITSDecContext pointed by header
> > + * @param avctx - AVCodec context
> > + * @param ptr - pointer to pointer to the data
> > + * @param header - pointer to the FITSDecContext
> > + * @param end - pointer to end of packet
> > + * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA
> > + */
> > +static int fits_read_header(AVCodecContext *avctx, const uint8_t
> **ptr, FITSDecContext * header,
> > +                            const uint8_t * end, AVDictionary **meta)
> > +{
> > +    const uint8_t *ptr8 = *ptr;
> > +    int lines_read = 0, i, dim_no, t, data_min_found = 0,
> data_max_found = 0, ret;
> > +    uint64_t size=1;
> > +    double d;
> > +    AVDictionary *metadata = NULL;
> > +    char keyword[10], value[72];
> > +
>
> > +    header->blank = LLONG_MIN;
>
> This is almost certainly not the value you want.
>

Actually it worked with this value. But, just to be on safe side i have
introduced a new variable blank_found. So, now BLANK pixels are checked on
when BLANK is present in header. It also makes the code more readable.
Also, added an option for the user to enter which value to be replaced for
BLANK pixels in data array.


>
> > +    header->bscale = 1.0;
> > +    header->bzero = 0;
> > +    header->rgb = 0;
> > +
>
> > +    if (end - ptr8 < 80)
> > +        return AVERROR_INVALIDDATA;
> > +
> > +    if (sscanf(ptr8, "SIMPLE = %c", &header->simple) == 1) {
>
> You did not check that there is a NUL character. A crafted input could
> cause an over-read, and relying on the padding is bad style. I suggest
> you check that there is at least one padding space at the end (and make
> sure the sscanf patterns do not). As a macro or a small helper function,
> for example.
>
> Another possibility, more strictly conforming to the standard, would be
> to copy the line into a 81-sized buffer.
>

Okay, i have created a new function which does proper checking while
reading a line. Please look at my new patch. It also makes the code more
modular and readable.


>
> > +        if (header->simple == 'F') {
> > +            av_log(avctx, AV_LOG_WARNING, "not a standard FITS file\n");
> > +            av_dict_set(&metadata, "SIMPLE", "F", 0);
>
> > +        } else if (header->simple != 'T') {
> > +            av_log(avctx, AV_LOG_ERROR, "invalid SIMPLE value, SIMPLE =
> %c\n", header->simple);
> > +            return AVERROR_INVALIDDATA;
> > +        } else {
> > +            av_dict_set(&metadata, "SIMPLE", "T", 0);
> > +        }
>
> I strongly suggest to swap the if and else clause.
>

Please look at the new patch.


>
> > +        header->xtension = 0;
> > +    } else if (!strncmp(ptr8, "XTENSION= 'IMAGE", 16)) {
> > +        header->xtension = 1;
> > +        av_dict_set(&metadata, "XTENSION", "'IMAGE   '", 0);
> > +    } else {
> > +        av_log(avctx, AV_LOG_ERROR, "missing SIMPLE keyword or invalid
> XTENSION\n");
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    ptr8 += 80;
> > +    lines_read++;
> > +
> > +    if (end - ptr8 < 80)
> > +        return AVERROR_INVALIDDATA;
> > +
> > +    if (sscanf(ptr8, "BITPIX = %d", &header->bitpix) != 1) {
> > +        av_log(avctx, AV_LOG_ERROR, "missing BITPIX keyword\n");
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    av_dict_set_int(&metadata, "BITPIX", header->bitpix, 0);
> > +    size = abs(header->bitpix) >> 3;
> > +    ptr8 += 80;
> > +    lines_read++;
> > +
> > +    if (end - ptr8 < 80)
> > +        return AVERROR_INVALIDDATA;
> > +
> > +    if (sscanf(ptr8, "NAXIS = %d", &header->naxis) != 1) {
> > +        av_log(avctx, AV_LOG_ERROR, "missing NAXIS keyword\n");
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    if (!header->naxis) {
> > +        av_log(avctx, AV_LOG_ERROR, "No image data found, NAXIS =
> %d\n", header->naxis);
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    if (header->naxis != 2 && header->naxis != 3) {
> > +        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions,
> NAXIS = %d\n", header->naxis);
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    av_dict_set_int(&metadata, "NAXIS", header->naxis, 0);
> > +    ptr8 += 80;
> > +    lines_read++;
> > +
> > +    for (i = 0; i < header->naxis; i++) {
> > +        if (end - ptr8 < 80)
> > +            return AVERROR_INVALIDDATA;
> > +
> > +        if (sscanf(ptr8, "NAXIS%d = %d", &dim_no, &header->naxisn[i])
> != 2 || dim_no != i+1) {
> > +            av_log(avctx, AV_LOG_ERROR, "missing NAXIS%d keyword\n",
> i+1);
> > +            return AVERROR_INVALIDDATA;
> > +        }
> > +
>
> > +        ret = snprintf(keyword, 10, "NAXIS%d", dim_no);
> > +        if (ret < 0 || ret >= 10) {
>
> Use sizeof(keyword).
>
> > +            return AVERROR_INVALIDDATA;
> > +        }
> > +
> > +        av_dict_set_int(&metadata, keyword, header->naxisn[i], 0);
>
> > +        size *= header->naxisn[i];
>
> This can overflow.
>

Added a check for overflow


>
> > +        ptr8 += 80;
> > +        lines_read++;
> > +    }
> > +
>
> > +    if (end - ptr8 < 80)
> > +            return AVERROR_INVALIDDATA;
>
> Strange indentation.


fixed


>


> > +
> > +    while (strncmp(ptr8, "END", 3)) {
> > +        if (sscanf(ptr8, "BLANK = %d", &t) == 1) {
> > +            header->blank = t;
> > +        } else if (sscanf(ptr8, "BSCALE = %lf", &d) == 1) {
> > +            header->bscale = d;
> > +        } else if (sscanf(ptr8, "BZERO = %lf", &d) == 1) {
> > +            header->bzero = d;
> > +        } else if (sscanf(ptr8, "DATAMAX = %lf", &d) == 1) {
> > +            data_max_found = 1;
> > +            header->data_max = d;
> > +        } else if (sscanf(ptr8, "DATAMIN = %lf", &d) == 1) {
> > +            data_min_found = 1;
> > +            header->data_min = d;
> > +        } else if (!strncmp(ptr8, "CTYPE3  = 'RGB", 14)) {
> > +            header->rgb = 1;
> > +            if (header->naxis != 3 || (header->naxisn[2] != 3 &&
> header->naxisn[2] != 4)) {
> > +                av_log(avctx, AV_LOG_ERROR, "File contains RGB image
> but NAXIS = %d and NAXIS3 = %d\n", header->naxis, header->naxisn[2]);
> > +                return AVERROR_INVALIDDATA;
> > +            }
> > +        }
> > +
> > +        if (ptr8[8] == '=') {
> > +            for (i = 0; i < 8 && ptr8[i] != ' '; i++) {
> > +                keyword[i] = ptr8[i];
> > +            }
> > +            keyword[i] = '\0';
> > +
> > +            t = 0;
> > +            i = 10;
> > +            while (i < 80 && ptr8[i] == ' ')
> > +                i++;
> > +
> > +            if (i < 80) {
> > +                value[t] = ptr8[i];
> > +                i++;
> > +                t++;
> > +                if (ptr8[i-1] == '\'') {
> > +                    while (i < 80 && ptr8[i] != '\'') {
> > +                        value[t] = ptr8[i];
> > +                        i++;
> > +                        t++;
> > +                    }
> > +                    value[t] = '\'';
> > +                    t++;
> > +                } else if (ptr8[i-1] == '(') {
> > +                    while (i < 80 && ptr8[i] != ')') {
> > +                        value[t] = ptr8[i];
> > +                        i++;
> > +                        t++;
> > +                    }
> > +                    value[t] = ')';
> > +                    t++;
> > +                } else {
> > +                    while (i < 80 && ptr8[i] != ' ' && ptr8[i] != '/') {
> > +                        value[t] = ptr8[i];
> > +                        i++;
> > +                        t++;
> > +                    }
> > +                }
> > +            }
> > +
> > +            value[t] = '\0';
> > +            av_dict_set(&metadata, keyword, value, 0);
> > +        }
> > +
> > +        ptr8 += 80;
> > +        lines_read++;
> > +
> > +        if (end - ptr8 < 80)
> > +            return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    if (!header->rgb && header->naxis != 2){
> > +        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions,
> NAXIS = %d\n", header->naxis);
> > +        return AVERROR_INVALIDDATA;
> > +    }
> > +
> > +    ptr8 += 80;
> > +    lines_read++;
> > +    lines_read %= 36;
> > +
> > +    t = ((36 - lines_read) % 36) * 80;
> > +    if (end - ptr8 < t)
> > +        return AVERROR_INVALIDDATA;
> > +    ptr8 += t;
> > +    *ptr = ptr8;
> > +
> > +    if (end - ptr8 < size)
> > +        return AVERROR_INVALIDDATA;
> > +
>
> > +    if (!header->rgb && (!data_min_found || !data_max_found)) {
> > +        if ((ret = fill_data_min_max(ptr8, header, end)) < 0) {
>
> If only one of the min/max is provided in the file, fill_data_min_max()
> will be called and override both. It is on purpose?
>

Yes, it is better to calculate min/max than to use the value in the header.
Calculated value is more reliable than the one present in header. If only
one is present in header, i will have to calculate the other. So, it is
better to find both and override both, isn't it ?


>
> > +            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n",
> header->bitpix);
>
> > +            return AVERROR_INVALIDDATA;
>
> return ret;
>

done


>
> > +        }
> > +    } else {
> > +        /*
> > +         * instead of applying bscale and bzero to every element, we
> can do inverse transformation on data_min and
> > +         * data_max
> > +         */
> > +        header->data_min = (header->data_min - header->bzero) /
> header->bscale;
> > +        header->data_max = (header->data_max - header->bzero) /
> header->bscale;
> > +    }
> > +
> > +    *meta = metadata;
> > +    return 1;
> > +}
> > +
> > +static int fits_decode_frame(AVCodecContext *avctx, void *data, int
> *got_frame, AVPacket *avpkt)
> > +{
> > +    AVFrame *p=data;
> > +    const uint8_t *ptr8 = avpkt->data, *end;
> > +    int16_t t16;
> > +    int32_t t32;
> > +    int64_t t64;
> > +    float   tflt;
> > +    double  tdbl;
> > +    int ret, i, j;
> > +    uint8_t *dst8;
> > +    uint16_t *dst16;
> > +    uint32_t *dst32;
> > +    uint64_t *dst64, size, r, g, b, a, t;
> > +    FITSDecContext * header = avctx->priv_data;
> > +
> > +    end = ptr8 + avpkt->size;
> > +    if ((ret = fits_read_header(avctx, &ptr8, header, end,
> &p->metadata)) < 0)
> > +        return ret;
> > +
> > +    size = (header->naxisn[0]) * (header->naxisn[1]);
> > +
> > +    if (header->rgb) {
> > +        if (header->bitpix == 8) {
> > +            avctx->pix_fmt = AV_PIX_FMT_RGB32;
> > +        } else if (header->bitpix == 16) {
> > +            avctx->pix_fmt = AV_PIX_FMT_RGBA64;
> > +        } else {
> > +            av_log(avctx, AV_LOG_ERROR, "unsupported BITPIX = %d\n",
> header->bitpix);
> > +            return AVERROR_INVALIDDATA;
> > +        }
> > +    } else {
> > +        if (header->bitpix == 8) {
> > +            avctx->pix_fmt = AV_PIX_FMT_GRAY8;
> > +        } else {
> > +            avctx->pix_fmt = AV_PIX_FMT_GRAY16;
> > +        }
> > +    }
> > +
> > +    if ((ret = ff_set_dimensions(avctx, header->naxisn[0],
> header->naxisn[1])) < 0)
> > +        return ret;
> > +
> > +    if ((ret = ff_get_buffer(avctx, p, 0)) < 0)
> > +        return ret;
> > +
> > +    if (header->rgb) {
>
> > +        if (header->bitpix == 8) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                /*
> > +                 * FITS stores images with bottom row first. Therefore
> we have
> > +                 * to fill the image from bottom to top.
> > +                 */
> > +                dst32 = (uint32_t *)(p->data[0] + (avctx->height-i-1)*
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    if (header->naxisn[2] == 4) {
> > +                        if (ptr8[size * 3] != header->blank)
> > +                            t = ptr8[size * 3] * header->bscale +
> header->bzero;
> > +                        a = t << 24;
> > +                    } else {
> > +                        a = (255 << 24);
> > +                    }
> > +
> > +                    if (ptr8[0] != header->blank)
> > +                        t = ptr8[0] * header->bscale + header->bzero;
> > +                    r = t << 16;
> > +
> > +                    if (ptr8[size] != header->blank)
> > +                        t = ptr8[size] * header->bscale + header->bzero;
> > +                    g = t << 8;
> > +
> > +                    if (ptr8[size * 2] != header->blank)
> > +                        t = ptr8[size * 2] * header->bscale +
> header->bzero;
> > +                    b = t;
> > +
> > +                    *dst32++ = ((uint32_t)a) | ((uint32_t)r) |
> ((uint32_t)g) | ((uint32_t)b);
> > +                    ptr8++;
> > +                }
> > +            }
> > +        } else if (header->bitpix == 16) {
> > +            // not tested ....
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst64 = (uint64_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +
> > +                    if (header->naxisn[2] == 4) {
> > +                        t = ((ptr8[size * 3] << 8) | ptr8[size * 3 +
> 1]);
> > +                        if (t != header->blank)
> > +                            t = t*header->bscale + header->bzero;
> > +                        a = t << 48;
> > +                    } else {
> > +                        a = 65535ULL << 48;
> > +                    }
> > +
> > +                    t = ptr8[0] << 8 | ptr8[1];
> > +                    if (t != header->blank)
> > +                        t = t*header->bscale + header->bzero;
> > +                    r = t << 32;
> > +
> > +                    t = ptr8[size] << 8 | ptr8[size + 1];
> > +                    if (t != header->blank)
> > +                        t = t*header->bscale + header->bzero;
> > +                    g = t << 16;
> > +
> > +                    t = ptr8[size * 2] << 8 | ptr8[size * 2 + 1];
> > +                    if (t != header->blank)
> > +                        t = t*header->bscale + header->bzero;
> > +                    b = t;
> > +
> > +                    *dst64++ = a | r | g | b;
> > +                    ptr8 += 2;
> > +                }
> > +            }
> > +        }
>
> These two almost identical blocks need to be factored.
>

done


>
> > +    } else {
> > +        if (header->bitpix == 8) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst8 = (uint8_t *) (p->data[0] + (avctx->height-i-1)*
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    if (ptr8[0] != header->blank) {
> > +                        *dst8++ = ((ptr8[0] - header->data_min) * 255)
> / (header->data_max - header->data_min);
> > +                    } else {
> > +                        *dst8++ = 0;
> > +                    }
> > +                    ptr8++;
> > +                }
> > +            }
> > +        } else if (header->bitpix == 16) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    t16 = AV_RB16(ptr8);
> > +                    if (t16 != header->blank) {
> > +                        t16 = ((t16 - header->data_min) * 65535) /
> (header->data_max - header->data_min);
> > +                    } else {
> > +                        t16 = 0;
> > +                    }
> > +                    *dst16++ = t16;
> > +                    ptr8 += 2;
> > +                }
> > +            }
> > +        } else if (header->bitpix == 32) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    t32 = AV_RB32(ptr8);
> > +                    if (t32 != header->blank) {
> > +                        t16 = ((t32 - header->data_min) * 65535) /
> (header->data_max - header->data_min);
> > +                    } else {
> > +                        t16 = 0;
> > +                    }
> > +                    *dst16++ = t16;
> > +                    ptr8 += 4;
> > +                }
> > +            }
> > +        } else if (header->bitpix == 64) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    t64 = AV_RB64(ptr8);
> > +                    if (t64 != header->blank) {
> > +                        t16 = ((t64 - header->data_min) * 65535) /
> (header->data_max - header->data_min);
> > +                    } else {
> > +                        t16 = 0;
> > +                    }
> > +                    *dst16++ = t16;
> > +                    ptr8 += 8;
> > +                }
> > +            }
>
> Ditto for these four blocks.
>
> > +        } else if (header->bitpix == -32) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    tflt = av_int2float(AV_RB32(ptr8));
> > +                    *dst16++ = ((tflt - header->data_min) * 65535) /
> (header->data_max - header->data_min);
> > +                    ptr8 += 4;
> > +                }
> > +            }
> > +        } else if (header->bitpix == -64) {
> > +            for (i = 0; i < avctx->height; i++) {
> > +                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) *
> p->linesize[0]);
> > +                for (j = 0; j < avctx->width; j++) {
> > +                    tdbl = av_int2double(AV_RB64(ptr8));
> > +                    *dst16++ = ((tdbl - header->data_min) * 65535) /
> (header->data_max - header->data_min);
> > +                    ptr8 += 8;
> > +                }
> > +            }
>
> And ditto for these two, or even possibly these six.
>

yes, done


>
> > +        } else {
> > +            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n",
> header->bitpix);
> > +            return AVERROR_INVALIDDATA;
> > +        }
> > +    }
> > +
> > +    p->key_frame = 1;
> > +    p->pict_type = AV_PICTURE_TYPE_I;
> > +
> > +    *got_frame = 1;
> > +
> > +    return avpkt->size;
> > +}
> > +
> > +AVCodec ff_fits_decoder = {
> > +    .name           = "fits",
> > +    .type           = AVMEDIA_TYPE_VIDEO,
> > +    .id             = AV_CODEC_ID_FITS,
> > +    .priv_data_size = sizeof(FITSDecContext),
> > +    .decode         = fits_decode_frame,
> > +    .capabilities   = AV_CODEC_CAP_DR1,
> > +    .long_name      = NULL_IF_CONFIG_SMALL("Flexible Image Transport
> System")
> > +};
> > diff --git a/libavcodec/version.h b/libavcodec/version.h
> > index 0661526..5b99785 100644
> > --- a/libavcodec/version.h
> > +++ b/libavcodec/version.h
> > @@ -28,8 +28,8 @@
> >  #include "libavutil/version.h"
> >
> >  #define LIBAVCODEC_VERSION_MAJOR  57
> > -#define LIBAVCODEC_VERSION_MINOR 100
> > -#define LIBAVCODEC_VERSION_MICRO 102
> > +#define LIBAVCODEC_VERSION_MINOR 101
> > +#define LIBAVCODEC_VERSION_MICRO 100
> >
> >  #define LIBAVCODEC_VERSION_INT  AV_VERSION_INT(LIBAVCODEC_VERSION_MAJOR,
> \
> >
>  LIBAVCODEC_VERSION_MINOR, \
> > diff --git a/libavformat/img2.c b/libavformat/img2.c
> > index 8432cc0..e405df8 100644
> > --- a/libavformat/img2.c
> > +++ b/libavformat/img2.c
> > @@ -80,6 +80,7 @@ const IdStrMap ff_img_tags[] = {
> >      { AV_CODEC_ID_XPM,        "xpm"      },
> >      { AV_CODEC_ID_XFACE,      "xface"    },
> >      { AV_CODEC_ID_XWD,        "xwd"      },
> > +    { AV_CODEC_ID_FITS,       "fits"     },
> >      { AV_CODEC_ID_NONE,       NULL       }
> >  };
> >
>
> Regards,
>
> --
>   Nicolas George
diff mbox

Patch

diff --git a/Changelog b/Changelog
index a8726c6..2c2bdec 100644
--- a/Changelog
+++ b/Changelog
@@ -26,6 +26,7 @@  version <next>:
   --x86asmexe=yasm to configure to restore the old behavior.
 - additional frame format support for Interplay MVE movies
 - support for decoding through D3D11VA in ffmpeg
+- FITS demuxer and decoder

 version 3.3:
 - CrystalHD decoder moved to new decode API
diff --git a/doc/general.texi b/doc/general.texi
index 8f582d5..c00ce32 100644
--- a/doc/general.texi
+++ b/doc/general.texi
@@ -591,6 +591,8 @@  following image formats are supported:
     @tab Digital Picture Exchange
 @item EXR          @tab   @tab X
     @tab OpenEXR
+@item FITS         @tab   @tab X
+    @tab Flexible Image Transport System
 @item JPEG         @tab X @tab X
     @tab Progressive JPEG is not supported.
 @item JPEG 2000    @tab X @tab X
diff --git a/libavcodec/Makefile b/libavcodec/Makefile
index b440a00..729e95e 100644
--- a/libavcodec/Makefile
+++ b/libavcodec/Makefile
@@ -291,6 +291,7 @@  OBJS-$(CONFIG_FFV1_DECODER)            += ffv1dec.o ffv1.o
 OBJS-$(CONFIG_FFV1_ENCODER)            += ffv1enc.o ffv1.o
 OBJS-$(CONFIG_FFWAVESYNTH_DECODER)     += ffwavesynth.o
 OBJS-$(CONFIG_FIC_DECODER)             += fic.o
+OBJS-$(CONFIG_FITS_DECODER)            += fitsdec.o
 OBJS-$(CONFIG_FLAC_DECODER)            += flacdec.o flacdata.o flac.o
 OBJS-$(CONFIG_FLAC_ENCODER)            += flacenc.o flacdata.o flac.o vorbis_data.o
 OBJS-$(CONFIG_FLASHSV_DECODER)         += flashsv.o
diff --git a/libavcodec/allcodecs.c b/libavcodec/allcodecs.c
index 0243f47..a4cfd80 100644
--- a/libavcodec/allcodecs.c
+++ b/libavcodec/allcodecs.c
@@ -192,6 +192,7 @@  static void register_all(void)
     REGISTER_ENCDEC (FFV1,              ffv1);
     REGISTER_ENCDEC (FFVHUFF,           ffvhuff);
     REGISTER_DECODER(FIC,               fic);
+    REGISTER_DECODER(FITS,              fits);
     REGISTER_ENCDEC (FLASHSV,           flashsv);
     REGISTER_ENCDEC (FLASHSV2,          flashsv2);
     REGISTER_DECODER(FLIC,              flic);
diff --git a/libavcodec/avcodec.h b/libavcodec/avcodec.h
index b697afa..8eba460 100644
--- a/libavcodec/avcodec.h
+++ b/libavcodec/avcodec.h
@@ -447,6 +447,7 @@  enum AVCodecID {
     AV_CODEC_ID_SRGC,
     AV_CODEC_ID_SVG,
     AV_CODEC_ID_GDV,
+    AV_CODEC_ID_FITS,

     /* various PCM "codecs" */
     AV_CODEC_ID_FIRST_AUDIO = 0x10000,     ///< A dummy id pointing at the start of audio codecs
diff --git a/libavcodec/codec_desc.c b/libavcodec/codec_desc.c
index cf1246e..0112517 100644
--- a/libavcodec/codec_desc.c
+++ b/libavcodec/codec_desc.c
@@ -1464,6 +1464,14 @@  static const AVCodecDescriptor codec_descriptors[] = {
                      AV_CODEC_PROP_LOSSLESS,
     },
     {
+        .id        = AV_CODEC_ID_FITS,
+        .type      = AVMEDIA_TYPE_VIDEO,
+        .name      = "fits",
+        .long_name = NULL_IF_CONFIG_SMALL("Flexible Image Transport System"),
+        .props     = AV_CODEC_PROP_INTRA_ONLY | AV_CODEC_PROP_LOSSY |
+                     AV_CODEC_PROP_LOSSLESS,
+    },
+    {
         .id        = AV_CODEC_ID_GIF,
         .type      = AVMEDIA_TYPE_VIDEO,
         .name      = "gif",
diff --git a/libavcodec/fitsdec.c b/libavcodec/fitsdec.c
new file mode 100644
index 0000000..4eaf3c8
--- /dev/null
+++ b/libavcodec/fitsdec.c
@@ -0,0 +1,580 @@ 
+/*
+ * FITS image decoder
+ * Copyright (c) 2017 Paras Chadha
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/**
+ * @file
+ * FITS image decoder
+ * It supports all 2-d images alongwith, bzero, bscale and blank keywords.
+ * RGBA images are supported as NAXIS3 = 3 or 4 i.e. Planes in RGBA order. Also CTYPE = 'RGB ' should be present.
+ * Also to interpret data, values are linearly scaled using min-max scaling but not RGB images.
+ */
+
+#include "avcodec.h"
+#include "internal.h"
+#include <float.h>
+#include "libavutil/intreadwrite.h"
+#include "libavutil/intfloat.h"
+#include "libavutil/dict.h"
+
+/**
+ * Structure to store the header keywords in FITS file
+ */
+typedef struct FITSContext {
+    char simple;
+    int bitpix;
+    int64_t blank;
+    int naxis;
+    int naxisn[3];
+    int rgb; /**< 1 if file contains RGB image, 0 otherwise */
+    int xtension;
+    double bscale;
+    double bzero;
+    double data_min;
+    double data_max;
+} FITSDecContext;
+
+/**
+ * function calculates the data_min and data_max values from the data.
+ * This is called if the values are not present in the header.
+ * @param ptr8 - pointer to the data
+ * @param header - pointer to the header
+ * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA
+ */
+static int fill_data_min_max(const uint8_t * ptr8, FITSDecContext * header, const uint8_t * end)
+{
+    int16_t t16;
+    int32_t t32;
+    int64_t t64;
+    float tflt;
+    double tdbl;
+    int i, j;
+
+    header->data_min = DBL_MAX;
+    header->data_max = DBL_MIN;
+    switch (header->bitpix) {
+        case -64:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    tdbl = av_int2double(AV_RB64(ptr8));
+                    if (tdbl > header->data_max)
+                        header->data_max = tdbl;
+                    if (tdbl < header->data_min)
+                        header->data_min = tdbl;
+                    ptr8 += 8;
+                }
+            }
+            break;
+        case -32:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    tflt = av_int2float(AV_RB32(ptr8));
+                    if (tflt > header->data_max)
+                        header->data_max = tflt;
+                    if (tflt < header->data_min)
+                        header->data_min = tflt;
+                    ptr8 += 4;
+                }
+            }
+            break;
+        case 8:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    if (ptr8[0] != header->blank) {
+                        if (ptr8[0] > header->data_max)
+                            header->data_max = ptr8[0];
+                        if (ptr8[0] < header->data_min)
+                            header->data_min = ptr8[0];
+                    }
+                    ptr8++;
+                }
+            }
+            break;
+        case 16:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    t16 = AV_RB16(ptr8);
+                    if (t16 != header->blank) {
+                        if (t16 > header->data_max)
+                            header->data_max = t16;
+                        if (t16 < header->data_min)
+                            header->data_min = t16;
+                    }
+                    ptr8 += 2;
+                }
+            }
+            break;
+        case 32:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    t32 = AV_RB32(ptr8);
+                    if (t32 != header->blank) {
+                        if (t32 > header->data_max)
+                            header->data_max = t32;
+                        if (t32 < header->data_min)
+                            header->data_min = t32;
+                    }
+                    ptr8 += 4;
+                }
+            }
+            break;
+        case 64:
+            for (i = 0; i < header->naxisn[1]; i++) {
+                for (j = 0; j < header->naxisn[0]; j++) {
+                    t64 = AV_RB64(ptr8);
+                    if (t64 != header->blank) {
+                        if (t64 > header->data_max)
+                            header->data_max = t64;
+                        if (t64 < header->data_min)
+                            header->data_min = t64;
+                    }
+                    ptr8 += 8;
+                }
+            }
+            break;
+        default:
+            return AVERROR_INVALIDDATA;
+    }
+    return 1;
+}
+
+/**
+ * function reads the fits header and stores the values in FITSDecContext pointed by header
+ * @param avctx - AVCodec context
+ * @param ptr - pointer to pointer to the data
+ * @param header - pointer to the FITSDecContext
+ * @param end - pointer to end of packet
+ * @return 1, if calculated successfully, otherwise AVERROR_INVALIDDATA
+ */
+static int fits_read_header(AVCodecContext *avctx, const uint8_t **ptr, FITSDecContext * header,
+                            const uint8_t * end, AVDictionary **meta)
+{
+    const uint8_t *ptr8 = *ptr;
+    int lines_read = 0, i, dim_no, t, data_min_found = 0, data_max_found = 0, ret;
+    uint64_t size=1;
+    double d;
+    AVDictionary *metadata = NULL;
+    char keyword[10], value[72];
+
+    header->blank = LLONG_MIN;
+    header->bscale = 1.0;
+    header->bzero = 0;
+    header->rgb = 0;
+
+    if (end - ptr8 < 80)
+        return AVERROR_INVALIDDATA;
+
+    if (sscanf(ptr8, "SIMPLE = %c", &header->simple) == 1) {
+        if (header->simple == 'F') {
+            av_log(avctx, AV_LOG_WARNING, "not a standard FITS file\n");
+            av_dict_set(&metadata, "SIMPLE", "F", 0);
+        } else if (header->simple != 'T') {
+            av_log(avctx, AV_LOG_ERROR, "invalid SIMPLE value, SIMPLE = %c\n", header->simple);
+            return AVERROR_INVALIDDATA;
+        } else {
+            av_dict_set(&metadata, "SIMPLE", "T", 0);
+        }
+        header->xtension = 0;
+    } else if (!strncmp(ptr8, "XTENSION= 'IMAGE", 16)) {
+        header->xtension = 1;
+        av_dict_set(&metadata, "XTENSION", "'IMAGE   '", 0);
+    } else {
+        av_log(avctx, AV_LOG_ERROR, "missing SIMPLE keyword or invalid XTENSION\n");
+        return AVERROR_INVALIDDATA;
+    }
+
+    ptr8 += 80;
+    lines_read++;
+
+    if (end - ptr8 < 80)
+        return AVERROR_INVALIDDATA;
+
+    if (sscanf(ptr8, "BITPIX = %d", &header->bitpix) != 1) {
+        av_log(avctx, AV_LOG_ERROR, "missing BITPIX keyword\n");
+        return AVERROR_INVALIDDATA;
+    }
+
+    av_dict_set_int(&metadata, "BITPIX", header->bitpix, 0);
+    size = abs(header->bitpix) >> 3;
+    ptr8 += 80;
+    lines_read++;
+
+    if (end - ptr8 < 80)
+        return AVERROR_INVALIDDATA;
+
+    if (sscanf(ptr8, "NAXIS = %d", &header->naxis) != 1) {
+        av_log(avctx, AV_LOG_ERROR, "missing NAXIS keyword\n");
+        return AVERROR_INVALIDDATA;
+    }
+
+    if (!header->naxis) {
+        av_log(avctx, AV_LOG_ERROR, "No image data found, NAXIS = %d\n", header->naxis);
+        return AVERROR_INVALIDDATA;
+    }
+
+    if (header->naxis != 2 && header->naxis != 3) {
+        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions, NAXIS = %d\n", header->naxis);
+        return AVERROR_INVALIDDATA;
+    }
+
+    av_dict_set_int(&metadata, "NAXIS", header->naxis, 0);
+    ptr8 += 80;
+    lines_read++;
+
+    for (i = 0; i < header->naxis; i++) {
+        if (end - ptr8 < 80)
+            return AVERROR_INVALIDDATA;
+
+        if (sscanf(ptr8, "NAXIS%d = %d", &dim_no, &header->naxisn[i]) != 2 || dim_no != i+1) {
+            av_log(avctx, AV_LOG_ERROR, "missing NAXIS%d keyword\n", i+1);
+            return AVERROR_INVALIDDATA;
+        }
+
+        ret = snprintf(keyword, 10, "NAXIS%d", dim_no);
+        if (ret < 0 || ret >= 10) {
+            return AVERROR_INVALIDDATA;
+        }
+
+        av_dict_set_int(&metadata, keyword, header->naxisn[i], 0);
+        size *= header->naxisn[i];
+        ptr8 += 80;
+        lines_read++;
+    }
+
+    if (end - ptr8 < 80)
+            return AVERROR_INVALIDDATA;
+
+    while (strncmp(ptr8, "END", 3)) {
+        if (sscanf(ptr8, "BLANK = %d", &t) == 1) {
+            header->blank = t;
+        } else if (sscanf(ptr8, "BSCALE = %lf", &d) == 1) {
+            header->bscale = d;
+        } else if (sscanf(ptr8, "BZERO = %lf", &d) == 1) {
+            header->bzero = d;
+        } else if (sscanf(ptr8, "DATAMAX = %lf", &d) == 1) {
+            data_max_found = 1;
+            header->data_max = d;
+        } else if (sscanf(ptr8, "DATAMIN = %lf", &d) == 1) {
+            data_min_found = 1;
+            header->data_min = d;
+        } else if (!strncmp(ptr8, "CTYPE3  = 'RGB", 14)) {
+            header->rgb = 1;
+            if (header->naxis != 3 || (header->naxisn[2] != 3 && header->naxisn[2] != 4)) {
+                av_log(avctx, AV_LOG_ERROR, "File contains RGB image but NAXIS = %d and NAXIS3 = %d\n", header->naxis, header->naxisn[2]);
+                return AVERROR_INVALIDDATA;
+            }
+        }
+
+        if (ptr8[8] == '=') {
+            for (i = 0; i < 8 && ptr8[i] != ' '; i++) {
+                keyword[i] = ptr8[i];
+            }
+            keyword[i] = '\0';
+
+            t = 0;
+            i = 10;
+            while (i < 80 && ptr8[i] == ' ')
+                i++;
+
+            if (i < 80) {
+                value[t] = ptr8[i];
+                i++;
+                t++;
+                if (ptr8[i-1] == '\'') {
+                    while (i < 80 && ptr8[i] != '\'') {
+                        value[t] = ptr8[i];
+                        i++;
+                        t++;
+                    }
+                    value[t] = '\'';
+                    t++;
+                } else if (ptr8[i-1] == '(') {
+                    while (i < 80 && ptr8[i] != ')') {
+                        value[t] = ptr8[i];
+                        i++;
+                        t++;
+                    }
+                    value[t] = ')';
+                    t++;
+                } else {
+                    while (i < 80 && ptr8[i] != ' ' && ptr8[i] != '/') {
+                        value[t] = ptr8[i];
+                        i++;
+                        t++;
+                    }
+                }
+            }
+
+            value[t] = '\0';
+            av_dict_set(&metadata, keyword, value, 0);
+        }
+
+        ptr8 += 80;
+        lines_read++;
+
+        if (end - ptr8 < 80)
+            return AVERROR_INVALIDDATA;
+    }
+
+    if (!header->rgb && header->naxis != 2){
+        av_log(avctx, AV_LOG_ERROR, "unsupported number of dimensions, NAXIS = %d\n", header->naxis);
+        return AVERROR_INVALIDDATA;
+    }
+
+    ptr8 += 80;
+    lines_read++;
+    lines_read %= 36;
+
+    t = ((36 - lines_read) % 36) * 80;
+    if (end - ptr8 < t)
+        return AVERROR_INVALIDDATA;
+    ptr8 += t;
+    *ptr = ptr8;
+
+    if (end - ptr8 < size)
+        return AVERROR_INVALIDDATA;
+
+    if (!header->rgb && (!data_min_found || !data_max_found)) {
+        if ((ret = fill_data_min_max(ptr8, header, end)) < 0) {
+            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header->bitpix);
+            return AVERROR_INVALIDDATA;
+        }
+    } else {
+        /*
+         * instead of applying bscale and bzero to every element, we can do inverse transformation on data_min and
+         * data_max
+         */
+        header->data_min = (header->data_min - header->bzero) / header->bscale;
+        header->data_max = (header->data_max - header->bzero) / header->bscale;
+    }
+
+    *meta = metadata;
+    return 1;
+}
+
+static int fits_decode_frame(AVCodecContext *avctx, void *data, int *got_frame, AVPacket *avpkt)
+{
+    AVFrame *p=data;
+    const uint8_t *ptr8 = avpkt->data, *end;
+    int16_t t16;
+    int32_t t32;
+    int64_t t64;
+    float   tflt;
+    double  tdbl;
+    int ret, i, j;
+    uint8_t *dst8;
+    uint16_t *dst16;
+    uint32_t *dst32;
+    uint64_t *dst64, size, r, g, b, a, t;
+    FITSDecContext * header = avctx->priv_data;
+
+    end = ptr8 + avpkt->size;
+    if ((ret = fits_read_header(avctx, &ptr8, header, end, &p->metadata)) < 0)
+        return ret;
+
+    size = (header->naxisn[0]) * (header->naxisn[1]);
+
+    if (header->rgb) {
+        if (header->bitpix == 8) {
+            avctx->pix_fmt = AV_PIX_FMT_RGB32;
+        } else if (header->bitpix == 16) {
+            avctx->pix_fmt = AV_PIX_FMT_RGBA64;
+        } else {
+            av_log(avctx, AV_LOG_ERROR, "unsupported BITPIX = %d\n", header->bitpix);
+            return AVERROR_INVALIDDATA;
+        }
+    } else {
+        if (header->bitpix == 8) {
+            avctx->pix_fmt = AV_PIX_FMT_GRAY8;
+        } else {
+            avctx->pix_fmt = AV_PIX_FMT_GRAY16;
+        }
+    }
+
+    if ((ret = ff_set_dimensions(avctx, header->naxisn[0], header->naxisn[1])) < 0)
+        return ret;
+
+    if ((ret = ff_get_buffer(avctx, p, 0)) < 0)
+        return ret;
+
+    if (header->rgb) {
+        if (header->bitpix == 8) {
+            for (i = 0; i < avctx->height; i++) {
+                /*
+                 * FITS stores images with bottom row first. Therefore we have
+                 * to fill the image from bottom to top.
+                 */
+                dst32 = (uint32_t *)(p->data[0] + (avctx->height-i-1)* p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    if (header->naxisn[2] == 4) {
+                        if (ptr8[size * 3] != header->blank)
+                            t = ptr8[size * 3] * header->bscale + header->bzero;
+                        a = t << 24;
+                    } else {
+                        a = (255 << 24);
+                    }
+
+                    if (ptr8[0] != header->blank)
+                        t = ptr8[0] * header->bscale + header->bzero;
+                    r = t << 16;
+
+                    if (ptr8[size] != header->blank)
+                        t = ptr8[size] * header->bscale + header->bzero;
+                    g = t << 8;
+
+                    if (ptr8[size * 2] != header->blank)
+                        t = ptr8[size * 2] * header->bscale + header->bzero;
+                    b = t;
+
+                    *dst32++ = ((uint32_t)a) | ((uint32_t)r) | ((uint32_t)g) | ((uint32_t)b);
+                    ptr8++;
+                }
+            }
+        } else if (header->bitpix == 16) {
+            // not tested ....
+            for (i = 0; i < avctx->height; i++) {
+                dst64 = (uint64_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+
+                    if (header->naxisn[2] == 4) {
+                        t = ((ptr8[size * 3] << 8) | ptr8[size * 3 + 1]);
+                        if (t != header->blank)
+                            t = t*header->bscale + header->bzero;
+                        a = t << 48;
+                    } else {
+                        a = 65535ULL << 48;
+                    }
+
+                    t = ptr8[0] << 8 | ptr8[1];
+                    if (t != header->blank)
+                        t = t*header->bscale + header->bzero;
+                    r = t << 32;
+
+                    t = ptr8[size] << 8 | ptr8[size + 1];
+                    if (t != header->blank)
+                        t = t*header->bscale + header->bzero;
+                    g = t << 16;
+
+                    t = ptr8[size * 2] << 8 | ptr8[size * 2 + 1];
+                    if (t != header->blank)
+                        t = t*header->bscale + header->bzero;
+                    b = t;
+
+                    *dst64++ = a | r | g | b;
+                    ptr8 += 2;
+                }
+            }
+        }
+    } else {
+        if (header->bitpix == 8) {
+            for (i = 0; i < avctx->height; i++) {
+                dst8 = (uint8_t *) (p->data[0] + (avctx->height-i-1)* p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    if (ptr8[0] != header->blank) {
+                        *dst8++ = ((ptr8[0] - header->data_min) * 255) / (header->data_max - header->data_min);
+                    } else {
+                        *dst8++ = 0;
+                    }
+                    ptr8++;
+                }
+            }
+        } else if (header->bitpix == 16) {
+            for (i = 0; i < avctx->height; i++) {
+                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    t16 = AV_RB16(ptr8);
+                    if (t16 != header->blank) {
+                        t16 = ((t16 - header->data_min) * 65535) / (header->data_max - header->data_min);
+                    } else {
+                        t16 = 0;
+                    }
+                    *dst16++ = t16;
+                    ptr8 += 2;
+                }
+            }
+        } else if (header->bitpix == 32) {
+            for (i = 0; i < avctx->height; i++) {
+                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    t32 = AV_RB32(ptr8);
+                    if (t32 != header->blank) {
+                        t16 = ((t32 - header->data_min) * 65535) / (header->data_max - header->data_min);
+                    } else {
+                        t16 = 0;
+                    }
+                    *dst16++ = t16;
+                    ptr8 += 4;
+                }
+            }
+        } else if (header->bitpix == 64) {
+            for (i = 0; i < avctx->height; i++) {
+                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    t64 = AV_RB64(ptr8);
+                    if (t64 != header->blank) {
+                        t16 = ((t64 - header->data_min) * 65535) / (header->data_max - header->data_min);
+                    } else {
+                        t16 = 0;
+                    }
+                    *dst16++ = t16;
+                    ptr8 += 8;
+                }
+            }
+        } else if (header->bitpix == -32) {
+            for (i = 0; i < avctx->height; i++) {
+                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    tflt = av_int2float(AV_RB32(ptr8));
+                    *dst16++ = ((tflt - header->data_min) * 65535) / (header->data_max - header->data_min);
+                    ptr8 += 4;
+                }
+            }
+        } else if (header->bitpix == -64) {
+            for (i = 0; i < avctx->height; i++) {
+                dst16 = (uint16_t *)(p->data[0] + (avctx->height-i-1) * p->linesize[0]);
+                for (j = 0; j < avctx->width; j++) {
+                    tdbl = av_int2double(AV_RB64(ptr8));
+                    *dst16++ = ((tdbl - header->data_min) * 65535) / (header->data_max - header->data_min);
+                    ptr8 += 8;
+                }
+            }
+        } else {
+            av_log(avctx, AV_LOG_ERROR, "invalid BITPIX, %d\n", header->bitpix);
+            return AVERROR_INVALIDDATA;
+        }
+    }
+
+    p->key_frame = 1;
+    p->pict_type = AV_PICTURE_TYPE_I;
+
+    *got_frame = 1;
+
+    return avpkt->size;
+}
+
+AVCodec ff_fits_decoder = {
+    .name           = "fits",
+    .type           = AVMEDIA_TYPE_VIDEO,
+    .id             = AV_CODEC_ID_FITS,
+    .priv_data_size = sizeof(FITSDecContext),
+    .decode         = fits_decode_frame,
+    .capabilities   = AV_CODEC_CAP_DR1,
+    .long_name      = NULL_IF_CONFIG_SMALL("Flexible Image Transport System")
+};
diff --git a/libavcodec/version.h b/libavcodec/version.h
index 0661526..5b99785 100644
--- a/libavcodec/version.h
+++ b/libavcodec/version.h
@@ -28,8 +28,8 @@ 
 #include "libavutil/version.h"

 #define LIBAVCODEC_VERSION_MAJOR  57
-#define LIBAVCODEC_VERSION_MINOR 100
-#define LIBAVCODEC_VERSION_MICRO 102
+#define LIBAVCODEC_VERSION_MINOR 101
+#define LIBAVCODEC_VERSION_MICRO 100

 #define LIBAVCODEC_VERSION_INT  AV_VERSION_INT(LIBAVCODEC_VERSION_MAJOR, \
                                                LIBAVCODEC_VERSION_MINOR, \
diff --git a/libavformat/img2.c b/libavformat/img2.c
index 8432cc0..e405df8 100644
--- a/libavformat/img2.c
+++ b/libavformat/img2.c
@@ -80,6 +80,7 @@  const IdStrMap ff_img_tags[] = {
     { AV_CODEC_ID_XPM,        "xpm"      },
     { AV_CODEC_ID_XFACE,      "xface"    },
     { AV_CODEC_ID_XWD,        "xwd"      },
+    { AV_CODEC_ID_FITS,       "fits"     },
     { AV_CODEC_ID_NONE,       NULL       }
 };