/*  
    C code from An Introduction to NURBS
    by David F. Rogers. Copyright (C) 2000 David F. Rogers,
    All rights reserved.

    Book reference: Chapter 7, Section 7.10, p. 248

    name: rbsurf.c
 
    programmer: David F. Rogers
 
    date: 25 December 1989, modified 16 April 1990 and July 2000
 
    history: Developed from the algorithm in MECG 2nd and the original
             Fortran algorithm called bsurf.for (See the ICCAS'82 paper
             by Rogers and Satterfield). Modified for correct incremental
             calculation of change in homogeneous coordinate (see the
             Rogers and Adlum paper in the Bezier 80th birthday issue
             of CAD.
 
    language: C
 
    purpose:  Subroutine to calculate a Cartesian product rational B-spline
              surface using open uniform knot vectors.
 
    method: During the first pass through the algorithm calculate and
            hold basis and sum function in arrays. Subsequent passes
            do not require recalculation of these functions.
 
            Use an incremental method to recalculate the surface while
            a single defining net point or homogeneous weight is changed
            dynamically.
 
            Note: this routine is 20 times faster than that in the
            pseudocode in the book.
 
 
    subroutines called: basis.c, knot.c, sumrbas.c
 
    header files: rbsurf.h
 
    variable list:
 
    b[]         = array containing the polygon net points
                  b[1] = x-component
                  b[2] = y-component
                  b[3] = z-component
                  b[4] = h-component
                  Note: Bi,j = b[] has dimensions of n*m*3 with j varying fastest
                        The polygon net is n x m
    itest       = value indicating whether the surface set-up parameters changed
                  p_itest = a pointer to itest
    k           = order in the u direction
    l           = order in the w direction
    mbasis[]    = array containing the nonrational basis functions for one value of w (see \eq{5--84})
                  it has dimensions of mpts
    mjlw[]      = array containing the nonrational basis functions for all values of w (see \eq{5--84})
                  it has dimensions of p2*mpts
    mpts        = the number of polygon vertices in w direction
    nbasis[]    = array containing the nonrational basis functions for one value of u (see \eq{5--84})
                  it has dimensions of npts
    niku[]      = array containing the nonrational basis functions all values of u (see \eq{5--84})
                  it has dimensions of p1*npts
    npts        = the number of polygon vertices in u direction
    p1          = number of parametric lines in the u direction
    p2          = number of parametric lines in the w direction
    q[]         = array containing the resulting surface
                  q[1] = x-component
                  q[2] = y-component
                  q[3] = z-component
                       for a fixed value of u the next m elements contain
                       the values for the curve q[u[sub i],w] q has dimensions
                       of p1*p2*3. The display surface is p1 x p2
    rsumij[i3] = array containing all the sum values for all u,w. it has 
                 dimensions of p1*p2
 
                 Note: the BSSD program should be set up to allow 6th order
                 rational B-splines is both u and w. The maximum values of
                 the appropriate parameters should be
                     npts =  50 ... max number of net points in the u direction
                     mpts =  20 ... max number of net points in the w direction
                     p1   = 150 ... max number of parametric lines in the u direction
                     p2   =  60 ... max number of parametric lines in the w direction 
                     k    =   6 ... max order in the u direction
                     l    =   6 ... max order in the w direction
                 Since the zeroth element of the array is ignored 1 must be
                 added to all array dimensions. Thus q[27001], b[4001],
                 nbasis[51],mbasis[21],niku[7501],mulw[1201],rsumij[9001],
                 x[57],y[27]
*/
 
#include <stdio.h>
#include "rbsurf.h"
 
void rbsurf(b,k,l,npts,mpts,p1,p2,p_itest,ibnum,bold,niku,mjlw,rsumij,savrsumij,q)
 
int k,l;
int npts,mpts;
int p1,p2;
int *p_itest;
int ibnum;
 
float b[],q[];
float bold[];
float niku[];
float mjlw[];
float rsumij[];
float savrsumij[];
 
