/* See {dg_plot_bezier.h}. */
/* Last edited on 2009-05-18 14:20:25 by stolfi */

#define dg_plot_func_C_COPYRIGHT "Copyright © 2007 by the State University of Campinas (UNICAMP)."

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <math.h>

#include <bool.h>
#include <jsmath.h>
#include <affirm.h>
#include <rn.h>
#include <vec.h>
#include <pswr.h>
#include <pswr_color.h>
#include <pswr_iso.h>

#include <mdg_grid.h>
#include <dg_spline.h>
#include <dg_plot.h>

#include <bz_basic.h>

#include <dg_plot_2D_func.h>

/* INTERNAL DEFINITIONS */

#define DDIM (dg_plot_func_DDIM)
  /* Shorter name for the domain dimension. */

#define MAX_VALUES (dg_plot_func_MAX_VALUES)
  /* Shorter name for the max function values (besides shape coords). */

#define MIN_RDIM (dg_plot_func_DDIM)
  /* Minimum dimension for the range of {F}. */

#define MAX_RDIM (dg_plot_func_DDIM + dg_plot_func_MAX_VALUES)
  /* Maximum dimension for the range of {F}. */

#define dg_plot_func_NF (9)
  /* Number of faces in a cell of the domain grid ({3^DDIM}). */

#define dg_plot_func_TOL_MM (0.5)
  /* Flatness tolerance (mm). */

#define dg_plot_func_DEBUG_SHAPE FALSE
  /* TRUE to debug {make_shape} instead of {dg_enum_faces}. */

/* IMPLEMENTATIONS */

