NurbsSurface.cpp

/*
 * NurbsSurface.cpp
 *
 *  Created on: Mar 8, 2010
 *      Author: bruenger
 */


#include <iostream>
#ifdef BIAS_HAVE_XML2
#include <Base/Common/XMLIO.hh>
#endif
#include <Base/Math/Random.hh>

//FIXME: Why is this necessary: 
#include <Base/Common/XMLIO.hh>

#include <MathAlgo/NurbsSurface.hh>

using namespace std;
using namespace BIAS;

NurbsSurface::NurbsSurface() {
    initialized_ = false;
}

NurbsSurface::~NurbsSurface() {
    // TODO Auto-generated destructor stub
}

void NurbsSurface::Init(int dimu, int dimv, int degree)
{
    initialized_ = true;

    dimu_ = dimu;
    dimv_ = dimv;
    degree_ = degree;

    if(dimu_ < degree_ || dimv_ < degree_){
        cout << "NurbsSurface::Init - The Dimensions have to be greater then the degree" << endl;
        BIASABORT;
    }

    knotsu_.resize(dimu_ + degree_ + 1);
    knotsv_.resize(dimv_ + degree_ + 1);

    int i;

    ctrlpoints_.resize(dimu_);
    for (i = 0; i < dimu_; i++){
        ctrlpoints_[i].resize(dimv_);
    }

    SetKnotsToDefault();
    SetCTRLPointsToRandom(Vector2<double>(0.0,0.0),Vector2<double>(10.0,10.0));
}

void NurbsSurface::SetKnotsToDefault()
{
    if(!initialized_){
        cout << "NurbsSurface::SetKnotsToDefault - Please initialize the Surface first" << endl;
        BIASABORT;
    }
    double tmp;
    for(int i = 0; i < (dimu_ + degree_ + 1); i++){
        if(i < degree_ + 1){
            tmp = (double)0.0;
        }else if(i > (dimu_ - 1)){
            tmp = (double)1.0;
        }else{
            tmp = (double)(1.0/(dimu_ - degree_))*(double)(i-degree_);
        }
        SetKnotU(i,tmp);
    }

    for(int i = 0; i < (dimv_ + degree_ + 1); i++){
        if(i < degree_ + 1){
            tmp = (double)0.0;
        }else if(i > (dimv_ - 1)){
            tmp = (double)1.0;
        }else{
            tmp = (double)(1.0/(dimv_ - degree_))*(double)(i-degree_);
        }
        SetKnotV(i,tmp);
    }
}


void NurbsSurface::SetCTRLPointsToRandom(Vector2<double> start, Vector2<double> ende)
{
    if(!initialized_){
        cout << "NurbsSurface::SetCTRLPointsToRandom - Please initialize the Surface first" << endl;
        BIASABORT;
    }

    Random r;
    r.Reset();
    SetKnotsToDefault();

    double xdiff, ydiff, xstart, ystart;

    if (start.GetX() < ende.GetX()){
        xdiff = (ende.GetX() - start.GetX())/(dimu_ - 1);
        xstart = start.GetX();
    }else{
        xdiff = (start.GetX() - ende.GetX())/(dimu_ - 1);
        xstart = ende.GetX();
    }

    if (start.GetY() < ende.GetY()){
        ydiff = (ende.GetY() - start.GetY())/(dimv_ - 1);
        ystart = start.GetY();
    }else{
        ydiff = (start.GetY() - ende.GetY())/(dimv_ - 1);
        ystart = ende.GetY();
    }

    HomgPoint3D tmp;
    for(int x = 0; x < dimu_; x++){
        for(int y = 0; y < dimv_; y++){
            tmp.Set(xstart + (x*xdiff),ystart + (y*ydiff),(double)r.GetUniformDistributed(300.0, 500.0),1.0);
            SetCtrlPoint(x,y,tmp);
        }
    }

}

