/* See dgplot.h */
/* Last edited on 2003-10-13 02:05:07 by stolfi */

#include <dgplot.h>
#include <dgtree.h>
#include <dgbezier.h>
#include <dgbasic.h>
#include <js.h>
#include <pswr.h>
#include <pswr_iso.h>
#include <pswr_color.h>
#include <vec.h>
#include <math.h>
#include <values.h>

/* INTERNAL PROTOTYPES */

/* IMPLEMENTATIONS */

void dg_enum_locus_leaves
  ( dg_Dim d,           /* Dimension of grid. */
    dg_Locus E,         /* Root locus. */
    dg_Rank r,          /* Rank of {E}. */
    bz_Patch *b,        /* Shape of {E}, or NULL. */
    dg_NodeStar *NE,    /* Node star of {E}. */
    dg_Rank rMax,       /* Truncate the recursion just below this rank. */
    dg_LocusProc visit  /* Client face-processing procedure. */  
  )
  { dg_Dim s = dg_axes_count(E.norm);  /* Co-dimension of {E}. */
    dg_Dim m = d - s;                  /* Dimension of {E}. */
    int nCells = ipow(2, s);   /* Number of cells in node star {*NE}. */
    bool leaf;
    
    affirm(NE->nel == nCells, "inconsistent node star count");
    affirm(dg_check_cell_rank(E.cell, r), "inconsistent cell rank");
    if (b != NULL) 
      { affirm(b->m == m, "inconsistent shape dimension"); }

    if (r <= rMax)
      { 
        fprintf(stderr, "+ dg_enum_locus_leaves(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r);

        /* Check whether splitting-and-recursion is needed. */ 
        /* The face must be split if it is not flat, or if any adjacent cell
          has a split along an axis parallel to the face: */
        if (r >= rMax)
          { /* A vertex never needs splitting. */
            leaf = TRUE;
          }
        else
          { int i;
            leaf = TRUE; 
            /* Could be improved --- should skip levels with parallel splits. */
            for (i = 0; leaf && (i < nCells); i++)
              { dg_Node *NEi = NE->el[i];
                if (NEi != NULL)
                  { leaf &= (NEi->ch[BLO] == NULL) && (NEi->ch[BHI] == NULL); }
              }
          }

        if (leaf) 
          { /* Leaf item, process it: */
            visit(d, E, r, b, NE);
          }
        else 
          { /* Recurse on children: */
            dg_Axis splax = dg_split_axis(d, r);  /* Split axis. */
            if (dg_axis_belongs(splax,E.norm))
              { /* The split axis is perpendicular to the item. */
                dg_AxisIndex splj = dg_find_axis(E.norm, splax);
                dg_Locus E1 = dg_locus(E.norm,2*E.cell);  /* Same item in next level. */
                /* Shrink node star around element: */
                dg_NodeStar NE1;
                dg_node_star_shrink(NE, splj, &NE1);
                /* Recurse on same item in the next level, with same shape. */
                dg_enum_locus_leaves(d, E1, r+1, b, &NE1, rMax, visit);
                free(NE1.el);
              }
            else
              { /* Axis {splax} is parallel to item {E}, split it: */
                /* Get the spanning axes of {E}: */
                dg_Axes SpnE = dg_axes_complement(d, E.norm);
                /* Find relative position {splj} of {splax} among the {Spn(E)}: */
                dg_AxisIndex splj = dg_find_axis(SpnE, splax);
                
                affirm(splj < d, "inconsistent axis sets");
                fprintf(stderr, "split on axis Spn(E)[%d] = %d\n", splj, splax);

                /* Compute normal axes of splitting item: */
                dg_Axes NrmMD = dg_axis_include(splax, E.norm); 
              
                /* Pieces of this locus: */
                dg_Locus ELO = dg_locus(E.norm, 2*E.cell);
                dg_Locus EHI = dg_locus(E.norm, 2*E.cell + 1);
                dg_Locus EMD = dg_locus(NrmMD, 2*E.cell + 1);
                
                /* Shapes of the pieces: */
                bz_Patch bMD = bz_Patch_facet_new(m, b->n, b->g, splj);
                bz_Patch bLO = bz_Patch_new(m, b->n, b->g);
                bz_Patch bHI = bz_Patch_new(m, b->n, b->g);

                 /* Node stars of the pieces: */
                dg_NodeStar NLO, NHI, NMD;
                
                /* Split the node stars: */
                dg_node_star_split(NE, splj, &NLO, &NHI, &NMD);

                /* Split the item's shape: */
                bz_split(b, splj, 0.5, &bLO, &bHI);
                bz_get_facet(&bHI, splj, BLO, &bMD);
                
                /* Recurse for each half of item {F}: */
                dg_enum_locus_leaves(d, ELO, r+1, &bLO, &NLO, rMax, visit);
                dg_enum_locus_leaves(d, EHI, r+1, &bHI, &NHI, rMax, visit);

                /* Recurse for the separating item {S} of {F}. */
                dg_enum_locus_leaves(d, EMD, r+1, &bMD, &NMD, rMax, visit);

                /* Reclaim temporary storage. */
                bz_Patch_free(bMD); free(NMD.el);
                bz_Patch_free(bLO); free(NLO.el);
                bz_Patch_free(bHI); free(NHI.el);
              }
          }

        fprintf(stderr, "- dg_enum_locus_leaves(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r);
      }
  }

/* POSTSCRIPT PLOTTING OF 2D MULTIGRID ITEMS */

void dg_plot_2D_vertex(PSStream *ps, double *c, bz_RDim n)
  {
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    pswr_set_fill_color(ps, 0.0,0.0,0.0);
    pswr_dot(ps, c[0], c[1], 0.25,  TRUE, FALSE);
  }

void dg_plot_2D_edge(PSStream *ps, double c0[], double c1[], bz_RDim n)
  {
    pswr_set_pen(ps, 0.0, 0.0, 0.0,   0.15, 0.0, 0.0);
    pswr_segment(ps, c0[0], c0[1], c1[0], c1[1]);
  }

#define dg_bilinear(x0,x1,c00,c01,c10,c11) \
  (1-(x0))*((1-(x1))*(c00) + (x1)*(c01)) + (x0)*((1-(x1))*(c10) + (x1)*(c11))
  /* Bilinear interpolation between values {c00,c01,c10,c11} with 
     parameters {x0,x1}, in [0 _ 1]. */

void dg_plot_2D_face
  ( PSStream *ps,
    double c00[], 
    double c01[],
    double c10[], 
    double c11[],
    bz_RDim n,
    int plotDepth,
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *R, double *G, double *B
  )
  {
    int i0, i1, j, pass;
    int steps = ipow(2, plotDepth);
    double p[BZ_MAX_RDIM], q[BZ_MAX_RDIM], r[BZ_MAX_RDIM], s[BZ_MAX_RDIM];
    double fMin = fStart + (kMin+0.5)*fStep;
    double fMax = fStart + (kMax+0.5)*fStep;
    for (pass = 0; pass < 2; pass++)
      {
        for (i0 = 1; i0 <= steps; i0++)
          { for (i1 = 0; i1 <= steps; i1++)
              { double u0 = ((double)i0)/((double)steps);
                double u1 = ((double)i1)/((double)steps);
                double v0 = ((double)i0-1)/((double)steps);
                double v1 = ((double)i1)/((double)steps);
                for (j = 0; j < n; j++)
                  { p[j] = dg_bilinear(u0,u1,c00[j],c01[j],c10[j],c11[j]);
                    r[j] = dg_bilinear(v0,v1,c00[j],c01[j],c10[j],c11[j]);

                  }
                if (i1 > 0) 
                  { if (pass == 0)
                      { if (n == 3)
                          { dg_bands_in_quadrilateral
                              (ps, p, q, r, s, n, fStart, fStep, kMin, kMax, R, G, B);
                          }
                        else
                          { dg_shade_quadrilateral(ps, p, q, r, s, n, fMin, fMax); }
                      }
                    else
                      { dg_isolines_in_quadrilateral
                          (ps, p, q, r, s, n, fStart, fStep, kMin, kMax);
                      }
                  }
                for (j = 0; j < n; j++) { q[j] = p[j]; s[j] = r[j]; }
              }
          }
      }
  }

void dg_shade_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_RDim n,      /* Dimension of range space. */
    double fMin,    /* Minimum function value. */
    double fMax     /* Maximum function value. */
  )
  {
    double ctr[BZ_MAX_RDIM];
    int j;
    for (j = 0; j < n; j++) 
      { ctr[j] = (p00[j]+p01[j]+p10[j]+p11[j])/4.0; }
    /* Paint triangles: */
    dg_shade_triangle(ps, p00, p01, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p01, p11, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p11, p10, ctr, n, fMin, fMax);
    dg_shade_triangle(ps, p10, p00, ctr, n, fMin, fMax);
  }

void dg_shade_triangle
  ( PSStream *ps,
    double a[],     /* First vertex. */
    double b[],     /* Second vertex. */
    double c[],     /* Third vertex. */
    bz_RDim n,      /* Dimension of range space. */
    double fMin,    /* Minimum function value. */
    double fMax     /* Maximum function value. */
  )
  { double R, G, B;
    double f[BZ_MAX_RDIM];
    int j;
    affirm (n >= 2, "invalid point coordinates");
    /* Compute mean function values normalized to [-1 _ +1]: */
    for (j = 0; j < n-2; j++)
      { /* Compute mean function value: */
        f[j] = (a[j+2]+b[j+2]+c[j+2])/3.0;
        /* Clip to {[fMin _ fMax]}: */
        if (f[j] < fMin) { f[j] = fMin; }
        if (f[j] > fMax) { f[j] = fMax; }
        /* Normalize to [-1 _ +1]: */
        if (f[j] > 0.0) 
          { f[j] /= fMax; }
        else if (f[j] < 0.0) 
          { f[j] /= -fMin; }
      }
    /* Compute color {R,G,B} for white-background plotting: */
    double R0 = 1.0, G0 = 0.342, B0 = 0.000;
    if (n == 2)
      { R = 1.0; G = 0.9; B = 0.8; }
    else if (n == 3)
      { pswr_color_scale_1(f[0], R0, G0, B0, &R, &G, &B); }
    else if (n == 4)
      { pswr_color_scale_2(f[0], f[1], R0, G0, B0, &R, &G, &B); }
    else if (n >= 5)
      { pswr_color_scale_3(f[0], f[1], f[2], &R, &G, &B); }
    pswr_set_fill_color(ps, R,G,B);
    pswr_triangle(ps, a[0], a[1], b[0], b[1], c[0], c[1], TRUE, FALSE);
  }

void dg_bands_in_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_RDim n,      /* Dim of point vectors (incl. function values). */
    double fStart,  /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax,       /* Maximum isoline index. */
    double *R, double *G, double *B
  )
  {
    if (n == 3)
      { pswr_bands_in_quadrilateral
          ( ps,
            p00[0], p00[1], p00[2],
            p01[0], p01[1], p01[2],
            p10[0], p10[1], p10[2],
            p11[0], p11[1], p11[2],
            fStart, fStep, kMin, kMax,
            R, G, B
          );
      }
    else
      { int j;
        double f00 = 0.0, f01 = 0.0, f10 = 0.0, f11 = 0.0;
        for (j = 2; j < n; j++)
          { f00 += p00[j]*p00[j]; 
            f01 += p01[j]*p01[j]; 
            f10 += p10[j]*p10[j]; 
            f11 += p11[j]*p11[j];
          }
        pswr_bands_in_quadrilateral
          ( ps, 
            p00[0], p00[1], sqrt(f00), 
            p01[0], p01[1], sqrt(f01), 
            p10[0], p10[1], sqrt(f10), 
            p11[0], p11[1], sqrt(f11), 
            fStart, fStep, kMin, kMax,
            R, G, B
          );
      }
  }

