// Last edited on 2013-01-07 11:01:26 by stolfilocal
package minetto.utils;

import minetto.utils.FloatImage; 

/* ----------------------------------------------------------------------
This class provides static functions for resizing a {FloatImage} using Lanczos filters.
Adapted to {FloatImage}s by J.Stolfi (UNICAMP) and R.Minetto (UTFPR) on 2013-01-04. */
public class FloatImageLanczos  {

    /* Filter support: */
    private static double SUPPORT = 3.0;  /* Interpolation filter support. */
    private static int MAX_DEN = 64;      /* Maximum denominator for fractional interpolation. */
    
    /* These attributes are computed for each image scaling job. */
    private int rw;   /* Weight table radius. */
    private int nw;   /* Weight table width (odd {= 2*rw + 1}). */
    private int den;  /* Denominator for fractional interpolation. */
    private double[][] weight;  /* Normalized weights for each phase. */
    /* When the scale factor is 
    
    /* Builds weight tables to rescale a signal from {n_old} samples to {n_new} samples: */
    private FloatImageLanczos(int n_old, int n_new, int shift, int den) {
    
        if (n_old < n_new) {
            /* Up-sampling (interpolation) weights: */
            
        } else if (n_old > n_new) {
            /* Down-sampling (splatting) weights: */
        } else {
            /* Trivial resampling: */
        }
        
        this.rw = (int) ((double) n_old * SUPPORT / (double) n_new);
        assert(this.rw >= 0);
        this.nw = this.rw * 2 + 1;

        /* Build weight tables: */
	this.weight = new double[nw];
	this.normWeight = new double[nw];
	this.tmpWeight = new double[nw];

        double tot_weight = 0;
        for (int i = this.rw; i > 0; i--) {
            double wi = compute_weight(i, n_old, n_new);
            this.weight[this.rw + i] = wi;
            this.weight[this.rw - i] = wi;
            tot_weight += 2*wi;
        }
        this.weight[this.rw] = 1.0;
        tot_weight += 1.0;
        for (int i = 0; i < this.nw; i++) { this.normWeight[i] = this.weight[i] / tot_weight; }
    }
    
    /* ----------------------------------------------------------------------
    Resizes the image {img} to {nx_new} columns and {ny_new} rows. */
    public static FloatImage resize(FloatImage img, int nx_new, int ny_new) {
        FloatImage img_scaled_x = resize_x(img, nx_new);
        FloatImage img_scaled_xy = resize_y(img_scaled_x, ny_new);
        return img_scaled_xy;
    }

    /* ----------------------------------------------------------------------
    Resizes the image {img} to {nx_new} columns, preserving the number of rows. */
    public static FloatImage resize_x(FloatImage img, int nx_new) {
    
        if (nx_new == img.nx) {  return img.copy(); }

        /* Build weight tables: */
        FloatImageLanczos L = new FloatImageLanczos(img.nx, nx_new);
        
        /* Allocate and fill the new image: */
        FloatImage out = new FloatImage(img.nc, nx_new, img.ny);
        for (int x = 0; x < nx_new; x++) {
            
            /* Old column corresponding to new column {x}: */
            int X = (int) (((double) x) * ((double) img.nx) / ((double) nx_new) + 0.5);

            /* Clip the filter window to the old image column range: */
            int startX = X - L.rw;
            int start = 0;
            if (startX < 0) {
                startX = 0;
                start = L.rw - X;
            }

            int stopX = X + L.rw;
            int stop = L.rw * 2;
            if (stopX > (img.nx - 1)) {
                stopX = img.nx - 1;
                stop = L.rw + (img.nx - 1 - X);
            }
            
            /* Compute new samples for all rows and channels in column {x}: */
            double[] pWeight = ((start > 0) || (stop < L.nw - 1) ? L.clip_weight_table(start, stop) : L.normWeight);
            for (int y = 0; y < img.ny; y++) {
                for (int c = 0; c < img.nc; c++) { 
                    float value = apply_x_filter(img, startX, stopX, start, stop, c, y, pWeight);
                    out.set_sample(c, x, y, value);
                }
            }
        }
        return out;
    }

