/* An area averaging resize method for PIL. 
 *  Useful for accurate downscaling.  
 *
 * history:
 * 2005-04-25 db  Created
 * 
 */


#include "Imaging.h"

#define PIXEL 4



/* add a row of pixels from the original image to a row
 * of floats that will become pixels in the new image. 
 * This version is for rows that cross a y-axis boundary in
 * the new image -- yfract indicates how much lies on this side
 * of the boundary (0 - 1).   
*/
static inline void
_average_row_fract(int width_in, int width_out, unsigned char *in, 
                   double *row, double step, double yfract){
    int i, x=0;
    int b = PIXEL * (int)step;
    double boundary = step;
    double fract, xfract;
    for (i = 0; i < width_in; i += PIXEL){      
        if (i == b){/* end of pixel  */
            /*finish this pixel */
            xfract = (boundary - floor (boundary));
            fract = yfract *  xfract;
            row[x] += (double)in[i] * fract;
            row[x+1] += (double)in[i+1] * fract;
            row[x+2] += (double)in[i+2] * fract;
            //row[x+3] += ((double)in[i+3]) * fract;

            /* now do the first bit of next one */
            fract = yfract * (1 - xfract);
            x += PIXEL;
            if (x >= width_out){
                /* rounding errors occasionally allow this */
                break;
            }
            row[x] += (double)in[i] * fract;
            row[x+1] += (double)in[i+1] * fract;
            row[x+2] += (double)in[i+2] * fract;
            // row[x+3] += ((double)in[i+3]) * fract;

            /*move the boundary onwards. */
            boundary += step;
            b = ((int)boundary) * PIXEL;            
        }
        else{
            row[x] += ((double)in[i]) * yfract;
            row[x+1] += ((double)in[i+1]) * yfract;
            row[x+2] += ((double)in[i+2]) * yfract;
            //row[x+3] += ((double)in[i+3]) * yfract;
        }
    }
}

/* 
 * for imIn pixels that don't cross an imOut pixel boundary.
 * same as above, but without the yfract (which is 1.0 in this case)
 * XXX you'd think GCC would spot the multiplication by a constant  
 * and optimise it out, but no.   
 */
static inline void
_average_row(int width_in, int width_out, unsigned char *in, 
                   double *row, double step){
    int i, x=0;
    int b = PIXEL * (int)step;
    double boundary = step;
    double fract;
    for (i = 0; i < width_in; i += PIXEL){      
        if (i == b){/* end of pixel  */
            /*finish this pixel */
            fract = (boundary - floor (boundary));
            row[x] += (double)in[i] * fract;
            row[x+1] += (double)in[i+1] * fract;
            row[x+2] += (double)in[i+2] * fract;
            //row[x+3] += ((double)in[i+3]) * fract;

            /* now do the first bit of next one */
            fract = (1 - fract);
            x += PIXEL;
            if (x >= width_out){
                /* rounding errors occasionally allow this */
                break;
            }
            row[x] += (double)in[i] * fract;
            row[x+1] += (double)in[i+1] * fract;
            row[x+2] += (double)in[i+2] * fract;
            //row[x+3] += ((double)in[i+3]) * fract;

            /*move the boundary onwards. */
            boundary += step;
            b = ((int)boundary) * PIXEL;            
        }
        else{
            row[x] += ((double)in[i]);
            row[x+1] += ((double)in[i+1]);
            row[x+2] += ((double)in[i+2]);
            //row[x+3] += ((double)in[i+3]);
        }
    }
}


Imaging
ImagingAreaAverage(Imaging imOut, Imaging imIn)
{
    ImagingSectionCookie cookie;
    int x, y;
    double *row;
    double xstep, ystep;
    xstep = (double) imIn->xsize / imOut->xsize;
    ystep = (double) imIn->ysize / imOut->ysize;

    double scaler = 1.0 / (xstep * ystep);
    int width_in = PIXEL * imIn->xsize;
    int width_out = PIXEL * imOut->xsize;

    if (imOut->xsize < 1 || imOut->ysize < 1){
        /* will segfault trying to scale to 0 width */
        return (Imaging) ImagingError_ValueError("Can't resize to zero width or height!");
    }
    if (xstep < 1.0 || ystep < 1.0){
        /* not downsampling, ignore for now */
        return (Imaging) ImagingError_ValueError("AREA_AVERAGE can only shrink images");
    }
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
        return (Imaging) ImagingError_Mismatch();

    ImagingCopyInfo(imOut, imIn);
    
    int rowsize = width_out * sizeof(double);
    /*  row of imOut where sums accumulate */
    row = (double *) malloc(rowsize);
    if (!row) {
        ImagingDelete(imOut);
        return (Imaging) ImagingError_MemoryError();
    }    

    memset(row, 0, rowsize);
    /* magic thread safety call? */
    ImagingSectionEnter(&cookie);

    double yfract = 0.0;
    double ybound = ystep;
    int yb = (int)ybound;
    int yout = 0;

    for (y = 0; y < imIn->ysize; y++){
        if(y == yb){
            /*new row, and end of old row.*/
            /*first do old row */
            yfract = ybound - floor(ybound);
            if (yfract){
                _average_row_fract(width_in, 
                                   width_out, 
                                   imIn->image[y], 
                                   row, xstep, 
                                   yfract);
            }
            for (x=0; x < width_out; x++){
                imOut->image[yout][x] = (unsigned int)(row[x] * scaler + 0.5);
            }
            memset(row, 0, rowsize);

            yout++;             
            if (yout >= imOut->ysize){
                /* in case it wanders too far, due to rounding */
                break;
            }
            
            ybound += ystep;
            yb = (int)ybound;
            
            /* start new row */
            yfract = 1.0 - yfract;
            if(yfract){
                _average_row_fract(width_in, 
                                   width_out, 
                                   imIn->image[y], 
                                   row, xstep, 
                                   yfract);
            }       
        }
        else {
            /*within row.*/
            _average_row(width_in, 
                               width_out, 
                               imIn->image[y], 
                               row, xstep);//, 1.0);
        }
    }
    /* may sometimes come up short, due to rounding? */
    if (yout < imOut->ysize){
        for ( x= 0; x < width_out; x++){
            imOut->image[yout][x] = (unsigned char)(row[x] * scaler);
        }
    }
    ImagingSectionLeave(&cookie);
    //printf("about to free(row)\n");fflush(stdout);            
    free(row);
    //printf("leaving AreaAverage");fflush(stdout);             
    return imOut;
}