/* See dg_plot.h */
/* Last edited on 2009-08-24 13:24:07 by stolfi */

#include <dg_plot.h>
#include <mdg_grid.h>
#include <dg_locus.h>
#include <dg_tree.h>

#include <bz_basic.h>
#include <affirm.h>
#include <set32.h>
#include <jsmath.h>
#include <pswr.h>
#include <pswr_iso.h>
#include <pswr_color.h>
#include <vec.h>

#include <math.h>
#include <assert.h>
#include <values.h>

/* INTERNAL PROTOTYPES */

/* IMPLEMENTATIONS */

void dg_enum_locus_leaves
  ( mdg_dim_t d,            /* Dimension of grid. */
    dg_locus_t E,           /* Root locus. */
    mdg_rank_t r,           /* Rank of {E}. */
    bz_patch_t *b,          /* Shape of {E}, or NULL. */
    mdg_grid_size_t psz[],  /* Size vector of the tree pack {NE}. */
    dg_tree_vec_t *NE,      /* Tree star of {E}. */
    mdg_rank_t rMax,        /* Truncate the recursion just below this rank. */
    dg_locus_proc_t visit   /* Client face-processing procedure. */  
  )
  { 
    int nCells = NE->ne; /* Number of cells in star. */
    mdg_dim_t s = set32_count(E.norm);  /* Co-dimension of {E}. */
    mdg_dim_t m = d - s;                  /* Dimension of {E}. */
    bool_t leaf;
    
    /* Paranoia: */
    dg_tree_pack_check(NE, d, psz, mdg_cell_rank(E.cell));
    
    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)
          { /* Pretend it is a leaf locus. */
            leaf = TRUE;
          }
        else
          { /* Check whether it is a leaf locus: */
            int i;
            /* dg_tree_pack_check(NE, d, psz, mdg_cell_rank(E.cell)); */
            leaf = TRUE; 
            /* Could be improved --- should skip levels with parallel splits. */
            for (i = 0; leaf && (i < nCells); i++)
              { dg_tree_node_t *NEi = NE->e[i];
                if (NEi != NULL)
                  { leaf &= (LOCH(*NEi) == NULL) || (HICH(*NEi) == NULL); }
              }
          }
        if (leaf) 
          { /* Leaf item, process it: */
            visit(d, E, r, b, NE);
          }
        else 
          { /* Recurse on children: */
            mdg_axis_t splax = mdg_split_axis(d, r);  /* Split axis. */
            if (set32_belongs(splax,E.norm))
              { /* The split axis is perpendicular to the item. */
                dg_locus_t E1 = dg_locus(E.norm,2*E.cell);  /* Same item in next level. */
                /* Shrink tree star around element: */
                dg_tree_vec_t NE1;
                assert(psz[splax] == 2);
                dg_tree_pack_child(d, psz, NE, splax, 1, &NE1);
                /* Recurse on same item in the next level, with same shape. */
                dg_enum_locus_leaves(d, E1, r+1, b, psz, &NE1, rMax, visit);
                free(NE1.e);
              }
            else
              { /* Axis {splax} is parallel to item {E}, split it: */
                /* Get the spanning axes of {E}: */
                mdg_axis_set_t SpnE = mdg_axis_set_complement(d, E.norm);
                /* Find relative position {splj} of {splax} among the {Spn(E)}: */
                mdg_axis_index_t splj = set32_index_from_elem(splax, SpnE);
                fprintf(stderr, "split on axis Spn(E)[%d] = %d\n", splj, splax);
                affirm((splj >= 0) &&(splj < d), "inconsistent axis sets");
                assert(psz[splax] == 1);

                /* Compute normal axes of splitting item: */
                mdg_axis_set_t NrmMD = set32_include(splax, E.norm); 
              
                /* Pieces of this locus: */
                dg_locus_t ELO = dg_locus(E.norm, 2*E.cell);
                dg_locus_t EHI = dg_locus(E.norm, 2*E.cell + 1);
                dg_locus_t EMD = dg_locus(NrmMD, 2*E.cell + 1);
                
                /* Shapes of the pieces: */
                bz_patch_t bMD = bz_patch_facet_new(m, b->n, b->g, splj);
                bz_patch_t bLO = bz_patch_new(m, b->n, b->g);
                bz_patch_t bHI = bz_patch_new(m, b->n, b->g);

                 /* Tree stars of the pieces: */
                dg_tree_vec_t NLO, NHI, NMD;
                
                /* Split the tree star: */
                assert(psz[splax] == 1);
                dg_tree_pack_child(d, psz, NE, splax, 0, &NLO);
                dg_tree_pack_child(d, psz, NE, splax, 1, &NHI);

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

                /* Recurse for the separating item {S} of {F}. */
                psz[splax] = 2;
                dg_tree_pack_split(d, psz, NE, splax, &NMD);
                dg_enum_locus_leaves(d, EMD, r+1, &bMD, psz, &NMD, rMax, visit);
                psz[splax] = 1;

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

        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_patch_rdim_t 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_patch_rdim_t 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]. */