// Last edited on 2013-01-05 01:25:27 by stolfilocal
package utils;

import java.awt.image.BufferedImage;
/**
 * 图形处理滤波算法
 * @author  The Jxva Framework Foundation
 * @since   1.0
 * @version 2008-11-27 14:13:19 by Jxva
 */
public class lanczos  {

    private static double support = 3.0;
    private int nHalfDots;   /* Half window width. */
    private int nDots;       /* Window width (odd {= 2*nHalfDots + 1}). */

    /* Weight tables: */
    private double[] contrib;
    private double[] normContrib;
    private double[] tmpContrib;
    
    /* Constructor (private): */
    private lanczos(int n_old, int n_new) {
    
        nHalfDots = (int) ((double) n_old * support / (double) n_new);
        nDots = nHalfDots * 2 + 1;

        /* Build weight tables: */
	contrib = new double[nDots];
	normContrib = new double[nDots];
	tmpContrib = new double[nDots];

        int center = nHalfDots;
        assert(center >= 0);

        contrib[center] = 1.0;

        double weight = 0.0;
        for (int i = 1; i <= center; i++) {
            contrib[center + i] = eval_lanczos(i, n_old, n_new, support);
            weight += contrib[center + i];
        }
        for (int i = center - 1; i >= 0; i--) {
            contrib[i] = contrib[center * 2 - i];
        }

        weight = weight * 2 + 1.0;
        for (int i = 0; i < nDots; i++) {
            normContrib[i] = contrib[i] / weight;
        }
    }
    
    public static FloatImage apply_lanczos (FloatImage img, int nx_new, int ny_new) {
        FloatImage img_scaled_x = HorizontalFiltering(img, nx_new);
        FloatImage img_scaled_xy = VerticalFiltering(img_scaled_x, ny_new);
        return img_scaled_xy;
    }

    /* 图片水平滤波 */
    private static FloatImage HorizontalFiltering(FloatImage img, int nx_new) {
        FloatImage out = new FloatImage(img.nc, nx_new, img.ny);
        
        /* Build weight tables: */
        lanczos L = new lanczos(img.nx, nx_new);
        
        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);

            int startX = X - L.nHalfDots;
            int start = 0;
            if (startX < 0) {
                startX = 0;
                start = L.nHalfDots - X;
            }

            int stopX = X + L.nHalfDots;
            int stop = L.nHalfDots * 2;
            if (stopX > (img.nx - 1)) {
                stopX = img.nx - 1;
                stop = L.nHalfDots + (img.nx - 1 - X);
            }
            
            double[] pContrib = ((start > 0) || (stop < L.nDots - 1) ? L.CalTempContrib(start, stop) : L.normContrib);
            for (int y = 0; y < img.ny; y++) {
                for (int c = 0; c < img.nc; c++) { 
                    float value = HorizontalFilter(img, startX, stopX, start, stop, c, y, pContrib);
                    out.set_sample(c, x, y, value);
                }
            }
        }
        return out;
    }

    /* 图片垂直滤波 */
    private static FloatImage VerticalFiltering(FloatImage img, int ny_new) {
        FloatImage out = new FloatImage(img.nc, img.nx, ny_new);

	/* Build weight tables: */
        lanczos L = new lanczos(img.ny, 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);
           
           int startY = Y - L.nHalfDots;
           int start = 0;
           if (startY < 0) {
               startY = 0;
               start = L.nHalfDots - Y;
           }

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

	   double[] pContrib = ((start > 0) || (stop < L.nDots - 1) ? L.CalTempContrib(start, stop) : L.normContrib);
           for (int x = 0; x < img.nx; x++) {
	       for (int c = 0; c < img.nc; c++) { 
                   float value = VerticalFilter(img, startY, stopY, start, stop, c, x, pContrib);
                    out.set_sample(c, x, y, value);
               }
           }
        }
        return out;
    }

    private static double eval_lanczos(int i, int inWidth, int outWidth, double support){
        final double PI = (double) 3.14159265358978;
        double x = (double) i * (double) outWidth / (double) inWidth;
        return Math.sin(x * PI) / (x * PI) * Math.sin(x * PI / support)/ (x * PI / support);
    }

    /* 处理边缘 */
    private double[] CalTempContrib(int start, int stop) {
        double weight = 0;
        int i = 0;
        for (i = start; i <= stop; i++) {
            weight += contrib[i];
        }
        for (i = start; i <= stop; i++) {
            tmpContrib[i] = contrib[i] / weight;
        }
        return tmpContrib;
    }


	/* 行水平滤波 */
    private static float HorizontalFilter(FloatImage img, int startX, int stopX, int start, int stop, int c, int y, double[] pContrib) {
        double sum = 0;
        for (int x = startX, j = start; x <= stopX; x++, j++) {
            double value = img.get_sample(c, x, y);
            sum += value*pContrib[j];
        }
	return (float)sum;
    }

    /* 行垂直滤波 */
    private static float VerticalFilter(FloatImage img, int startY, int stopY, int start, int stop, int c, int x, double[] pContrib) {
        double sum = 0;
        for (int y = startY, j = start; y <= stopY; y++, j++) {
            double value = img.get_sample(c, x, y);
            sum += value*pContrib[j];
        }
	return (float)sum;
    }
}