/* Last edited on 2013-03-19 21:12:11 by stolfilocal */
/* Generic parent class for multichannel 2D procedural images with {float}-valued samples. */

package minetto.utils;

import java.lang.Math;
import java.lang.Float;
import java.util.Arrays;

import java.io.File;
import java.io.IOException;

import java.awt.image.BufferedImage;
import java.awt.image.DataBufferInt;
import java.awt.image.DataBufferByte;
import java.awt.image.DataBufferUShort;
import java.awt.image.AffineTransformOp;
import java.awt.geom.AffineTransform;
import javax.imageio.ImageIO; 

public class FloatImageProc {

    /* === REPRESENTATION =========================================== */
    
    /* An instance of this class is a procedural color or monochromatic image. */
    
    public final int nc;   // Number of channels.
        
    /* === BASIC MANIPULATION =========================================== */
    
    /* ----------------------------------------------------------------------
    Evaluates channel {c} at the point with coordinates {(x,y)}.*/
    public float eval_channel (int c, double x, double y);
        
    /* ----------------------------------------------------------------------
    Stores into {val[0..nc]} the image values from all channels at the point with 
    coordinates {(x,y)}.*/
    public void eval (double x, double y, float[] val);
        
    /* ----------------------------------------------------------------------
    Extracts (as a float vector) the pixel in column {ix} and row {iy}. 
    The pixel is assumed to be a unit square {[ix_ix+1]×[iy_iy+1]}.  The sample 
    value {get_sample(c,ix,iy)} is some average of {interpolate(c,x,y)} over that
    square, or in a fuzzy neighborhood of its center {(ix+1/2,iy+1/2)}. from specified column and row. Note the index order. */
    public void get_pixel (int ix, int iy, double[] wx, double[] wy, float[] val) {
        if ((wx == null) && (wy == null)) {
            /* Evaluate just once at pixel center: */
            this.eval(ix + 0.5, iy + 0.5, val);
        } else {
            double dx = 1.0/(1.0 + mx);  /* Subsampling step in {X}. */
            double dy = 1.0/(1.0 + my);  /* Subsampling step in {Y}. */

            double sum[] = new double[this.nc];
            for (int c = 0; c < this.nc; c++) { sum[c] = 0; }
            float smp[] = new float[this.nc];
            for (int kx = -mx; kx <= mx; kx++) {
                double sx = ix + 0.5 + dx*kx
                for (int ky = -my; ky <= my; ky++) {
                    double sy = iy + 0.5 + dy*ky
                    this.eval(sx, sy, smp);
                    for (int c = 0; c < this.nc; c++) { sum[c] += wx[kx+mx]*wy[ky+my]*smp[c]; }
                }
            }
        }
    }
        
    /* ----------------------------------------------------------------------
    Returns a copy of a sub-image specified by index ranges
    {ixmin..ixmax} and {iymin,iymax}. */
    public FloatImage get_sub_image (int ixmin, int iymin, int ixmax, int iymax, int mx, int my) {

        double[] wx = (mx <= 0 ? null : sample_weights(mx));
        double[] wy = (mx == my ? wx : (my <= 0 ? null : sample_weights(my)));
                
        /* Allocate result image: */
        int nc = this.nc;
        int nx = xmax - xmin + 1;
        int ny = ymax - ymin + 1;
        FloatImage res = new FloatImage (nc, nx, ny);

        /* Copy samples: */
        float[] smp = new float[];
        for (int iy = 0; iy < ny; iy++) {
            for (int ix = 0; ix < nx; ix++) { 
                this.get_pixel(ix+xmin, iy+ymin, smp);
                res.set_pixel(ix, iy, smp);
            }
        }
        return res;
    }

    /* ----------------------------------------------------------------------
    Returns a {float[2]} array {[vmin,vmax]} with the minimum and maximum sample
    values in channels {cmin..cmax} of the given image, as returned by {get_sample},
    ignoring invalid samples (infinite values and {NaN}s) and non-existant channels.
    
    Will return an empty range with {vmin > vmax} if there are no valid samples
    in those channels. Will return  a trivial range with {vmin == vmax} if 
    there are some valid samples but they all have the same value. */
    public float[] get_sample_range (int cmin, int cmax) {
        if (cmin < 0) { cmin = 0; }
        if (cmax >= nc) { cmax = nc-1; }
        
        float smax = -Float.MAX_VALUE;
        float smin = +Float.MAX_VALUE;
                
        for (int y = 0; y < this.ny; y++) {
            for (int x = 0; x < this.nx; x++) {   
                for (int c = 0; c < this.nc; c++) {   
                    float val = this.get_sample(nc, x+xmin, y+ymin);
                    if ((! Float.isNaN(val)) && (! Float.isInfinite(val))) {
                        if (val < smin) { smin = val; }
                        if (val > smax) { smax = val; }
                    } 
                }
            }
        }
        return new float[]{ smin, smax };
    }

    /* ----------------------------------------------------------------------
    Returns a {float[nc][2]} with the minimum and maximum values in
    each channel of the given image, ignoring infinite values. */
    public float[][] get_channel_ranges () {
        float[][] ranges = new float[nc][];           
        for (int c = 0; c < nc; c++) { ranges[c] = this.get_sample_range(c, c); } 
        return ranges;
    }

    /* ----------------------------------------------------------------------
    Simple gradient of channel {c} at column {x} and row {y}.
    Currently assumes that {x} and {y} are at least one pixel 
    away from the image border. */
    public double[] gradient_simple (int c, int x, int y) {
        double[] grad = new double[2];
        assert((x >= 1) && (x < nx - 1));
        assert((y >= 1) && (y < ny - 1));
        int pos = (y * nx + x)*nc + c;
        int dx = nc;    /* Position increment for next column. */
        int dy = nx*nc; /* Position increment for next row. */

        double vxm = smp[pos - dx];
        double vxp = smp[pos + dx];
        double vym = smp[pos - dy];
        double vyp = smp[pos + dy];

        grad[0] = (vxp - vxm)/2;
        grad[1] = (vyp - vym)/2;
        return grad;
    }
 
}