/* Last edited on 2013-03-21 20:43:27 by stolfilocal */
/* Faints geometric figures in {FloatImage}. */

package minetto.utils;

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

import minetto.utils.FloatImage; 

public class FloatImagePaint {

    private static boolean debug = false;

    /* ----------------------------------------------------------------------
    Paints an axis-aligned rectangle into {img}, spanning coordinate ranges
    {[xmin_xmax]} and {ymin_ymax]}. If {thck} is positive and less than
    {min(xmax-xmin,ymax-ymin)/2}, paints only a frame {thck} wide
    inside the rectangle. */
    public static void paint_rectangle (
        FloatImage img, 
        double xmin,  
        double xmax, 
        double ymin, 
        double ymax, 
        double thck, 
        float[] rgb
    ) {
        if (debug) {
            System.err.println("    ractangle X = [" + String.format("%6.1f",xmin) + " _ " + String.format("%6.1f",xmax) + "]");
            System.err.println("    rectangle Y = [" + String.format("%6.1f",ymin) + " _ " + String.format("%6.1f",ymax) + "]");
        }
        
        /* We model a pixel as a fuzzy square with radius {rp} centered at half-integer coords. */
        double rp = 1.0;
        
        int[] xrange = pixel_index_range(xmin, xmax, rp);
        int[] yrange = pixel_index_range(ymin, ymax, rp);
        double xctr = (xmax + xmin)/2;
        double yctr = (ymax + ymin)/2;
               
        for(int y = yrange[0]; y <= yrange[1]; y++) {
            for(int x = xrange[0]; x <= xrange[1]; x++) {
                if ((x >= 0) && (x < img.nx) && (y >= 0) && (y < img.ny)) {
                    /* Coords of pixel center relative to rectangle axes: */
                    double xp = x + 0.5;
                    double yp = y + 0.5;
                    /* Compute coverage fraction of pixel by rectangle: */
                    double cov_ex = 
                        interval_interval_coverage(xp, rp, xmin, xmax) *
                        interval_interval_coverage(yp, rp, ymin, ymax); 
                    double cov_in = 0.0;
                    if ((thck > 0) && (thck < Math.min(xmax-xmin,ymax-ymin)/2)) {
                        cov_in = 
                            interval_interval_coverage(xp, rp, xmin+thck, xmax-thck) *
                            interval_interval_coverage(yp, rp, ymin+thck, ymax-thck); 
                    }
                    double cov = Math.min(1, Math.max(0, cov_ex - cov_in));
                    if (cov > 1.0e-6) {
                        float[] pix = img.get_pixel(x, y);
                        for (int c = 0; c < img.nc; c++) { pix[c] = (float)((1-cov)*pix[c] + cov*rgb[c]); }
                        img.set_pixel(x, y, pix);
                    }
                }
            }
        }
    }
    
    /* ----------------------------------------------------------------------
    Fraction of the fuzzy X-interval {[xp-rp _ xp+rp]} that is covered by {[xmin _ xmax]}. */
    public static double interval_interval_coverage(double xp, double rp, double xmin, double xmax) {
        assert(xmin <= xmax);
        assert(rp >= 0);
        if ((xp + rp <= xmin) || (xp - rp >= xmax)) {
            return 0;
        } else if ((xp - rp >= xmin) && (xp + rp <= xmax)) {
            return 1;
        } else {
            double clo = sigmoid((xmin - xp)/rp);
            double chi = sigmoid((xmax - xp)/rp);
            return chi - clo;
        }
    }
    
    /* ----------------------------------------------------------------------
    A sigmoid that rises smoothly from 0 to 1 when {z} ranges in {[-1 _ +1]}.
    It is the integral of the Hann window in {[-1 _ +1]}. */
    public static double sigmoid(double z) {
        if (z <= -1) {
            return 0;
        } else if (z > +1) {
            return 1;
        } else {
            return (1 + z + Math.sin(Math.PI*z)/Math.PI)/2;
        }
    }
    