    /* ----------------------------------------------------------------------
    Resizes the image {img} to {ny_new} rows, preserving the number of columns. */
    public static FloatImage resize_y(FloatImage img, int ny_new) {
        
        if (ny_new == img.ny) {  return img.copy(); }
        
	/* Build weight tables: */
        FloatImageLanczos L = new FloatImageLanczos(img.ny, ny_new);
        
        /* Allocate and fill the new image: */
        FloatImage out = new FloatImage(img.nc, img.nx, ny_new);
        for (int y = 0; y < ny_new; y++) {

           /* Old row corresponding to new row {y}: */
           int Y = (int) (((double) y) * ((double) img.ny) / ((double) ny_new) + 0.5);
           
           /* Clip the filter window to the old image row range: */
            int startY = Y - L.rw;
           int start = 0;
           if (startY < 0) {
               startY = 0;
               start = L.rw - Y;
           }

           int stopY = Y + L.rw;
           int stop = L.rw * 2;
           if (stopY > (img.ny - 1)) {
               stopY = img.ny - 1;
               stop = L.rw + (img.ny - 1 - Y);
           }

	   /* Compute new samples for all columns and channels in row {y}: */
           double[] pWeight = ((start > 0) || (stop < L.nw - 1) ? L.clip_weight_table(start, stop) : L.normWeight);
           for (int x = 0; x < img.nx; x++) {
	       for (int c = 0; c < img.nc; c++) { 
                   float value = apply_y_filter(img, startY, stopY, start, stop, c, x, pWeight);
                    out.set_sample(c, x, y, value);
               }
           }
        }
        return out;
    }

    /* ----------------------------------------------------------------------
    Computes the Lanczos filter weight for position {i} in the window, assuming a 
    signal with {n_old} samples is being resized to {n_new} samples: */
    private static double compute_weight(int i, int n_old, int n_new){
        final double PI = (double) 3.14159265358978;
        double x = ((double) i) * ((double) n_new) / ((double) n_old);
        return Math.sin(x * PI) / (x * PI) * Math.sin(x * PI / SUPPORT)/ (x * PI / SUPPORT);
    }

    /* ----------------------------------------------------------------------
    This method of a {FloatImageLanczos} object computes and returns a
    weight table to use when only part of the filter window fits within
    the old image: namely, {weight[start..stop]} instead of
    {weight[0..nw-1]}. */
    private double[] clip_weight_table(int start, int stop) {
        double tot_weight = 0;
        for (int i = start; i <= stop; i++) { tot_weight += weight[i]; }
        for (int i = start; i <= stop; i++) { tmpWeight[i] = weight[i] / tot_weight; }
        return tmpWeight;
    }


    /* ----------------------------------------------------------------------
    Computes one sample in channel {c} and row {y} of the result image when resizing image {img} in the {x} diretion.
    The new sample is obtained by applying the filter weights {weight[start..stop]} to columns {startX..stopX}
    of {img}. */
    private static float apply_x_filter(FloatImage img, int startX, int stopX, int start, int stop, int c, int y, double[] weight) {
        double sum = 0;
        for (int x = startX, j = start; x <= stopX; x++, j++) {
            double value = img.get_sample(c, x, y);
            sum += value*weight[j];
        }
	return (float)sum;
    }

    /* ----------------------------------------------------------------------
    Computes one sample in channel {c} and column {x} of the result image when resizing image {img} in the {y} diretion.
    The new sample is obtained by applying the filter weights {weight[start..stop]} to rows {startY..stopY}
    of {img}. */
    private static float apply_y_filter(FloatImage img, int startY, int stopY, int start, int stop, int c, int x, double[] weight) {
        double sum = 0;
        for (int y = startY, j = start; y <= stopY; y++, j++) {
            double value = img.get_sample(c, x, y);
            sum += value*weight[j];
        }
	return (float)sum;
    }
}