BSplineCurve.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009, 2005    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#include "BSplineCurve.hh"
#include <Base/Common/BIASpragma.hh>

using namespace BIAS;
using namespace std;

BSplineCurve::BSplineCurve(){
#ifdef BIAS_DEBUG
  initializedQuadratic_=false;
  initializedCubic_=false; 
#endif
}

void BSplineCurve::InitQuadraticUniformBSpline(){
#ifdef BIAS_DEBUG
  if(cPnts_.size()<3){
    if(cPnts_.size()==0){
      BIASERR("BSpline: no control points given - call ");
      BIASERR("SetControlPoints(cPnts)!\r\n");
      exit(-1);
    }
    else{
      BIASERR("BSpline: less then 3 control points given - call ");
      BIASERR("SetControlPoints(cPnts) with 3 or more control points!\r\n");
      exit(-1);
    }
  }
#endif
  Matrix<double> coeffM(3,3);
  Vector<double> xItem(3);
  Vector<double> yItem(3);
  Vector<double> xTmp(3);
  Vector<double> yTmp(3);


  coeffM[0][0]=1;
  coeffM[0][1]=-2;
  coeffM[0][2]=1;
  coeffM[1][0]=-2;
  coeffM[1][1]=2;
  coeffM[1][2]=0;
  coeffM[2][0]=1;
  coeffM[2][1]=1;
  coeffM[2][2]=0;

  xCoeff_.clear();
  yCoeff_.clear();

  //compute coeff matrix for all segments
  for(unsigned int i=0; i<cPnts_.size()-2; i++){
    for(unsigned int j=0; j<3;j++){
      xTmp[j]=cPnts_[i+j][0];
      yTmp[j]=cPnts_[i+j][1];
    }

    coeffM.Mult(xTmp, xItem);
    coeffM.Mult(yTmp, yItem);
    xCoeff_.push_back(xItem);
    yCoeff_.push_back(yItem);
  }
#ifdef BIAS_DEBUG
  initializedQuadratic_=true;
  initializedCubic_=false;
#endif
  
}

void BSplineCurve::InitCubicUniformBSpline(){
#ifdef BIAS_DEBUG
  if(cPnts_.size()<4){
    if(cPnts_.size()==0){
      BIASERR("BSpline: no control points given - call ");
      BIASERR("SetControlPoints(cPnts)!\r\n");
      exit(-1);
    }
    else{
      BIASERR("BSpline: less then 4 control points given - call ")
      BIASERR("SetControlPoints(cPnts) with 4 or more control points!\r\n");
      exit(-1);
    }
  }
#endif
  Matrix<double> coeffM(4,4);
  Vector<double> xItem(4);
  Vector<double> yItem(4);
  Vector<double> xTmp(4);
  Vector<double> yTmp(4);

  coeffM[0][0]=-1;
  coeffM[0][1]=3;
  coeffM[0][2]=-3;
  coeffM[0][3]=1;
  coeffM[1][0]=3;
  coeffM[1][1]=-6;
  coeffM[1][2]=3;
  coeffM[1][3]=0;
  coeffM[2][0]=-3;
  coeffM[2][1]=0;
  coeffM[2][2]=3;
  coeffM[2][3]=0;
  coeffM[3][0]=1;
  coeffM[3][1]=4;
  coeffM[3][2]=1;
  coeffM[3][3]=0;

  xCoeff_.clear();
  yCoeff_.clear();

  //compute coeffmatrix for all segments
  for(unsigned int i=0; i<cPnts_.size()-3;i++){
    for(unsigned int j=0; j<4; j++){
      xTmp[j]=cPnts_[i+j][0];
      yTmp[j]=cPnts_[i+j][1];
    }
    coeffM.Mult(xTmp,xItem);
    coeffM.Mult(yTmp,yItem);
    xCoeff_.push_back(xItem);
    yCoeff_.push_back(yItem);
  }
#ifdef BIAS_DEBUG
  initializedCubic_=true;
  initializedQuadratic_=false;
#endif
}

  
void BSplineCurve::QuadraticUniformBSpline(Vector2<double> &res, double t){
#ifdef BIAS_DEBUG
  if(!initializedQuadratic_){
    BIASERR("BSpline: call InitQuadraticUniformBSpline() at first!\r\n");
    exit(-1);
  }
#endif
  unsigned int actSegment;
  
  //compute spline segment and new t for the result segment
  if(t<=0.){
    actSegment=0;
    t=0.;
  } else if(t>=1.){
    actSegment=cPnts_.size()-3;
    t=1.;
  } else {
    actSegment= (unsigned int)(t * (double)(cPnts_.size()-2));
    t= t* (double)(cPnts_.size()-2) - (double) actSegment;
  }
  
  Vector<double> tVec(3);
  tVec[0]=t*t;
  tVec[1]=t;
  tVec[2]=1;

  res[0]=0.5 * tVec.ScalarProduct(xCoeff_[actSegment]);
  res[1]=0.5 * tVec.ScalarProduct(yCoeff_[actSegment]);
}

void BSplineCurve::CubicUniformBSpline(Vector2<double>  &res, double t){
#ifdef BIAS_DEBUG
  if(!initializedCubic_){
    BIASERR("BSpline: call InitCubicUniformBSpline() at first!\r\n");
      exit(-1);
  }
#endif
  unsigned int actSegment;

  //compute spline segment and new t for the result segment
  if(t<=0.){
    actSegment=0;
    t=0.;
  } else if(t>=1.){
    actSegment=cPnts_.size()-4;
    t=1.;
  } else {
    actSegment= (unsigned int) (t * (double)(cPnts_.size()-3));
    t= t * (double)(cPnts_.size()-3) - (double) actSegment;
  }

  Vector<double> tVec(4);
  tVec[0]=t*t*t;
  tVec[1]=t*t;
  tVec[2]=t;
  tVec[3]=1;
  
  res[0]=(1./6.) * tVec.ScalarProduct(xCoeff_[actSegment]);
  res[1]=(1./6.) * tVec.ScalarProduct(yCoeff_[actSegment]);
}