/*

                              DISCLAIMER
                              ==========

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.   

    If the software is modified by someone else and passed on, we, the authors
    want its recipients to know that what they have is not the original, so
    that any problems introduced by others will not reflect on the original
    authors' reputations.
*/                                            

#include "nurbh.h"
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_torus(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Create a torus.
*/
{
  Pmatrix3 xx;
  PR_nurb *nrb1;
  Pvector3 v1;
  nrb1 = NULL;
  nrb_allocatenurb(&nrb1);
  /*create a circle*/
  nrb_circle(nrb1);
  /*shift in X*/
  /*rotate to XZ plane*/
  v1 = ptk_vector3(2.0, 0.0, 0.0);
  ptk_shift3(&v1 , PREPLACE, xx);
  ptk_rotate3(90.0, 1, PPRECONCATENATE, xx);

  nrb_xform(nrb1, xx);

  nrb_revolve(nrb1, nrb);
  nrb_deallocatenurb(&nrb1);

}  /*nrb_torus*/
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_polyline(C(Pint) npts, C(Ppoint3 *)pp, C(PR_nurb *)nrb)
PreANSI(Pint npts) 
PreANSI(Ppoint3 *pp) 
PreANSI(PR_nurb *nrb)
/*
Converts a POLYLINE/POLYMARKER to a NURB Curve - each point is a double knot
*/
{
  Pint i;

  /*nrb_polyline*/
  nrb->pf_u.pf_k = 2;   /*Second Order = Linear*/
  nrb->pf_u.pf_n = npts * 2;
  nrb->pf_u.pf_nt = nrb->pf_u .pf_n + nrb->pf_u.pf_k;

  nrb_allocateknots(&nrb->pf_u);

  nrb->pf_v.pf_k = 0;   /*Zero Order in Y direction*/
  nrb->pf_v.pf_n = 1;
  nrb->pf_v.pf_nt = 1;
  nrb_allocateknots(&nrb->pf_v);
  nrb->pf_v.pf_kk->knots[0] = 0.0;

  nrb_allocatepts(npts * 2, &nrb->pf_ppp);

  for (i = 1; i <= npts; i++) {   
    nrb->pf_ppp->pts[i * 2 - 2] = ptk_pt3topt4(&pp[i - 1]);
	/*Duplicate each point*/
    nrb->pf_ppp->pts[i * 2 - 1] = ptk_pt3topt4(&pp[i - 1]);
    nrb->pf_u.pf_kk->knots[i * 2 - 2] = i - 1.0;
	/*Duplicate Knots at each point*/
    nrb->pf_u.pf_kk->knots[i * 2 - 1] = i - 1.0;
  }


}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_square(C(PR_nurb *)nrb1)
PreANSI(PR_nurb *nrb1)
/*
    produce a square in the xy plane about the origin
    centre=(0.0,0), length of side = 1
*/
{
  nrb1->pf_u.pf_k = 2;
  nrb1->pf_v.pf_k = 0;
  nrb1->pf_u.pf_n = 8;
  nrb1->pf_v.pf_n = 1;
  nrb1->pf_u.pf_nt = 10;
  nrb1->pf_v.pf_nt = 1;

  nrb_allocateknots(&nrb1->pf_u);
  nrb_allocateknots(&nrb1->pf_v);
  nrb_allocatepts(nrb1->pf_u.pf_n * nrb1->pf_v.pf_n, &nrb1->pf_ppp);

  nrb1->pf_ppp->pts[0] = tp_litv4(0.5, 0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[7] = nrb1->pf_ppp->pts[0];

  nrb1->pf_ppp->pts[1] = tp_litv4(-0.5, 0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[2] = nrb1->pf_ppp->pts[1];

  nrb1->pf_ppp->pts[3] = tp_litv4(-0.5, -0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[4] = nrb1->pf_ppp->pts[3];

  nrb1->pf_ppp->pts[5] = tp_litv4(0.5, -0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[6] = nrb1->pf_ppp->pts[5];

  nrb1->pf_u.pf_kk->knots[0] = 0.0;
  nrb1->pf_u.pf_kk->knots[1] = 0.0;
  nrb1->pf_u.pf_kk->knots[2] = 0.25;
  nrb1->pf_u.pf_kk->knots[3] = 0.25;
  nrb1->pf_u.pf_kk->knots[4] = 0.5;
  nrb1->pf_u.pf_kk->knots[5] = 0.5;
  nrb1->pf_u.pf_kk->knots[6] = 0.75;
  nrb1->pf_u.pf_kk->knots[7] = 0.75;
  nrb1->pf_u.pf_kk->knots[8] = 1.0;
  nrb1->pf_u.pf_kk->knots[9] = 1.0;

  nrb1->pf_v.pf_kk->knots[0] = 0.0;   
}  /* of nrb_square */
/*-------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_triangle(C(PR_nurb *)nrb1)
PreANSI(PR_nurb *nrb1)
/*
    produce a triangle in the xy plane about the origin
*/
{
  nrb1->pf_u.pf_k = 2;
  nrb1->pf_v.pf_k = 0;
  nrb1->pf_u.pf_n = 6;
  nrb1->pf_v.pf_n = 1;
  nrb1->pf_u.pf_nt = 8;
  nrb1->pf_v.pf_nt = 1;

  nrb_allocateknots(&nrb1->pf_u);
  nrb_allocateknots(&nrb1->pf_v);
  nrb_allocatepts(nrb1->pf_u.pf_n * nrb1->pf_v.pf_n, &nrb1->pf_ppp);

  nrb1->pf_ppp->pts[0] = tp_litv4(0.0, 0.366, 0.0,1.0);
  nrb1->pf_ppp->pts[5] = nrb1->pf_ppp->pts[0];

  nrb1->pf_ppp->pts[1] = tp_litv4(-0.5, -0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[2] = nrb1->pf_ppp->pts[1];

  nrb1->pf_ppp->pts[3] = tp_litv4(0.5, -0.5, 0.0,1.0);
  nrb1->pf_ppp->pts[4] = nrb1->pf_ppp->pts[3];

  nrb1->pf_u.pf_kk->knots[0] = 0.0;
  nrb1->pf_u.pf_kk->knots[1] = 0.0;
  nrb1->pf_u.pf_kk->knots[2] = 1.0 / 3.0;
  nrb1->pf_u.pf_kk->knots[3] = 1.0 / 3.0;
  nrb1->pf_u.pf_kk->knots[4] = 2.0 / 3.0;
  nrb1->pf_u.pf_kk->knots[5] = 2.0 / 3.0;
  nrb1->pf_u.pf_kk->knots[6] = 1.0;
  nrb1->pf_u.pf_kk->knots[7] = 1.0;

  nrb1->pf_v.pf_kk->knots[0] = 0.0;   
}  /* of nrb_triangle */
/*-------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_vase(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Produce a vase.
*/
{
  PR_nurb *nrb1;
  Ppoint4 temp;

  /*nrb_vase*/
  nrb1 = NULL;
  nrb_allocatenurb(&nrb1);

  nrb1->pf_u.pf_k = 3;
  nrb1->pf_v.pf_k = 1;
  nrb1->pf_u.pf_n = 9;
  nrb1->pf_v.pf_n = 1;
  nrb1->pf_v.pf_nt = 1;
  nrb1->pf_u.pf_nt = 12;

  nrb_allocateknots(&nrb1->pf_u);
  nrb_allocateknots(&nrb1->pf_v);
  nrb_allocatepts(nrb1->pf_u.pf_n * nrb1->pf_v.pf_n, &nrb1->pf_ppp);

  nrb1->pf_u.pf_kk->knots[0] = 0.0;
  nrb1->pf_u.pf_kk->knots[1] = 0.0;
  nrb1->pf_u.pf_kk->knots[2] = 0.0;
  nrb1->pf_u.pf_kk->knots[3] = 1.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[4] = 2.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[5] = 3.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[6] = 4.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[7] = 5.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[8] = 6.0 / 7.0;
  nrb1->pf_u.pf_kk->knots[9] = 1.0;
  nrb1->pf_u.pf_kk->knots[10] = 1.0;
  nrb1->pf_u.pf_kk->knots[11] = 1.0;

  nrb1->pf_v.pf_kk->knots[0] = 0.0;

  temp = tp_litv4(0.0, 0.0, 0.0,1.0);

  nrb1->pf_ppp->pts[0] = temp;

  temp = tp_litv4(0.4, 0.0, 0.0,1.0);
  nrb1->pf_ppp->pts[1] =temp;

  temp = tp_litv4(0.4, 0.0, 0.0,1.0);
  nrb1->pf_ppp->pts[2] =temp;

  temp = tp_litv4(0.1, 0.0, 0.25,1.0);
  nrb1->pf_ppp->pts[3] =temp;

  temp = tp_litv4(0.4, 0.0, 0.5,1.0);
  nrb1->pf_ppp->pts[4] =temp;

  temp = tp_litv4(0.4, 0.0, 1.0,1.0);
  nrb1->pf_ppp->pts[5] =temp;

  temp = tp_litv4(0.35, 0.0, 1.0,1.0);
  nrb1->pf_ppp->pts[6] =temp;

  temp = tp_litv4(0.35, 0.0, 0.5,1.0);
  nrb1->pf_ppp->pts[7] =temp;

  temp = tp_litv4(0.0, 0.0, 0.5,1.0);   
  nrb1->pf_ppp->pts[8] =temp;

  nrb_revolve(nrb1, nrb);
  nrb_deallocatenurb(&nrb1);

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_makehjknots(C(PR_nurb *)nrbin)
PreANSI(PR_nurb *nrbin)
/*
Make Harley-Judd knots.
*/
{
  Pint jj, ii, nn, k;
  Pfloat delta, sum;
  PR_dir *tdirts;
  Pint mm1;
  PR_pts *tpts1;

  tdirts = &nrbin->pf_u;
  nn = tdirts->pf_n + 1;
  k = tdirts->pf_k;
  sum = 0.0;
  delta = 0.0;
  nrb_allocateknots(&nrbin->pf_u);
  tdirts->pf_kk->knots[k - 1] = 0.0;
  for (ii = k + 1; ii <= nn; ii++) {
    delta = 0.0;
    mm1 = ii - k + 2;
    for (jj = ii - k; jj <= mm1; jj++) {
      tpts1 = nrbin->pf_ppp;
      delta += tp_modv4(tp_sub4(tpts1->pts[jj], tpts1->pts[jj - 1]));
    }
    sum += delta;
    tdirts->pf_kk->knots[ii - 1] = tdirts->pf_kk->knots[ii - 2] + delta;
  }

  /*Now renormalise to 0..1*/
  for (ii = k; ii < nn; ii++)
    tdirts->pf_kk->knots[ii] /= sum;

  /*Force Multiplicty at ends*/
  for (ii = 0; ii <= k - 2; ii++) {   
    tdirts->pf_kk->knots[ii] = tdirts->pf_kk->knots[k - 1];
    tdirts->pf_kk->knots[nn + ii] = tdirts->pf_kk->knots[nn - 1];
  }

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_interpolate(C(PR_nurb *)nrbin, C(Pfloat) tol, 
                            C(PR_nurb *)nrbout)
PreANSI(PR_nurb *nrbin) 
PreANSI(Pfloat tol) 
PreANSI(PR_nurb *nrbout)
/*
Interpolate a NURB to a given tolerance through a set of data points.
*/
{
  Pint iter, i, ll;
  Pfloat delta,xx;
  PR_pts *dd;
  Ppoint4 t1, t2, d;
  Pint nxin, nyin;
  boolean converged;
  PR_dir *tdirts;
  Pint mm, mm1;

#define NRBPT(i,ll) nrbin->pf_ppp->pts[nrb_index((i), (ll), nxin) - 1]

  /*Copy Data*/
  nrb_copy(nrbin, nrbout);

  /*Make Harley-Judd Knots*/
  nrb_makehjknots(nrbout);

  /*Copy Knots Y -> X*/
  nrb_interchange(nrbout);
  nrb_deallocateknots(&nrbout->pf_v);
  nrbout->pf_v.pf_kk = nrbout->pf_u.pf_kk;
  nrbout->pf_u.pf_kk = NULL;

  /*Copy Y knots -> X*/
  nrb_allocateknots(&nrbout->pf_u);
  tdirts = &nrbout->pf_u;

  mm = tdirts->pf_nt;
  for (i = 0; i < mm; i++)
    tdirts->pf_kk->knots[i] = nrbin->pf_v.pf_kk->knots[i];
  dd = NULL;

  nrb_allocatepts(nrbout->pf_u.pf_n, &dd);
  nrb_allocatepts(nrbout->pf_u.pf_n * nrbout->pf_v.pf_n, &nrbout->pf_ppp);

  nxin = nrbin->pf_u.pf_n;
  nyin = nrbin->pf_v.pf_n;

/* Macro for approximate (and hopefully converged) solution point */

#define NRBAPP(i,ll) nrbout->pf_ppp->pts[nrb_index((ll), (i), nyin) - 1] 

  for (ll = 1; ll <= nyin; ll++) {   /*Loop over each Y row - store in X*/
    for (i = 1; i <= nxin; i++) {   /*Initial Approximation*/
       NRBAPP(i,ll) = NRBPT(i,ll);
    }

    converged = FALSE;
    iter = 0;

    while (!converged)
    {  /*Another Iteration - Note first and last points unmodified*/
      delta = 0.0;   /*Maximum Change*/
      iter++;
      for (i = 2; i < nxin; i++) {
	d  = tp_sub4(NRBPT(i,ll),    NRBAPP(i,ll));
	t1 = tp_linsum4(0.5,NRBAPP(i-1,ll), 0.5, NRBAPP(i+1,ll));
	t2 = tp_linsum4(0.5, NRBPT(i,ll),-0.5, t1);
	d = tp_add4(d, t2);
	xx = tp_dotv4(d, d);
	delta = P_max(delta, xx);
	dd->pts[i - 1] = d;
      }

      /*Now update control points*/
      for (i = 2; i < nxin; i++) 
          NRBAPP(i,ll) = tp_add4(NRBAPP(i,ll),dd->pts[i - 1]);

      /*Converged?*/
      converged = (delta < tol || iter > 100);
      fprintf(stderr,"nrb_interpolate:%3ld %12.4f %12.4f\n",iter,delta,tol);
    }
  }
  /*y loop*/
  nrb_deallocatepts(&dd);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_spli2d(C(PR_nurb *)nrbin, C(Pint) mm, C(Pfloat *)tpts, 
                       C(PR_nurb *)nrbout)
PreANSI(PR_nurb *nrbin) 
PreANSI(Pint mm) 
PreANSI(Pfloat *tpts) 
PreANSI(PR_nurb *nrbout)
/*
Documentation - to be completed
*/
{
  Pfloat Q[PC_MaxNetpts], work[PC_MaxNetpts];
  Pint iflag;

  /*Copy Data*/
  nrb_copy(nrbin, nrbout);

  /*Make Harley-Judd Knots*/
  nrb_makehjknots(nrbout);

  /*Copy Knots X -> Y*/
  nrb_interchange(nrbout);
  nrb_deallocateknots(&nrbout->pf_v);
  nrbout->pf_v.pf_kk = nrbout->pf_u.pf_kk;

  nrbout->pf_u.pf_kk = NULL;

  /*Copy Y knots -> X*/
  nrb_copyknots(&nrbin->pf_v, &nrbout->pf_u);

  spli2d(tpts, (Pfloat*)nrbin->pf_ppp->pts,
	 nrbout->pf_v.pf_kk->knots, nrbout->pf_v.pf_n, nrbout->pf_v.pf_k, mm,
	 4, nrbout->pf_u.pf_n, work, Q, (Pfloat*) nrbout->pf_ppp->pts, &iflag);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_tessalate(C(PR_nurb *)nrb1, C(Pint) nx, C(PR_nurb *)nrb2)
PreANSI(PR_nurb *nrb1)
PreANSI(Pint nx) 
PreANSI(PR_nurb *nrb2)
/*
Increase the refinement of the nrb1 by nx by inserting new knots.
*/
{
  Pfloat tmin, tmax;


  tmin = nrb1->pf_u.pf_kk->knots[0];
  tmax = nrb1->pf_u.pf_kk->knots[nrb1->pf_u .pf_nt - 1];

  nrb2->pf_u.pf_nt = P_max(nx, nrb1->pf_u .pf_k * 2 + 1);
  nrb2->pf_u.pf_k = nrb1->pf_u.pf_k;
  nrb2->pf_u.pf_n = nrb1->pf_u.pf_nt - nrb2->pf_u.pf_k;

  /*Copy Y knots*/
  nrb_copyknots(&nrb1->pf_v, &nrb2->pf_v);

  nrb_makeknots(tmin, tmax, &nrb2->pf_u);
  nrb_unionknots(&nrb1->pf_u, &nrb2->pf_u , &nrb2->pf_u);

  nrb_osloc(nrb1, nrb2);   


}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_drawcnet(C(PR_nurb *)nrb, C(Pint) start, C(PE_dir) xory)
PreANSI(PR_nurb *nrb) 
PreANSI(Pint start) 
PreANSI(PE_dir xory)
/*
Draw the control polygon of a curve.
*/
{
  /*Draw the control Net*/
  /* Should add a parameter to control
     1) Symbols at the points
     2) annotation at the points
     3) Join points by lines
   */
  Ppoint3 pts[PC_MaxNetpts];
  Pint stride, ii, npts, i;

  if (xory == PVudir) {
    npts = nrb->pf_u.pf_n;
    stride = 1;
  } else {
    npts = nrb->pf_v.pf_n;
    stride = nrb->pf_u.pf_n;
  }

  for (i = 0; i < npts; i++) {
    ii = start + i * stride;
    pts[i] = ptk_pt4topt3(&(nrb->pf_ppp->pts[ii - 1]));
  }
#ifdef __TURBOC__
#elif SUN
#else
  PPolyline3(npts, pts);
  PPolyMarker3(npts, pts);   
#endif
}  /* nrb_drawcnet */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_drawc(C(PR_nurb *)nrb, C(Pfloat) t, C(PE_dir) xory)
PreANSI(PR_nurb *nrb) 
PreANSI(Pfloat t) 
PreANSI(PE_dir xory)
/*
Draw a NURB curve together with its control polygon.
*/
{
  /*Draw a curve
    This will calculate the size of the refined polygon
    create a refined polygon and draw it*/
  PR_nurb *nrbR, *nrbS;
  Pfloat tmin, tmax;

  nrbR = NULL;
  nrb_allocatenurb(&nrbR);
  nrbS = NULL;
  nrb_allocatenurb(&nrbS);

  nrb_copy(nrb, nrbR);

  if (nrb->pf_v.pf_n <= 1)  /*Its already a curve!*/
    nrb_copy(nrb, nrbS);
  else {  /*Make it into a curve*/
    if (xory == PVudir)
      nrb_transpose(nrbR);
    nrb_evaluate(nrbR, t, nrbS);
  }

  /*draw control net first */
  /*penenq(icol, width, itype);*/
  /*pensel(1,width*3,3);*/

  nrb_drawcnet(nrbS, 1, xory);   /*refine to maximum size knot vector*/


  tmin = nrbS->pf_u.pf_kk->knots[0];
  tmax = nrbS->pf_u.pf_kk->knots[nrbS->pf_u .pf_nt - 1];

  nrbR->pf_u.pf_nt = nrb_numsubdivs(nrbS, PVudir, 0.1);

  fprintf(stderr,"%s%4ld%s\n", "NRB_DRAWC - refined to ",nrbR->pf_u.pf_nt," points");

  nrb_makeknots(tmin, tmax, &nrbR->pf_u);
  nrb_unionknots(&nrbS->pf_u, &nrbR->pf_u , &nrbR->pf_u);
  nrb_osloc(nrbS, nrbR);
  nrb_interchange(nrbR);   

  /*pensel(2,width,3);*/

  nrb_drawcnet(nrbR, 1, xory);

  /*nrb_dump(stdout,nrbr^);*/

  nrb_deallocatenurb(&nrbR);
  nrb_deallocatenurb(&nrbS);
}  /* nrb_drawc */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_drawsnet(C(PR_nurb *)cont_grid)
PreANSI(PR_nurb *cont_grid)
/*
Draw the control polygon of a surface.
*/
{
  /*
    draw a quadrilateral mesh to display control points
   */

  Pint i, j, mm;

  /*printf("NRB_DRAWSNET\n");
   nrb_dump(stdout,cont_grid);
  */

  if (cont_grid->pf_v.pf_nt > 1) {
    mm = cont_grid->pf_u.pf_n;
    for (i = 1; i <= mm; i++)
      nrb_drawcnet(cont_grid, i, PVvdir);
  }

  mm = cont_grid->pf_v.pf_n;
  for (j = 1; j <= mm; j++) {   
    i = ni(1, j, cont_grid->pf_u.pf_n);
    nrb_drawcnet(cont_grid, i, PVudir);
  }

}  /* of nrb_drawsnet */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_draws(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Draw a NURB surface together with its control polygon.
*/
{
  PR_nurb *nrbR, *nrbs;
  Pint nx, ny;

  nrbR = NULL;
  nrbs = NULL;

  nrb_allocatenurb(&nrbR);
  nrb_allocatenurb(&nrbs);

  /* draw control net first */
  /*penenq(icol, width, itype);*/
  /*pensel(1,width*3,3);*/

  nrb_drawsnet(nrb);

  /*refine to full size*/

  nx = nrb_numsubdivs(nrb, PVudir, 10.0);
  ny = nrb_numsubdivs(nrb, PVvdir, 10.0);

  nrb_tessalate(nrb, nx, nrbR);
  nrb_tessalate(nrbR, ny, nrbs);

  /*pensel(2,width,3);*/
  nrb_drawsnet(nrbR);

  nrb_deallocatenurb(&nrbR);
  nrb_deallocatenurb(&nrbs);
}
/*---------------------------------------------------------------------------*/