#ifdef BIAS_HAVE_XML2
int NurbsSurface::LoadFromXML(std::string xmlfile)
{
    if(initialized_){
        cout << "NurbsSurface::LoadFromXML - Your initialization will be overwritten by the XML" << endl;
    }

    if (xmlfile == "") {
        cout << "NurbsSurface::LoadFromXML - No XML-File given" << endl;
        return -1;
        //BIASABORT;
    }

    XMLIO myXML;
    xmlNodePtr rootNode;
    xmlNodePtr currentNode;
    xmlNodePtr ctrlpointsNode;
    xmlNodePtr currentCtrlPoint;

    rootNode = myXML.read(xmlfile);

    if (rootNode == NULL) {
        cout << "NurbsSurface::LoadFromXML - Error while parsing XML-File" << endl;
        BIASABORT;
    }

    currentNode = myXML.getChild(rootNode, "dimu");
    int dimu = myXML.getNodeContentInt(currentNode);
    currentNode = myXML.getChild(rootNode, "dimv");
    int dimv = myXML.getNodeContentInt(currentNode);
    currentNode = myXML.getChild(rootNode, "degree");
    int degree = myXML.getNodeContentInt(currentNode);


    Init(dimu,dimv,degree);


    ctrlpointsNode = myXML.getChild(rootNode, "ctrlpoints");
    currentCtrlPoint = myXML.getFirstChild(ctrlpointsNode);

    HomgPoint3D vec;

    for (int i = 0; i < dimu_; i++){
        for (int j = 0; j < dimv_; j++){
            vec.Set((double)myXML.getAttributeValueDouble(currentCtrlPoint,"x"),(double)myXML.getAttributeValueDouble(currentCtrlPoint,"y"),(double)myXML.getAttributeValueDouble(currentCtrlPoint,"z"), 1.0);
            SetCtrlPoint(i,j,vec);
            currentCtrlPoint = myXML.getNextChild();
        }
    }
  return 0;
}


void NurbsSurface::SaveToXML(std::string xmlfile)
{
  
  if(!initialized_){
    cout << "NurbsSurface::SaveToXML - Cannot save uninitialized NURBS surface" << endl;
  }
  
  if (xmlfile == "") {
    cout << "NurbsSurface::SaveToXML - No XML-File given" << endl;
    BIASABORT;
  }

  XMLIO myXML;
  xmlNodePtr rootNode;
  //xmlNodePtr currentNode;
  xmlNodePtr ctrlpointsNode;
  xmlNodePtr currentCtrlPoint;

  myXML.create("Surface");
  rootNode = myXML.getRootNode();
  myXML.addContent(myXML.addChildNode(rootNode,"dimu"),dimu_);
  myXML.addContent(myXML.addChildNode(rootNode,"dimv"),dimv_);
  myXML.addContent(myXML.addChildNode(rootNode,"degree"),degree_);
  //(myXML.getChild(rootNode,"dimu"))->addContent(dimu);
  //(myXML.getChild(rootNode,"dimv"))->addContent(dimv);
  //(myXML.getChild(rootNode,"degree"))->addContent(degree);
  
  myXML.addChildNode(rootNode, "ctrlpoints");

  ctrlpointsNode = myXML.getChild(rootNode, "ctrlpoints");
  //currentCtrlPoint = myXML.getFirstChild(ctrlpointsNode);
  

  for (int i = 0; i < dimu_; i++){
    for (int j = 0; j < dimv_; j++){
      currentCtrlPoint = myXML.addChildNode(ctrlpointsNode,"ctrlpoint");
      HomgPoint3D cPoint = GetCtrlPoint(i,j);
      myXML.addAttribute(currentCtrlPoint,"x",cPoint[0]);
      myXML.addAttribute(currentCtrlPoint,"y",cPoint[1]);
      myXML.addAttribute(currentCtrlPoint,"z",cPoint[2]);
      //vec.Set((double)myXML.getAttributeValueDouble(currentCtrlPoint,"x"),(double)myXML.getAttributeValueDouble(currentCtrlPoint,"y"),(double)myXML.getAttributeValueDouble(currentCtrlPoint,"z"), 1.0);
      //SetCtrlPoint(i,j,vec);
      //currentCtrlPoint = myXML.getNextChild();
    }
  }
  
  myXML.write(xmlfile);

}