    /* ----------------------------------------------------------------------
    Paints an axis-aligned ellipse into {img}, given the center
    {xctr,yctr} and the main radii {xrad,yrad}. If {thck} is positive
    and less than {min(xrad,yrad)}, excludes the inner ellipse with main
    radii {xrad-thck,yrad-thck}. */
    public static void paint_axial_ellipse (
        FloatImage img, 
        double xctr,  
        double yctr, 
        double xrad, 
        double yrad, 
        double thck, 
        float[] rgb
    ) {
        if (debug) {
            System.err.println("    ellipse ctr   = " + String.format("%6.1f",xctr) + " " + String.format("%6.1f",yctr));
            System.err.println("    ellipse radii = " + String.format("%6.1f",xrad) + " " + String.format("%6.1f",yrad));
        }
        
        /* We model a pixel as a fuzzy disk with radius {rp} centered at half-integer coords. */
        double rp = 1.0;
                
        int[] xrange = pixel_index_range(xctr - xrad, xctr + xrad, rp);
        int[] yrange = pixel_index_range(yctr - yrad, yctr + yrad, rp);
               
        for(int y = yrange[0]; y <= yrange[1]; y++) {
            for(int x = xrange[0]; x <= xrange[1]; x++) {
                if ((x >= 0) && (x < img.nx) && (y >= 0) && (y < img.ny)) {
                    /* Coords and distance of pixel center relative to ellipse center: */
                    double xp = x + 0.5 - xctr;
                    double yp = y + 0.5 - yctr;
                    double dp = Math.hypot(xp, yp);
                    /* Find radii of inner and outer ellipse in direction {(xp,yp)}: */
                    double re_ex = ellipse_radius(xrad, yrad, xp, yp);
                    double re_in = (thck <= 0 ? 0.0 : ellipse_radius(xrad-thck, yrad-thck, xp, yp));
                    /* Compute coverage fraction of pixel by ellipse: */
                    double cov_ex = ellipse_pixel_coverage(dp, rp, re_ex); 
                    double cov_in = (thck <= 0 ? 0.0 : ellipse_pixel_coverage(dp, rp, re_in)); 
                    double cov = Math.min(1, Math.max(0, cov_ex - cov_in));
                    if (cov > 1.0e-6) {
                        float[] pix = img.get_pixel(x, y);
                        for (int c = 0; c < img.nc; c++) { pix[c] = (float)((1-cov)*pix[c] + cov*rgb[c]); }
                        img.set_pixel(x, y, pix);
                    }
                }
            }
        }
    }
    
    /* ----------------------------------------------------------------------
    Paints a dot onto {img}, given the center {xctr,yctr} radius {rad}.
    If {thck} is positive and less than {rad}, excludes the inner circle
    with radius {rad-thck}. */
    public static void paint_dot (
        FloatImage img, 
        double xctr,  
        double yctr, 
        double rad, 
        double thck, 
        float[] rgb
    ) {
        if (debug) {
            System.err.println("    dot ctr   = " + String.format("%6.1f",xctr) + " " + String.format("%6.1f",yctr));
            System.err.println("    dot radius = " + String.format("%6.1f",rad));
        }
        
        /* We model a pixel as a fuzzy disk with radius {rp} centered at half-integer coords. */
        double rp = 1.0;
                
        int[] xrange = pixel_index_range(xctr - rad, xctr + rad, rp);
        int[] yrange = pixel_index_range(yctr - rad, yctr + rad, rp);
               
        for(int y = yrange[0]; y <= yrange[1]; y++) {
            for(int x = xrange[0]; x <= xrange[1]; x++) {
                if ((x >= 0) && (x < img.nx) && (y >= 0) && (y < img.ny)) {
                    /* Coords and distance of pixel center relative to dot center: */
                    double xp = x + 0.5 - xctr;
                    double yp = y + 0.5 - yctr;
                    double dp = Math.hypot(xp, yp);
                    /* Find radii of inner and outer circle in direction {(xp,yp)}: */
                    double rc_ex = rad;
                    double rc_in = (thck <= 0 ? 0.0 : rad-thck);
                    /* Compute coverage fraction of pixel by dot: */
                    double cov_ex = circle_pixel_coverage(dp, rp, rc_ex); 
                    double cov_in = (thck <= 0 ? 0.0 : circle_pixel_coverage(dp, rp, rc_in)); 
                    double cov = Math.min(1, Math.max(0, cov_ex - cov_in));
                    if (cov > 1.0e-6) {
                        float[] pix = img.get_pixel(x, y);
                        for (int c = 0; c < img.nc; c++) { pix[c] = (float)((1-cov)*pix[c] + cov*rgb[c]); }
                        img.set_pixel(x, y, pix);
                    }
                }
            }
        }
    }
    
