/* General tensor-style Bézier patches. */
/* Last edited on 2011-09-26 13:58:14 by stolfilocal */

#ifndef bz_patch_H
#define bz_patch_H

#include <psp_basic.h>
#include <psp_grid.h>
#include <bz_basic.h>
#include <double_array.h>

#include <interval.h>
#include <interval_array.h>
#include <box.h>
#include <stdint.h>

/* 
  BÉZIER REPRESENTATION FOR MULTIVALUED MULTIVARIATE POLYNOMIALS
  
  A /{d}-dimensional tensor-style Bézier patch/ {P} is a function from {\RR^d} into
  some real tensor space. Specifically, the patch {P} associates to each point
  {x[0..d-1]} of {\RR^d} a multi-index array {P(x:B)} of real numbers, whose
  elements are polynomial functions of the coordinates {x[0],..x[d-1]}.

  The map {P()} is specified by a collection of /Bézier coefficients/, /coeffs/
  for short, stored in an array {C} of real numbers. This array has {d+r}
  indices, for some {r} greater than 0. Each Bézier coefficient is a slice of
  the {C} array, with {r} indices, selected by the first {d} indices. Thus, for
  example, if {d=2} and {r=3}, the domain of {P} is {\RR^2}, each coefficient of
  {P} is an array with 3 indices, and so is any value {P(x:B)} of the patch. The
  coeff array {C} then has 5 indices, and element {C[p,q,r,s,t]} is element
  {[r,s,t]} of Bézier coefficient {[p,q]}.

  The map {x --> P(x:B)} is defined as a linear combination of /tensor-type
  Bézier basis elements/, each multiplied by one of the coeffs stored in {C}.
  Each basis element is the product of {d} univariate polynomials. Factor {i} of
  this product is polynomial of some degree {g[i]} on the coordinate {x[i]} of
  the argument. More precisely, the factor is a Bernstein-Bézier polynomial of
  degree {g[i]} and some index {e[i]} on the coordinate {x[i]}. The degree
  {g[i]} is one less than the size of the array {C} along its index {i}, and the
  index {e[i]} ranges over {0..g[i}}. The index vector {e[0..d-1]} identifies
  the Bézier basis element and selects the corresponding Bézier coeff in the {C}
  array.

  More formally, the value of {P(x)} on a point {x} of {R^d} is the sum
  of

    { C'[e] * PROD{ H_{e[i]}^{g[i]}(x[i]) : i = 0..d-1 } }

  over every {d}-index tuple {e} that is a prefix of a valid
  {d+r}-index tuple of {C}; where {C'[e]} is the slice of {C}
  whose first {d} indices is the tuple {e}, and 

    {H_k^m(z) == z^k*(1-z)^{m-k}*\choose(m,k)/2^m}

  is the Bernstein-Bezier polynomial of index {k} and degree {m}.

  In practice, the value of {P(x)} can be computed more efectively by
  successive interpolations along each coordinate (the DeCasteljau algorithm).

  Although the map {P()} is defined over all {R^d}, for many applications it is
  convenient to define the patch proper as the restriction of {P} to the unit
  {d}-cube {U^d = [0 _ 1]^d} of {R^d}. In particular, a Bezier coefficient
  {P.coef[e]} where every subscript {e[i]} is ether {0} or {g[i]} is a `corner'
  of the patch (the image of a corner of the cube {U^d}). The remaining coeffs
  (which exist only if some {g[i]} is greater than 1) affect the shape of the
  patch between those corners.

  An important example is the bicubic Bézier quadrilateral surface patch, widely
  used for freeform graphics modeling in 3D. It can be represented as a
  {bz_patch_t} {P} with {d=2}, {r=1}, and a coeff array {C} with three indices
  and size vector {(4,4,3)}. An argument {x} for {P} is a vector with two real
  numbers, the surface parameters {u} and {v}. Each value {P(x)} of the patch,
  and each Bézier coefficient, is a point of 3D space; that is, a vector
  (single-index array) with 3 elements, the spatial coordinates {X,Y,Z}.

  The gradient of such a surface patch can be represented by a {bz_patch_t}
  {db}, also with {d=2} but with {r=2}, so that the coeff array {C} has four
  indices. Each coeff, and the value {db(x)} of {db} at any point {x}, is a
  matrix with 2 indices and size vector {(3,2)}, containing the partial
  derivatives of the components of {P} at {x}; namely, {((dX/du, dX/dv), (dY/du,
  dY/dv), (dZ/du, dZ/dv))}. */

