/***************************************************************************
 * This file is a part of CADS/UVS fits2jpeg conversion software           *
 *   Copyright (C) 2012 by CADS/UV Software Team,                          *
 *                         Indian Institute of Astrophysics                *
 *                         Bangalore 560034                                *
 *                         cads_AT_iiap.res.in                             *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program 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 General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

/*Header Definitions*/
#include "fits2jpeg.h"

void scale_image(int scale, int npixels,
                 float *data, JSAMPLE *image_buffer)
{
    unsigned int i = 0;
    int JMAXVAL = 255;
    float datamax = 0.0, datamin = 0.0, tmp = 0.0;
    float hist[256] = {0.0}, cumhist[256] = {0.0};
    float scl_data = 0.0;


    /* first find min & max in data */
    datamax = -1.0 * FLT_MAX;
    datamin = FLT_MAX;
    for (i = 0; i < npixels; ++i)
    {
        if (data[i] > datamax) datamax = data[i];
        if (data[i] < datamin) datamin = data[i];
    } /*endfor*/

    /* Convert data into bytscaled values for jpeg file                 */
    /* the dynamic range is reduced to 255 for jpeg                     */
    scl_data = (datamax - datamin)/(float)JMAXVAL;
    for (i = 0; i < npixels; ++i)
        data[i] = (data[i] - datamin)/scl_data;

    /* All data is now squeezed into the range 0 - 255                   */
    /* NOTE: At this point onwards min & max is 0 and 255 respectively   */
    datamax = (float)JMAXVAL;
    datamin = 0.0;

    /* initialize image histogram. ensure all are zeroes in hist[]       */
    /*-------------------------------------------------------------------*/
    for (i = 0; i <= JMAXVAL; ++i) hist[i] = 0;

    /* construct the image histogram */
    tmp = 1.0/(float)npixels;
    for (i = 0; i <= npixels; ++i)
        hist[(int)floor(data[i])] += tmp;

    /* And the cumulative histogram */
    cumhist[0] = hist[0];
    for (i = 1; i <= JMAXVAL; ++i)
        cumhist[i] += cumhist[i - 1] + hist[i];

    /* Linear scale (min-max) : This is the default scaling
     * if we dont generate image_buffer here, histo-eq will fail */
    for (i = 0; i < npixels; ++i)
        image_buffer[i] = (int)(data[i]);

    /*-------------------------------------------------------------------*/


    switch (scale)
    {
        case 1 :                                              /* Square root */
            printinfo("Using square-root scale");
            scl_data = sqrt((float)JMAXVAL)/(float)JMAXVAL;
            for (i = 0; i < npixels; ++i)
                image_buffer[i] = (int)(sqrt(data[i])/scl_data);
            break;

        case 2 :                                                   /* Square */
            printinfo("Using quadratic scale");
            scl_data = pow((float)JMAXVAL,2)/(float)JMAXVAL;
            for (i = 0; i < npixels; ++i)
                image_buffer[i] = (int)abs((pow(data[i],2) - 1.0)/scl_data);
            break;

        case 3 :                                                    /* Cubic */
            printinfo("Using cubic scale");
            scl_data = pow((float)JMAXVAL,3)/(float)JMAXVAL;
            for (i = 0; i < npixels; ++i)
                image_buffer[i] = (int)abs((pow(data[i],3) - 1.0)/scl_data);
            break;

        case 4 :                                                      /* log */
            printinfo("Using log scale");
            scl_data = log(1.0 + (float)JMAXVAL)/(float)JMAXVAL;
            for (i = 0; i < npixels; ++i)
                image_buffer[i] = (int)((log(abs(data[i]) + 1.0))/scl_data);
            break;

        case 5 :
            /* contrast stretch */
            printinfo("Performing histogram stretch (normalization)");

            /* We need to go through the cumulative histogram to pick the
             *  appropriate values for datamin and datamax               */
            i = 0;
            while (i < JMAXVAL)
            {
                if (cumhist[i] >= 0.01)
                {
                    datamin = (float) i;
                    break;
                }
                i++;
            }
            i = JMAXVAL;
            while (i > 0)
            {
                if (cumhist[i] <= 0.99)
                {
                    datamax = (float) i;
                    break;
                }
                i--;
            }
            scl_data = (datamax - datamin)/(float)JMAXVAL;
            for (i = 0; i < npixels; ++i)
            {
                if (image_buffer[i] >= datamax)
                    image_buffer[i] = JMAXVAL;
                else if (image_buffer[i] <= datamin)
                    image_buffer[i] = 0;
                else
                    image_buffer[i] = (int) abs((image_buffer[i] - datamin)/scl_data);
            }
            break;

        case 6 :
            /* histogram equalization */
            printinfo("Performing Histogram Equalization");
            for (i = 0; i <  npixels; ++i)
                image_buffer[i] = cumhist[image_buffer[i]] * 255;
            break;
        default :
            printinfo("Using linear scale");
            break;
    }

}