void NurbsSurface::LoadFromXMLWeighted(std::string xmlfile,float weight){
  HomgPoint3D* oldCtrlPoints = new HomgPoint3D[dimu_*dimv_];
  for (int u = 0; u < dimu_; u++){
    for (int v = 0; v < dimv_; v++){
      oldCtrlPoints[u+v*dimu_] = GetCtrlPoint(u,v);
    }
  }
  LoadFromXML(xmlfile);
  
  for (int u = 0; u < dimu_; u++){
    for (int v = 0; v < dimv_; v++){
      HomgPoint3D pNew = GetCtrlPoint(u,v);
      pNew[0] = pNew[0] * (weight) + oldCtrlPoints[u+v*dimu_][0] * (1.0 - weight);
      pNew[1] = pNew[1] * (weight) + oldCtrlPoints[u+v*dimu_][1] * (1.0 - weight);
      pNew[2] = pNew[2] * (weight) + oldCtrlPoints[u+v*dimu_][2] * (1.0 - weight);
      SetCtrlPoint(u,v,pNew);
    }
  }  
    
    
}

#endif

int NurbsSurface::FindSpan(int n, int p, double u, std::vector<double> U){
    if(u == U[n+1]) return n;
    int low = p;
    int high = n+1;
    int mid = (low+high)/2;
    while(u < U[mid] || u >= U[mid+1]){
        if(u < U[mid])
            high = mid;
        else
            low = mid;
        mid= (low+high)/2;
    }
    return mid;
}

std::vector<double> NurbsSurface::BasisFuns(int i, double u, int p, std::vector<double> U){
    std::vector<double> N, left, right;
    N.resize(p+1);
    left.resize(p+1);
    right.resize(p+1);
    double temp, saved;
    N[0] = 1.0f;
    for(int j = 1; j <= p; j++){
        left[j] = u-U[i+1-j];
        right[j] = U[i+j]-u;
        saved = 0.0f;
        for(int r = 0; r < j; r++){
            temp = N[r]/(right[r+1]+left[j-r]);
            N[r] = saved+right[r+1]*temp;
            saved = left[j-r]*temp;
        }
        N[j] = saved;
    }

    return N;
}


Vector3<double> NurbsSurface::GetSurfacePoint(double u, double v){
    if(u == 1.0){
        u = 0.999f;
    }
    if(v == 1.0){
        v = 0.999f;
    }
    int uspan = FindSpan(dimu_,degree_,u,knotsu_);
    std::vector<double> Nu = BasisFuns(uspan,u,degree_,knotsu_);
    int vspan = FindSpan(dimv_,degree_,v,knotsv_);
    std::vector<double> Nv = BasisFuns(vspan,v,degree_,knotsv_);
    int uind = uspan - degree_;
    Vector3<double> temp, S;
    S.SetZero();
    
    for(int l = 0; l <= degree_; l++){
        temp.SetZero();
        int vind = vspan-degree_+l;
        for(int k = 0; k <= degree_; k++){
          double* ptr =  GetCtrlPoint(uind+k,vind).GetData();
            temp[0] += Nu[k] * ptr[0];
            temp[1] += Nu[k] * ptr[1];
            temp[2] += Nu[k] * ptr[2];
        }
        S.Set(S.GetData()[0] + Nv[l]*temp[0],S.GetData()[1] + Nv[l]*temp[1],S.GetData()[2] + Nv[l]*temp[2]);
    }
    return S;
}