void dg_plot_2D_func
  ( PSStream *ps, 
    bz_patch_t *bz,
    bz_patch_rdim_t nf,
    dg_plot_2D_func_t *func,
    dg_tree_t G,
    mdg_rank_t plotDp,
    bool_t isolines
  )
  {
    /* Check {bz} and provide the default, if needed: */
    bz_patch_t bz_can;  /* Patch that describes the canonical box, or NULL. */
    bz_can.c = NULL;
    if (bz == NULL)
      { /* Provide a default Bézier patch for the canonical rectangle: */
        interval_t box[DDIM];
        mdg_cell_box_canonical(DDIM, mdg_ROOT_CELL, box);
        bz_can = bz_patch_from_box(DDIM, box);
        bz = &bz_can;
      }
    demand(bz->m == DDIM, "Bezier patch has wrong domain dimension");
    demand(bz->n >= DDIM, "Bezier patch has bad range dimension");
    
    /* Check {nf}: */
    if (func == NULL) 
      { demand(nf == 0, "null function should have zero range dimension"); }
    else
      { demand(nf >= 0, "function has bad range dimension"); }
    
    /* Compute the total range dimension of {F}: */
    int nt = bz->n + nf;  /* Total range dimension (including shape coords). */
    int nv = nt - DDIM;   /* Number of function values (excluding shape coords). */
    demand(nt > DDIM, "total range dimension too low");
    demand(nv <= MAX_VALUES, "too many function values");
    
    /* Number of function values: */
    if (isolines)
      { demand(nv == 1, "isolines/bands require exactly one function value"); }
    else
      { demand(nv <= 3, "too many function values for shading"); }

    /* Choose the nominal range for function values: */
    double fMax = 1.0;
    double fMin = -fMax;

    /* Choose the isoline spacing: */
    double fStep = 0.2;
    double fSync = fStep/2;

    /* Compute the number of isolines: */
    int kMin = pswr_inf_isoline(fSync, fStep, -fMax);
    int kMax = pswr_sup_isoline(fSync, fStep, +fMax);

    /* Build the color band tables: */
    int N;
    double *Rtb, *Gtb, *Btb;
    pswr_make_color_table(
      fSync, fStep, kMin, kMax,
      0.000, 0.500, 1.000,
      1.000, 1.000, 1.000,
      1.000, 0.167, 0.000,
      &N, &Rtb, &Gtb, &Btb);

    auto void plot_2D_item(mdg_dim_t d, dg_locus_t E, mdg_rank_t r, bz_patch_t *b, dg_tree_vec_t *na);
      /* Procedure called to plot a locus {E} from level {r} of the grid {G}. The patch {b}
        describes the shape of the item. The vector {na} contains the tree pack 
        corresponding to the star of {E}. */ 

    void plot_2D_item(mdg_dim_t d, dg_locus_t E, mdg_rank_t r, bz_patch_t *b, dg_tree_vec_t *na)
      { int mb = b->m; /* Dimension of item. */
        int nb = b->n; /* Dimension of range (coords). */
        /* fprintf(stderr, "+ plot_2D_item(E = <%d:%lld> r = %d)\n", E.norm, E.cell, r); */
        /* fprintf(stderr, "item shape\n"); */
        /* bz_patch_print(stderr, b, "%6.2f"); */
        affirm(mb <= 2, "bad domain dimension");
        /* Get root-relative coordinates of locus: */
        interval_t EB[DDIM];
        dg_locus_box_root_relative(DDIM, E, EB);
        /* Compute the number of distinct corners {nc}: */
        int nc = ipow(2, mb);
        /* Get the domain coords of corners {x[k*DDIM+i]}, for {i=0..DDIM-1}, {k=0..nc-1}: */
        double x[nc*DDIM];
        int i, j, k;
        fprintf(stderr, "  mb = %d", mb);
        for (k = 0; k < nc; k++)
          { double *xk = &(x[k*DDIM]);
            fprintf(stderr, "  x[%d] = (", k);
            for (j = 0; j < mb; j++)
              { int i = dg_locus_spanning_axis(DDIM, E, mb-1-j);
                int side = ((k >> j) & 1);
                xk[j] = EB[i].end[side];
                fprintf(stderr, " %d:%d %7.4f", i, side, xk[j]); 
              }
            fprintf(stderr, " )\n");
          }
        /* Get the corner positions and values {c[k*nt+i]}, for {i=0..nt-1}: */
        double c[nc*nt];
        for (k = 0; k < nc; k++)
          { double *ck = &(c[k*nt]);
            /* Extract corner {k} of the Bézier patch, save in {ck[0..nb-1]}: */
            double *bck = &(b->c[k*nb]);
            for (i = 0; i < nb; i++) { ck[i] = bck[i]; }
            /* Evaluate {func} at the domain point {xk[0..DDIM-1], save in {ck[nb..nt-1]}: */
            double *xk = &(x[k*DDIM]);
            func(xk, DDIM, &(ck[nb]), nf);
          }
        /* Now plot the item: */
        switch(mb)
        { case 0: 
            { /* Plot the vertex: */
              dg_plot_2D_vertex(ps, c, nt);
            }
            break;
          case 1: 
            affirm(b->g[0] == 1, "bad patch degree");
            { /* Plot the edge: */
              double *c0 = &(c[0*nt]);
              double *c1 = &(c[1*nt]);
              dg_plot_2D_edge(ps, c0, c1, nt); 
            }
            break;
          case 2: 
            affirm((b->g[0] == 1) && (b->g[1] == 1), "bad patch degree");
            { /* Choose the plot subdivision depth for this face: */
              int subDp = plotDp - r + 1;
              if (subDp < 0) { subDp = 0; }
              if (subDp > 3) { subDp = 3; }
              /* Plot the face: */
              double *c00 = &(c[0*nt]);
              double *c01 = &(c[1*nt]);
              double *c10 = &(c[2*nt]);
              double *c11 = &(c[3*nt]);
              if (isolines)
                { dg_plot_2D_face_bands
                    ( ps, c00, c01, c10, c11, nt, subDp, fSync, fStep, kMin, kMax, Rtb,Gtb,Btb );
                  dg_plot_2D_face_isolines
                    ( ps, c00, c01, c10, c11, nt, subDp, fSync, fStep, kMin, kMax );
                }
              else
                { dg_plot_2D_face_shade
                    ( ps, c00, c01, c10, c11, nt, subDp, fMin, fMax );
                }
            }
            break;
          default: affirm(FALSE, "bad item dimension"); 
        }
        /* fprintf(stderr, "- plot_2D_item\n"); */
      }

    dg_locus_t E0 = dg_locus(0, mdg_ROOT_CELL);
    dg_tree_vec_t na = dg_tree_vec_new(1);

    fprintf(stderr, "grid cells = %d\n", G->nodes);

    mdg_grid_size_t psz[DDIM];
    { int k; for (k = 0; k < DDIM; k++) { psz[k] = 1; } } 
    
    na.e[0] = G;
    dg_enum_locus_leaves(DDIM, E0, 0, bz, psz, &na, plotDp, plot_2D_item);
      
    free(Rtb); free(Gtb); free(Btb);
    if (bz_can.c != NULL) { bz_patch_free(bz_can); }
  }