// Last edited on 2013-03-19 18:04:43 by stolfilocal

/* ----------------------------------------------------------------------
Consider reserving the last channel for an alpha channel, then these images
can be composited over other images. */

    /* ----------------------------------------------------------------------
    Returns an image with {nc} channels, {nx} columns, and {ny} rows,
    containing an image of an ellipse on a background of all zeros.*/
    public static FloatImage ellipse( 
        int nx,
        int ny,
        int c,          /* Channels. */
        double[] ctr;   /* Center. */
        double[] u;     /* Major radius vector. */
        double[] v;     /* Major radius vector. */
        double hwd,     /* Radius of pen tip. */
        float[] vfill,  /* Ink value for the interior, or {null}. */  
        float[] vdraw,  /* Ink value for stroking the outline, of {null}. */  
        int m           /* Subsampling parameter. */
    ) {
        FloatImage res = new FloatImage(nc,nx,ny);
        res.set_all_samples(0f);
        
        /* If the drawing ink is invisible, ignore the line width: */
        if (vdraw = null) { hwd = 0.0; }

        /* If no fill and no draw, there is nothing to do: */
        if ((vfill = null) && (hwd <= 0)) { /* Nothing to do: */ return res; }

        /* Get bounding box with some skosh: */
        int[][] bbox = ellipse_int_bbox(ctr, u, v, hwd + 2.0, nx, ny);

        auto float[] func_ellip(double x, double y);
          /* Returns {vdraw} if point {x,y} of {A} is 
            close to the border, else {vfill} if inside the ellipse,
            else {NAN}. */

        float_image_paint_samples(A, c, bbox[0][0], bbox[0][1], bbox[1][0], bbox[1][1], &func_ellip, NULL, m);

        /* IMPLEMENTATION OF INTERNAL PROCS */

        float func_ellip(double x, double y)
          { 
            /* Compute point coordinates {p} relative to center: */
            r2_t p = (r2_t){{ x - E->ctr.c[0], y - E->ctr.c[1] }};
            float ink;
            if (hwd <= 0)
              { /* No outline; just check inside/outside: */
                ink = (ellipse_ouv_inside(&F, &p) ? vfill : NAN);
              }
            else
              { /* Has outline; check whether we are close to the border: */
                double f = ellipse_ouv_border_position(&F, hwd, &p);
                if (f >= +1.0)
                  { ink = NAN; }
                else if (f <= -1.0)
                  { ink = vfill; }
                else
                  { ink = vdraw; }
              }
            return ink;
          }
      }

    
void float_image_paint_samples
  ( float_image_t *A, 
    int c, 
    int xLo, 
    int xHi,
    int yLo,
    int yHi, 
    float_image_func_t *func, 
    float_image_func_t *mask,
    int m 
  )
  {
    int nx = A->sz[1];
    int ny = A->sz[2];
    /* Clip the arguments to the valid range: */
    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; }
    /* Paint each pixel: */
    int ix, iy;
    for (ix = xLo; ix <= xHi; ix++)
      { for (iy = yLo; iy <=yHi; iy++)
          { float_image_paint_sample(A, c, ix, iy, func, mask, m); }
      }
  }

void float_image_paint_sample
  ( float_image_t *A, 
    int c, 
    int ix, 
    int iy, 
    float_image_func_t *func, 
    float_image_func_t *mask, 
    int m
  )
  {
    int nx = A->sz[1];
    int ny = A->sz[2];
    
    /* Clip the arguments to the valid range: */
    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_wp = 0.0;  /* Sum of antialias weights times mask values. */
    double sum_wpv = 0.0; /* Sum of antialias weights times mask values times paint 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. */
            /* Get the overlay's opacity {p} at the sample point: */
            double p = (isnan(v) ? 0.0 : (mask == NULL ? 1.0 : mask(xp + xs, yp + ys)));
            assert(! isnan(p));
            p = fmax(0.0, fmin(1.0, p));
            if (p != 0.0)
              { sum_wp += w*p;
                sum_wpv += w*p*v;
              }
            sum_w += w;
          }
      }
    
    if (sum_wp > 0)
      { /* Overlay is not totally transparent, mix it into the image sample: */
        assert(sum_w > 0.0);
        double w_old = 1.0 - sum_wp/sum_w;   
        float *sp = float_image_get_sample_address(A, c, ix, iy);
        double v_old = (*sp);
        (*sp) = w_old * v_old + sum_wpv/sum_w;
      }
  }