    /* ----------------------------------------------------------------------
    Assumes that a pixel with column index {ix} is a fuzzy subset of 
    the plane that is zero at any coordinate {x} such that {abs(ix+0.5-x) >= rp}
    Returns a pair of integers {[ixmin,ixmax]]} such that 
    only the pixels with column index {ix} in {ixmin .. ixmax} 
    will intersect the open interval {(xmin _ xmax)}. */
    public static int[] pixel_index_range(double xmin, double xmax, double rp) {
        int ixmin = (int)Math.ceil(xmin - 0.5 - rp) + 1;
        int ixmax = (int)Math.floor(xmax - 0.5 + rp) - 1;
        return new int[]{ ixmin, ixmax};
    }
    
    
    /* ----------------------------------------------------------------------
    Returns the distance {re} from ellipse center
    to ellipse's edge in the direction of the vector {xp.yp}. */
    private static double ellipse_radius(double xrad, double yrad, double xp, double yp) {
        double xm = xp/xrad;
        double ym = yp/yrad;
        double dm = Math.hypot(xm,ym);
        if (dm == 0) { return Math.min(xrad, yrad); }
        return Math.hypot(xrad*xm/dm, yrad*ym/dm);
    }
    
    /* ----------------------------------------------------------------------
    Fraction of pixel that is covered by ellipse, given distance {d} from ellipse center to pixel center,
    the nominal radius {rp} of the pixel, and distance {re} from ellipse center
    to ellipse's edge in the direction of the pixel center. */
    private static double ellipse_pixel_coverage(double d, double rp, double re) {
        assert(2*re > rp); /* Too complicated if ellipse is too small. */
        assert(d >= 0);
        if (re < d - rp) {
            return 0.0; 
        } else if (re > d + rp) { 
            return 1.0;
        } else {
            double z = (re - d)/rp;
            double cov = sigmoid(z);
            return cov;
        }
    }
    
    /* ----------------------------------------------------------------------
    Fraction of pixel that is covered by circle, given distance {d} from circle center to pixel center,
    the nominal radius {rp} of the pixel, and the nominal radius {rc} of the circle. */
    private static double circle_pixel_coverage(double d, double rp, double rc) {
        assert(d >= 0);
        assert(rc >= 0);
        assert(rp >= 0);
        if (rc + rp < d) {
            /* Disjoint */
            return 0.0; 
        } else if (rc > d + rp) { 
            /* Pixel is inside circle */
            return 1.0;
        } else if (rp > d + rc) { 
            /* Circle is inside pixel */
            return (rc/rp)*(rc/rp);
        } else {
            /* Circles intersect but neither is contained in the other: */
            double covmax = (rc >= rp ? 1.0 : (rc/rp)*(rc/rp));
            double dmax = rc + rp;
            double dmin = (rc >= rp ? rc - rp : rp - rc);
            double dmed = (dmax + dmin)/2;
            double drad = (dmax - dmin)/2;
            double z = (dmed - d)/drad;
            double cov = covmax*sigmoid(z);
            return cov;
        }
    }
                           
}