typedef psp_dim_t bz_patch_ddim_t;
  /* Dimension of the domain of a Bézier patch. */
 
#define bz_patch_MAX_DDIM (double_array_MAX_AXES)
  /* Maximum dimension for the domain of a Bézier patch. Note that some vectors
    in the representation are preallocated with this size, and that the number
    of coeffs grows exponentially with the domain dimension. */

typedef psp_dim_t bz_patch_rinds_t;
  /* Number of indices in the value (and in each coefficient) of a Bézier patch.
    A {bz_patch_rinds_t} of 0 means that the patch values and coefficients are
    scalars. */
 
#define bz_patch_MAX_RINDS (double_array_MAX_AXES)
  /* Maximum number of indices in the values of a Bézier patch. */
 
typedef struct bz_patch_t
  { bz_patch_ddim_t d;      /* Dimension of domain. */  
    bz_patch_rinds_t r;     /* Indices in values. */  
    double_array_t C;       /* The Bézier coefficients of the patch. */
  } bz_patch_t;
  /* A Bézier patch. */
  
const ix_size_t *bz_patch_get_dsize ( bz_patch_t *P );
  /* Returns the address of an array {dsz[0..P.d-1]} with the number of
    coefficients of patch {P} along each domain axis (which is 1 more than the degree of
    the corresponding argument coordinate). Namely, {dsz[i]} is the size of {P.C} along axis
    {i}, for {i} in {0..P.d-1}. */

const ix_size_t *bz_patch_get_rsize ( bz_patch_t *P );
  /* Returns the address of an array {rsz[0..P.r-1]} with the size of the values
    (and coefficients) of patch {P}, along each indexing axis. Namely, {rsz[i]}
    is the size of {P.C} along axis {P.d + i}, for {i} in {0..P.r-1}. */

void bz_patch_get_degrees ( bz_patch_t *P, bz_degree_t g[] );
  /* Stores in {g[0..P.d-1]} the degree vector of {P}; namely, stores
    into {g[i]} the degree of patch {P} along domain axis {i}. */

#define bz_patch_MAX_DEGREE (6)
  /* Maximum degre of a Bézier patch along each coordinate,
    just to be safe. Could be increased with modest harm. */
    
/*
  ALLOCATION */

bz_patch_t bz_patch_new ( bz_patch_ddim_t d, const bz_degree_t g[], bz_patch_rinds_t r, const ix_size_t rsz[] );
  /* Returns a Bézier patch {P} with {d}-dimensional domain and
    degrees {g[0..d-1]} along the {d} axes of the domain. The value of
    the patch at any point, as well as each Bézier coefficient, will
    be an array with {r} indices and {rsz[j]} elements along each
    index. If {rsz} is NULL, assumes {rsz[j] = 1} for all {j}. The
    coefficient array {P.C} will be newly allocated. */

bz_patch_t bz_patch_uniform_new ( bz_patch_ddim_t d, bz_degree_t g, bz_patch_rinds_t r, const ix_size_t rsz[] );
  /* Returns a Bézier patch {P} with {d}-dimensional domain and the
    same degree {g} along each axis of the domain. The parameters {r}
    and {rsz[0..r-1]} have the same meaning as in {bz_patch_new}. The
    coefficient array {P.C} will be newly allocated. */

bz_patch_t bz_patch_copy ( bz_patch_t *P );
  /* Returns a bewly-allocated Bézier patch with same
    attributes and coefficients as {P}, but not shared with any
    existing patch. */

void bz_patch_free ( bz_patch_t *P );
  /* De-allocates the coefficient vector and other internal storage
    areas of {P} (but not {P} itself). */

typedef ix_index_t bz_patch_cpindex_t;
  /* Each coeff has a {bz_patch_cpindex_t} along each domain axis,
    which ranges from 0 to the corresponding degree, inclusive. */

double_array_t bz_patch_value_new ( bz_patch_t *P );
  /* Returns a descriptor to a new array {F} where one can store 
    a value (or a coeff) of patch {P}. */

