/* See {float_image_splat.h}. */
/* Last edited on 2020-10-11 02:48:18 by jstolfi */

#define _GNU_SOURCE
#include <math.h>
#include <limits.h>
#include <assert.h>
#include <string.h>
 
#include <bool.h>
#include <r2.h>
#include <jsmath.h>
#include <affirm.h>
#include <gauss_table.h>
#include <float_image.h>
#include <float_image_func.h>
#include <float_image_splat.h>

#define Pr fprintf
#define Er stderr

void float_image_splat_sample
  ( float_image_t *A, 
    int c, 
    double *vTot,     
    int ix, 
    int iy, 
    float_image_func_t *func, 
    int m
  )
  {
    if ((A == NULL) && (vTot == NULL)) { return; }
    
    int32_t NX, NY;
    if (A == NULL) 
      { /* Compute {vTot} without clipping/folding: */
        NX = 0; NY = 0;
      }
    else
      { NX = (int)A->sz[1];
        NY = (int)A->sz[2];
        /* Skip painting (and leave {vTot} alone) if outside domain: */
        if ((ix < 0) || (ix >= NX)) { return; }
        if ((iy < 0) || (iy >= NY)) { return; }
      }
    
    if (m < 0) { m = 0; }

    /* Get coordinates {xp,yp} of pixel center: */
    double xp = ix + 0.5;
    double yp = iy + 0.5;
        
    /* Sample the new image {func} over the pixel: */
    double sum_w = 0.0;   /* Sum of antialias weights of all subsamples. */
    double sum_wv = 0.0; /* Sum of antialias weights times func value. */
    int dx, dy;
    for (dx = -m; dx <= +m; dx++)
      { double xs = ((double)dx)/(m+1); /* X displacement of subpixel sample. */
        double xa = fabs(xs);
        double wx = (xa < 0.5 ? 1 - 2*xa*xa : 2*(1 - xa)*(1 - xa)); /* X weight factor. */
        for (dy = -m; dy <= m; dy++)
          { double ys = ((double)dy)/(m+1); /* Y displacement of subpixel sample. */
            double ya = fabs(ys);
            double wy = (ya < 0.5 ? 1 - 2*ya*ya : 2*(1 - ya)*(1 - ya)); /* Y weight factor. */
            double w = wx*wy; /* Subsample weight. */
            assert(w > 0.0);
            double v = func(xp + xs, yp + ys);  /* New image value. */
            sum_wv += w*v;
            sum_w += w;
          }
      }
    assert(sum_w > 0);
    
    /* Add to image and/or to {*vTot}: */
    if ((! isnan(sum_wv)) && (sum_wv != 0.0))
      { double v_new = sum_wv/sum_w;
        if (A != NULL)
          { /* Add to the image sample: */
            float *sp = float_image_get_sample_address(A, c, ix, iy);
            double v_old = (*sp);
            (*sp) = (float)(v_old + v_new);
          }
        if(vTot != NULL) { (*vTot) += v_new; }
      }
  }

void float_image_splat_samples
  ( float_image_t *A, 
    int c, 
    double *vTot,     
    int xLo, 
    int xHi,
    int yLo,
    int yHi, 
    float_image_func_t *func,
    int m 
  )
  {
    if ((A == NULL) && (vTot == NULL)) { return; }

    int32_t NX, NY;
    if (A == NULL)
      { /* Update {vTot} without clipping/folding: */
        NX = 0; NY = 0;
      }
    else
      { NX = (int)A->sz[1];
        NY = (int)A->sz[2];
        /* Clip the index ranges to the image domain: */
        if (xLo < 0) { xLo = 0; }
        if (xHi >= NX ) { xHi = NX-1; }
        if (yLo < 0) { yLo = 0; }
        if (yHi >= NY ) { yHi = NY-1; }
      }
    if (m < 0) { m = 0; }
    /* Splat each pixel: */
    int ix, iy;
    for (ix = xLo; ix <= xHi; ix++)
      { for (iy = yLo; iy <= yHi; iy++)
          { float_image_splat_sample(A, c, vTot, ix, iy, func, m); }
      }
  }

void float_image_splat_dot
  ( float_image_t *A, 
    int c, 
    double *vTot,     
    double xctr, 
    double yctr, 
    double rad,
    float v,
    int m
  )
  {
    if (((A == NULL) && (vTot == NULL)) || isnan(v) || (v == 0.0f))
      { /* No splatting to do, leave {vTot} unchanged: */ 
        return;
      }
    
    int32_t NX, NY;
    if (A == NULL)
      { /* Update {vTot} without clipping/folding: */
        NX = 0; NY = 0;
      }
    else
      { /* Clip/fold to image domain: */
        NX = (int)A->sz[1];
        NY = (int)A->sz[2];
      }
    
    /* The sign of {rad} is irrelevant: */
    rad = fabs(rad);
    
    /* Determine the indices of affected columns and rows: */
    int xLo, xHi; float_image_func_get_index_range(xctr, rad, NX, &xLo, &xHi);
    int yLo, yHi; float_image_func_get_index_range(yctr, rad, NY, &yLo, &yHi);
    
    float func_dot(double x, double y)
      { double dx = x - xctr;
        double dy = y - yctr;
        double dr2 = dx*dx + dy*dy;
        if (dr2 <= rad) 
          { return v; }
        else
          { return NAN; }
      }

    float_image_splat_samples(A, c, vTot, xLo, xHi, yLo, yHi, &func_dot, m);
  }
  
void float_image_splat_smudge
  ( float_image_t *A, 
    int c,          /* Channel. */                                
    double *vTot,     
    double xctr,    /* Center's X coordinate. */                  
    double yctr,    /* Center's Y coordinate. */                  
    double xdev,    /* Standard deviation in X direction. */   
    double ydev,    /* Standard deviation in Y direction. */                      
    float v,        /* Value at center. */ 
    int m           /* Subsampling parameter. */
  )
  {
    if (((A == NULL) && (vTot == NULL)) || isnan(v) || (v == 0.0f))
      { /* No splatting to do, leave {vTot} unchanged: */ 
        return;
      }
    
    int32_t NX, NY;
    if (A == NULL)
      { /* Update {vTot} without clipping/folding: */
        NX = 0; NY = 0;
      }
    else
      { /* Clip/fold to image domain: */
        NX = (int)A->sz[1];
        NY = (int)A->sz[2];
      }
    
    /* Determine the affected ranges : */
    double max_devs = gauss_table_BIG_ARG; /* Max devs that is worth considering. */
    int xLo, xHi; float_image_func_get_index_range(xctr, max_devs*xdev, NX, &xLo, &xHi);
    int yLo, yHi; float_image_func_get_index_range(yctr, max_devs*ydev, NY, &yLo, &yHi);
    
    auto float func_bell(double x, double y);
      /* Procedural image with gaussian bell shape and value {v} at center. */

    float func_bell(double x, double y)
      { double gx = gauss_bell_eval(x, xctr, xdev);
        double gy = gauss_bell_eval(y, yctr, ydev);
        return (float)(gx*gy*v);
      }

    float_image_splat_samples(A, c, vTot, xLo, xHi, yLo, yHi, &func_bell, m);
  }