{
 
/*   allows for 20 data points with basis function of order 5 */
 
    int i,ibas,j,j1,jbas;
    int icount,scount,tcount,i1,i2,i3,i1inc;
    int uinc,winc;
    int ninc,minc;
    int nplusc,mplusc;
    int iindex,jindex;
    int x[30],y[30];
 
    float nbasis[30],mbasis[30];
    float pbasis;
    float sum,sumrbas();
    float u,w;
    float stepu,stepw;
    float bx,by,bz,bh;
    float sumnew,sumold,sumratio;
 
/*  int stack_probe;
    fprintf(stderr, "rbsurf:\t\t&stack_probe = %lu\n", (long)&stack_probe); */
/*
    printf("in rbsurf.c \n");
*/
    nplusc = npts + k;
    mplusc = mpts + l;
/*  
    check if a new surface needs to be generated.
     
    itest must be equal to zero initially to get the first surface.  If the
    number of net points, the order, or the number of parametric lines
    change in either the u or w parametric directions then a new surface
    must be generated. Otherwise the basis functions and the summation
    functions do not change when a defining polygon vertex is moved. When a
    homogeneous weight factor is changed the sum function must also be
    recalculated. However, the non-rational basis functions need not be
    recalculated. The assumption is that only one thing will change at a
    time. itest is a measure of the change.  If this is not the case then
    the calling program must force itest to zero to correctly generate the
    surface.
*/
 
    if (*p_itest != (nplusc + mplusc + p1 + p2)){
/* 
    printf("calculating basis sum functions \n");
*/ 
/*  zero the arrays */
 
        for (i = 1; i <= nplusc; i++){
            x[i] = 0.;
        }
        for (i = 1; i <= mplusc; i++){
            y[i] = 0.;
        }
        for (i = 1; i <= npts; i++){
            nbasis[i] = 0.;
        }
        for (i = 1; i <= mpts; i++){
            mbasis[i] = 0.;
        }
 
/*   generate the open uniform knot vectors */
 
        knot(npts,k,x);       /*  calculate u knot vector */
        knot(mpts,l,y);       /*  calculate w knot vector */
 
/*
        printf("u-knot vector = ");
        for (i = 1; i <= nplusc; i++){
            printf("  %d",x[i]);
        }
        printf("\n\n");
 
        printf("w-knot vector = ");
        for (i = 1; i <= mplusc; i++){
            printf("  %d",y[i]);
        }
        printf("\n\n");
*/
        icount = 1;
 
/*  calculate and store the basis functions */
 
        stepu = (float)x[nplusc]/(float)(p1-1);
        stepw = (float)y[mplusc]/(float)(p2-1);
/*
        printf("stepu = %f \n",stepu);
        printf("stepw = %f \n\n",stepw);
*/
        i1 = 1;
        i2 = 1;
        i3 = 1;
 
/*  calculate the Nik's at each u parametric value */
 
        u = 0.;
        for (uinc = 1; uinc <= p1; uinc++){
            if ((float)x[nplusc] - u < 5e-6){
                u = (float)x[nplusc];
            }
            basis(k,u,npts,x,nbasis);    /* basis function for this value of u */
/*
            printf("basis function for u = %f",u);
            for (i = 1; i <= npts; i++){
                printf("  %f",nbasis[i]);
            }
            printf("\n");
*/
            for (i = 1; i <= npts; i++){
                niku[i1] = nbasis[i];
                i1 = i1 + 1;
            }
            u = u + stepu;
        }
/*
        printf("\nu-basis functions\n\n");
        printf("I1,p1,npts = %d %d %d \n",i1-1, p1, npts);
 
        for (i = 1; i <= p1; i++){
            for (j = 1; j <= npts; j++){
                printf(" %f ",niku[(i-1)*npts + j]);
            }
            printf("\n");
        }
        printf("\n");
*/
 
/*
        getchar();
*/
 
/*  calculate the Mjl's at each w parametric value */
 
        w = 0.;
        for (winc = 1; winc <= p2; winc++){
            if ((float)y[mplusc] - w < 5e-6){
                w = (float)y[mplusc];
            }
            basis(l,w,mpts,y,mbasis);    /* basis function for this value of w */
/*
            printf("basis function for w = %f",w);
            for (i = 1; i <= mpts; i++){
                printf("  %f",mbasis[i]);
            }
            printf("\n");
*/
            for (j = 1; j <= mpts; j++){
                mjlw[i2] = mbasis[j];
                i2 = i2 + 1;
            }
            w = w + stepw;
        }
/*
        printf("\nw-basis functions\n\n");
        printf("I2 = %d \n",i2-1);
 
        for (i = 1; i <= p2; i++){
            for (j = 1; j <= mpts; j++){
            printf(" %f ",mjlw[(i-1)*mpts + j]);
            }
            printf("\n");
        }
        printf("\n");
*/
 
/*
        getchar();
*/
 
/*  calculate the sum function at each parametric value */
 
        for (uinc = 1; uinc <= p1; uinc++){
            for (i = 1; i <= npts; i++){
                ibas = (uinc-1)*npts + i;
                nbasis[i] = niku[ibas];
            }
            for (winc = 1; winc <= p2; winc++){
                for (j = 1; j <= p2; j++){
                    jbas = (winc-1)*mpts + j;
                    mbasis[j] = mjlw[jbas];
                }
                sum = sumrbas(b,nbasis,mbasis,npts,mpts);
/*              printf("returned from sumrbas. uinc, winc, u,w, sum = %d %d %f %f %f \n",uinc, winc, u,w,sum);*/
                if (sum is 0.){
                    printf("Error - Sum of the basis functions = 0 This is most likely caused\n"); 
                    printf("by all zero homogeneous weighting factors. Program stops.\n"); 
                    exit(0);
                }
                rsumij[i3] = (1/sum);
                i3 = i3 + 1;
            }
        }
/*
        printf("\nsum function\n\n");
        printf("I3 = %d \n",i3-1);
 
        for (j = 1; j <= p2; j++){
            for (i = 1; i <= p1; i++){
                printf(" %f ",rsumij[(j-1)*p1 + i]);
            }
            printf("\n");
        }
        printf("\n\n");
*/
        *p_itest = 0;
        }
 
 
/*  generate the complete rational B-spline surface */
 
    if (*p_itest is 0){             /* is is equivalent to ==  see rbsurf.h */
/*
        printf("calculating complete surface\n\n");
*/
        for (i = 1; i <= 3*p1*p2; i++){
            q[i] = 0.;
        }
 
        icount = 1;
    
        for (uinc = 1; uinc <= p1; uinc++){
            for (winc = 1; winc <= p2; winc++){
                scount = (uinc-1)*p2 + winc;
                for (i = 1; i <= npts; i++){
                    jbas = 4*mpts*(i-1);
                    ninc = (uinc-1)*npts + i;
                    if (niku[ninc] != 0.){
                        for (j = 1; j <= mpts; j++){
                            j1 = jbas +4*(j-1) + 1;
                            minc = (winc-1)*mpts + j;
                            if (mjlw[minc] !=0.){
                                pbasis = b[j1+3]*niku[ninc]*mjlw[minc]*rsumij[scount];
                                q[icount] = q[icount]+b[j1]*pbasis;  /* calculate surface point */
                                q[icount+1] = q[icount+1]+b[j1+1]*pbasis;
                                q[icount+2] = q[icount+2]+b[j1+2]*pbasis;
/*
                                printf("uinc,winc,i,jbas,ninc,j,j1,minc,scount\n");
                                printf("%d %d %d %d %d %d %d %d %d \n",uinc,winc,i,jbas,ninc,j,j1,minc,scount);
                                printf("niku[ninc],mjlw[minc],rsumij[scount],pbasis\n");
                                printf("%f %f %f %f \n\n",niku[ninc],mjlw[minc],rsumij[scount],pbasis);
                                printf("building surface %f %f %f \n\n",q[icount],q[icount+1],q[icount+2]);
*/                                                              
                            }
                        }
                    }
                }
/*              printf("surface point %d %d %f %f %f \n\n",uinc,winc,q[icount],q[icount+1],q[icount+2]);*/
                icount = icount + 3;
             }
        }
        *p_itest = nplusc + mplusc + p1 + p2;
        return;
    }
 
 
/*  calculate the incremental change to the surface */
 
        bx = b[ibnum] - bold[1];
        by = b[ibnum+1] - bold[2];
        bz = b[ibnum+2] - bold[3];
        bh = b[ibnum+3] - bold[4];
 
    if ((bx != 0.) || (by != 0.) || (bz != 0.) || (bh != 0.)){
 
/*  if true surface unchanged -- no calculation is needed
    if not true surface has changed -- do incremental calculation */
 
 
/*  printf("calculating incemental change -- bold,bx, by, bz, bh = \n %f %f %f %f %f %f %f %f \n",bold[1],bold[2],bold[3],bold[4],bx,by,bz,bh);*/
 
 
/*  calculate the i,j index for Bij from ibnum, where ibnum is
    assumed to be the lineal number of the x-component of Bij      */
 
        iindex = (ibnum/4/mpts) + 1;     /* depends on integer arithmatic to work */
        jindex=(ibnum-4*mpts*(iindex-1))/4 + 1;
 
/*  printf("ibnum,iindex,jindex %d %d %d \n",ibnum,iindex,jindex);*/
 
/*  for the special case of the homogeneous weighting factor changing
    the sum function must be recalculated for each value of u,w        */
 
    if (bh != 0.){
/*
        printf("save the old sum function \n");
*/
        for (i1 = 1; i1 <= p1*p2; i1++){
            savrsumij[i1] = rsumij[i1];
        }
/*
        printf("\nthe save sum function is \n\n");
 
        for (j = 1; j <= p2; j++){
                        for (i = 1; i <= p1; i++){
                                printf(" %f ",savrsumij[(j-1)*p1+i]);
                        }
                        printf("\n");
                        }
                printf("\n");
*/
 
/*
        printf("calculating new sum function \n");
*/
                tcount = 1;
        for (uinc = 1; uinc <= p1; uinc++){
            ninc = (uinc -1)*npts + iindex;
            for (winc = 1; winc <= p2; winc++){
                minc = (winc-1)*mpts + jindex;
/*
                printf("tcount, rsumij[tcount] = %d %f \n",tcount,rsumij[tcount]);
                printf("uinc,winc,ninc,niku[ninc],minc,mjlw[minc] = \n %d %d %d %f %d %f \n",uinc,winc,ninc,niku[ninc],minc,mjlw[minc]);
*/
                if ((niku[ninc] != 0.) && (mjlw[minc] != 0.)){
                    sumold = 1./(rsumij[tcount]);
                    sumnew = sumold + niku[ninc]*mjlw[minc]*(b[ibnum+3]-bold[4]);
                    rsumij[tcount] = 1./sumnew;
                }
/*
                printf("tcount,sumold,sumnew,rsumij = %d %f %f %f \n", tcount,sumold,sumnew,rsumij[tcount]);
*/
                tcount = tcount + 1;
            }
        }
    }
/*
    printf("\nsum function after incremental change \n\n");
 
    for (j = 1; j <= p2; j++){
        for (i = 1; i <= p1; i++){
            printf(" %f ",rsumij[(j-1)*p1+i]);
        }
        printf("\n");
    }
    printf("\n");
*/
 
/*  calculate the change in the surface for each u,w */
/* 
        printf("incremental surface calculation \n\n");
*/ 
        icount = 1;
        for (uinc = 1; uinc <= p1; uinc++){
            ninc = (uinc -1)*npts + iindex;
            if (niku[ninc] != 0.){
                for (winc = 1; winc <= p2; winc++){
                    minc = (winc-1)*mpts + jindex;
                    if (mjlw[minc] !=0.){
                        scount = (uinc-1)*p2 + winc;
                        if (bh is 0.){
/*                          printf("space coordinate changed \n\n");*/
                            pbasis = b[ibnum+3]*niku[ninc]*mjlw[minc]*rsumij[scount];
/*
                            printf("uinc,winc,ninc,niku[ninc],minc,mjlw[minc] = %d %d %d %f %d %f \n",uinc,winc,ninc,niku[ninc],minc,mjlw[minc]);
                            printf("scount,rsumij[scount],pbasis = %d %f %f \n",scount,rsumij[scount],pbasis);
*/                                                      
                            q[icount] = q[icount]+bx*pbasis;  /* calculate surface point */
                            q[icount+1] = q[icount+1]+by*pbasis;
                            q[icount+2] = q[icount+2]+bz*pbasis;
/*
                            printf("icount, qx qy qz = %d %f %f %f \n",icount,q[icount],q[icount+1],q[icount+2]);
*/                                                      
                        }
                        else{
/*
                            printf("homogeneous coordinate changed \n\n");
*/
                            pbasis = niku[ninc]*mjlw[minc]*rsumij[scount];
                            sumratio = rsumij[scount]/savrsumij[scount];
/*
                            printf("uinc,winc,ninc,niku[ninc],minc,mjlw[minc] = %d %d %d %f %d %f \n",uinc,winc,ninc,niku[ninc],minc,mjlw[minc]);
                            printf("scount,savrsumij[scount],rsumij[scount],pbasis = %d %f %f %f \n",scount,savrsumij[scount],rsumij[scount],pbasis);
*/
                            q[icount] = q[icount]*sumratio+bh*b[ibnum]*pbasis;  /* calculate surface point */
                            q[icount+1] = q[icount+1]*sumratio+bh*b[ibnum+1]*pbasis;
                            q[icount+2] = q[icount+2]*sumratio+bh*b[ibnum+2]*pbasis;
/*
                            printf("icount, qx qy qz = %d %f %f %f \n",icount,q[icount],q[icount+1],q[icount+2]);
*/
                        }
                    }
                    icount = icount + 3;
                }
            }
            else{
                icount = icount + 3*p2;
            }
        }
    }
    bold[1] = b[ibnum];
    bold[2] = b[ibnum+1];
    bold[3] = b[ibnum+2];
    bold[4] = b[ibnum+3];
}