Nurbs3D.cpp

/*
 * Nurbs3D.cpp
 *
 *  Created on: Jul 1, 2010
 *      Author: jordt
 */

#include "Nurbs3D.hh"
#include <Base/Math/Vector2.hh>
#include <Base/Math/Vector4.hh>
#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Geometry/HomgPoint3D.hh>
#include <Base/Math/Operators.hh>

#include "NurbsTemplateInit.hh"

using namespace std;
using namespace BIAS;


// This function allocates and initializes the control point- and weight-arrays and
// calculates the sphere-points (control points in the vicinity of the unit sphere)
template <class ControlPointsT>
Nurbs3D<ControlPointsT>::Nurbs3D(int knotsPerDim, int degree ) : NurbsBase(knotsPerDim,degree,3){

  // Allocate and initialize
  controlPoints_ = new ControlPointsT**[numCtrlPoints_];
  weights_ = new double**[numCtrlPoints_];
  for (int x = 0; x < numCtrlPoints_; x++){
    controlPoints_[x] = new ControlPointsT*[numCtrlPoints_];
    weights_[x] = new double*[numCtrlPoints_];
    for (int y = 0; y < numCtrlPoints_; y++){
      controlPoints_[x][y] = new ControlPointsT[numCtrlPoints_];
      weights_[x][y] = new double[numCtrlPoints_];
      for (int z = 0; z < numCtrlPoints_; z++){
        TemplateInit::Init(controlPoints_[x][y][z]);
        weights_[x][y][z] = 1.0;

        // Calculate the sphere points
        double xSphere = ((double)x/((double)numCtrlPoints_-1.0) - 0.5) * 2.0;
        double ySphere = ((double)y/((double)numCtrlPoints_-1.0) - 0.5) * 2.0;
        double zSphere = ((double)z/((double)numCtrlPoints_-1.0) - 0.5) * 2.0;
        Vector3<double> prm = Vector3<double>(xSphere,ySphere,zSphere);
        if (prm.Length() > 0.0001) prm.Normalize();
        prm[0] = (prm[0] + 1.0) / 2.0;
        prm[1] = (prm[1] + 1.0) / 2.0;
        prm[2] = (prm[2] + 1.0) / 2.0;
        double testResult;
        EvalBaseFunc(x,y,z,prm,testResult);
        if (testResult > 0.0) sphereCtrlPoints_.push_back(Vector3<int>(x,y,z));

        //if (fabs(1.0 - sqrt(xSphere*xSphere+ySphere*ySphere+zSphere*zSphere)) <= (double)degree/(double)(numCtrlPoints_-1) )
        
      }
    }
  }
  
}

template <class ControlPointsT>
Nurbs3D<ControlPointsT>::~Nurbs3D(){
  for (int x = 0; x < numCtrlPoints_; x++){
    for (int y = 0; y < numCtrlPoints_; y++){
      delete[] controlPoints_[x][y];
      delete[] weights_[x][y];
    }
    delete[] controlPoints_[x];
    delete[] weights_[x];
  }
  delete[] controlPoints_;
  delete[] weights_;


}


template <class ControlPointsT>
void Nurbs3D<ControlPointsT>::Eval(const Vector3<double> &nurbsParam, ControlPointsT &result) const{

  BIASASSERT(nurbsParam[0] >= 0.0 && nurbsParam[0] <= 1.0 &&
             nurbsParam[1] >= 0.0 && nurbsParam[1] <= 1.0 &&
             nurbsParam[2] >= 0.0 && nurbsParam[2] <= 1.0)

  // Initialize the return value
  TemplateInit::Init(result);

  // Get the relevant control points
  double weightSum = 0.0;
  Vector3<int> relevantBegin;
  Vector3<int> relevantEnd;
  GetRelevantBaseCtrlPoints(nurbsParam, relevantBegin, relevantEnd);
  double* baseFuncCache = new double[degree_+1];
  
  // Evaluate the relevant base functions and use them to wheight the
  // control points
  for (int x = relevantBegin[0]; x <= relevantEnd[0]; x++)
    for (int y = relevantBegin[1]; y <= relevantEnd[1]; y++)
      for (int z = relevantBegin[2]; z <= relevantEnd[2]; z++){
        double weight = (BaseFunc(0,x,degree_,nurbsParam[0],baseFuncCache)*
                   BaseFunc(1,y,degree_,nurbsParam[1],baseFuncCache)*
                   BaseFunc(2,z,degree_,nurbsParam[2],baseFuncCache))*
                   weights_[x][y][z];
        weightSum += weight;
  
        result += (weight * controlPoints_[x][y][z]);
      }

  if (weightSum != 0.0)
    result = result * (1.0 / weightSum);

  delete[] baseFuncCache;
}

template <class ControlPointsT>
void Nurbs3D<ControlPointsT>::Init(const ControlPointsT minX,const ControlPointsT maxX,
                                   const ControlPointsT minY,const ControlPointsT maxY,
                                   const ControlPointsT minZ,const ControlPointsT maxZ){
  for (int x = 0; x < numCtrlPoints_; x++)
    for (int y = 0; y < numCtrlPoints_; y++)
      for (int z = 0; z < numCtrlPoints_; z++){
        double ratioX = (double)x/(double)(numCtrlPoints_-1);
        double ratioY = (double)y/(double)(numCtrlPoints_-1);
        double ratioZ = (double)z/(double)(numCtrlPoints_-1);
        controlPoints_[x][y][z] = ratioX * maxX + (1.0 - ratioX) * minX
                                + ratioY * maxY + (1.0 - ratioY) * minY
                                + ratioZ * maxZ + (1.0 - ratioZ) * minZ;
      }

}

namespace BIAS {

  // Implement NURBS3D for all common types of vectors:
  template class Nurbs3D<double >;
  template class Nurbs3D<Vector<double> >;
  template class Nurbs3D<Vector2<double> >;
  template class Nurbs3D<Vector3<double> >;
  template class Nurbs3D<Vector4<double> >;
  template class Nurbs3D<HomgPoint2D >;
  template class Nurbs3D<HomgPoint3D >;

}