void dg_isolines_in_quadrilateral
  ( PSStream *ps,
    double p00[],   /* Coords and values at corner [0,0]. */
    double p01[],   /* Coords and values at corner [0,1]. */
    double p10[],   /* Coords and values at corner [1,0]. */
    double p11[],   /* Coords and values at corner [1,1]. */
    bz_RDim n,      /* Dimension of range space. */
    double fStart,   /* Synchronize isolines with this level. */
    double fStep,   /* Isoline spacing. */
    int kMin,       /* Minimum isoline index. */
    int kMax        /* Maximum isoline index. */
  )
  {
    if (n == 3)
      { pswr_isolines_in_quadrilateral
          ( ps,
            p00[0], p00[1], p00[2],
            p01[0], p01[1], p01[2],
            p10[0], p10[1], p10[2],
            p11[0], p11[1], p11[2],
            fStart, fStep, kMin, kMax
          );
      }
    else
      { int j;
        double f00 = 0.0, f01 = 0.0, f10 = 0.0, f11 = 0.0;
        for (j = 2; j < n; j++)
          { f00 += p00[j]*p00[j]; 
            f01 += p01[j]*p01[j]; 
            f10 += p10[j]*p10[j]; 
            f11 += p11[j]*p11[j];
          }
        pswr_isolines_in_quadrilateral
          ( ps, 
            p00[0], p00[1], sqrt(f00), 
            p01[0], p01[1], sqrt(f01), 
            p10[0], p10[1], sqrt(f10), 
            p11[0], p11[1], sqrt(f11), 
            fStart, fStep, kMin, kMax
          );
      }
  }