double_array_t bz_patch_get_coeff ( bz_patch_t *P, const bz_patch_cpindex_t e[] );
  /* Given a Bézier patch {P} of domain dimension {d = P.d} from {R^d} to {R},
    and a vector {e[0],..e[d-1]} of indices, returns the address of the coeff
    {P.coef[e]}. Note that the coeff itself is neither allocated nor copied;
    what is returned is just a descriptor for a slice of the control
    array {P.C}. */

void bz_patch_coeff_position ( bz_patch_t *P, const ix_index_t ix[], double x[] );
  /* Sets {x[0..P.d-1]} to the nominal position in {U^P.d}
    of the Bézier coefficient with indices {ix[0..P.d-1]}. */

void bz_patch_eval(bz_patch_t *P, const double x[], double_array_t *F);
  /* Stores into the array {*F} the value {P(x)} of the patch {P} at the point
    {x[0..d-1]} of {R^P.d}. The array {*F} must be allocated by the caller, with
    {P.r} indices and the proper size along each axis. */

/* 
   FACES AND FACETS */

bz_patch_t bz_patch_get_facet ( bz_patch_t *P, ix_axis_t i, interval_side_t dir );
  /* Given a Bézier patch {P} of domain dimension {d}, returns
    a patch {T} which is the restriction of {P} to the {d-1}-dimensional face
    (facet) of its nominal box {D} that is perpendicular to axis {i}, and is 
    located in the {dir} direction of that axis with respect to {D}. 
    
    The coeffs of patch {T} will be shared with patch {P}. */

bz_patch_t bz_patch_get_face ( bz_patch_t *P, const box_signed_dir_t dir[] );
  /* Returns a Bezier patch {T} which is the restriction of {P} to a 
    specified face {F} of the unit hypercube {U^d}.  The face is selected
    by its signature vector {dir[i]}, as in {box.h}. 
    
    The coeffs of patch {T} will be shared with patch {P}.
    
    The patch {T} must be allocated by the caller, with range dimension
    {T->r = P->r}. Its domain dimension {T->d} must be the 
    dimension of the face, i.e. the number of elements {dir[i]}
    which are equal to {SMD}; and its degree vector {T->g}
    must contain the elements of {P->g} that correspond to 
    those axes. */

bz_degree_t bz_patch_max_degree ( bz_patch_t *P );
  /* The maximum degree of {P} along any domain axis. */

void bz_patch_raise_degree ( bz_patch_t *P, bz_patch_t *T );
  /* 
    Stores in {T} a Bézier patch of degree vector {T->g} which coincides
    with the patch {P}.  The patch {T} must be allocated by the caller.
    It must have {T.d == P.d}, {T.r == P.r}, and must satisfy {T.sz[i] >= P.sz[i] > 0}
    for all {i} in {0..P.d-1}, and {T.sz[i] == P.sz[i]} for all {i} in {P.d..P.d+P.r-1},. */

void bz_patch_write ( FILE *wr, bz_patch_t *P, char *cmt );
  /* Writes the patch {P} to file {wr}, in a format that 
    can be read back with {bz_patch_read}.  The {cmt} string
    is written after the header line, each of its lines 
    prefixed with "# ". */

bz_patch_t bz_patch_read ( FILE *rd, char **cmtP ); 
  /* Reads from {rd} a bezier patch {P}. If {cmtP} is not mull,
    the comment after the header line, with leading "# "s 
    stripped, is stored in {*cmtP}. */

void bz_patch_debug ( FILE *wr, char *pref, bz_patch_t *P, char *suff ); 
  /* Prints to {wr} the patch {P}, preceded by {pref} and followed by {suff}. */

void bz_patch_arg_debug ( FILE *wr, char *pref, bz_patch_t *P, const double x[], char *suff );
  /* Prints to {wr} the argument {x[0..P.d-1]}, preceded by {pref} and followed by {suff}. */

/* AUXILIARY PROCEDURES */

bz_patch_t bz_patch_make_test ( bz_patch_ddim_t d, const bz_degree_t g[], bz_patch_rinds_t r, const ix_size_t rsz[] );
  /* Creates a Bézier patch of domain dimension {d}, degree vector {g[0..d-1]},
    {r} indices in value, and range size vector {rsz[0..r-1]}. */
    
#endif