ContourBSplineData.cpp

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

Copyright (C) 2003-2009    (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 "ContourBSplineData.hh"

using namespace BIAS;
using namespace std;
//public:

ContourBSplineData* ContourBSplineData::Register(const unsigned int order,
                                       const ContourBSplineType::Type bType,
                                       const std::vector<unsigned int>& mPnts){
  ContourBSplineData* res; //return value

  //look if already registered
  res=LookUpData_(order, bType, mPnts);
  cout<<"Lookup data:"<<endl;
  if(res==NULL){//not registered
    res = new ContourBSplineData;
    res->order_=order;
    res->bType_=bType;
    res->mPnts_=mPnts;
    res->numBasePolynoms_=mPnts.size();
    cout<<"before RegisterPeriodicBSpline"<<endl;
    //compute numBasePolynoms,numSpans, placedSpanMatrices for BSpline
    if(bType==ContourBSplineType::Open){
      res->RegisterAperiodicBSpline_();
    }
    else{ //bType==Closed
      res->RegisterPeriodicBSpline_();
    }
    res->numSpansVec_.push_back(res->numSpans_);
    cout<<"RegisterPeriodicBSpline"<<endl;
    //compute data matrices
    res->ComputeSplineMetricMatrix_();
    res->ComputeSplineMetricMatrixParametric_();
    res->ComputeAreaCoefficientMatrix_();

    //save pointer and init reference counter for pointer to 1
    (res->reference_)=1;
    pntr_.push_back(res);
    cout<<"ComputeAreaCoefficientMatrix:"<<endl;
  }else{

   cout<<"already registerd:"<<endl;
  }
  return res;
}


ContourBSplineData* ContourBSplineData::Register(
                         const std::vector<ContourBSplineData*>& bSplinesData){
  unsigned int i,j,k,l; //counter variables

  ContourBSplineData* res;
  /**
   * no lookup here ... cause composition of periodic and aperiodic bsplines
   * possible
   */
  res = new ContourBSplineData;

  //test order for all splines (must be same at the moment)
  for (i=1; i<bSplinesData.size(); i++){
    if ((bSplinesData[i])->order_!= (bSplinesData[0])->order_ ){
      BIASERR("Clustered ContourBSplines have to be of same order.\r\n");
    }
  }
  res->order_=(bSplinesData[0])->order_;
  /**
   * type Cluster isnt meaningful at all ... clustered splines can be
   * composition of periodic and aperiodic bsplines
   */
  res->bType_=ContourBSplineType::Cluster;

  std::vector<unsigned int> mPnts;
  std::vector<BIAS::Matrix<double> > placedSpanMatrices;
  unsigned int numBasePolynoms=0;
  unsigned int numSpans=0;

  //compute new mPnt vector, numSpansVec, numBasePolynoms, numSpans
  for(i=0; i<bSplinesData.size(); i++){
    for(j=0; j<((bSplinesData[i])->mPnts_).size(); j++){
      mPnts.push_back(((bSplinesData[i])->mPnts_)[j]);
    }
    numBasePolynoms+=(bSplinesData[i])->numBasePolynoms_;
    numSpans+=(bSplinesData[i])->numSpans_;
    res->numSpansVec_.push_back((bSplinesData[i])->numSpans_);
  }
  res->mPnts_=mPnts;
  res->numBasePolynoms_=numBasePolynoms;
  res->numSpans_=numSpans;

  //correct/build placedSpanMatrices
  unsigned int numBasePolynomsTmp = 0;
  for(i=0;i<bSplinesData.size(); i++){
    for(j=0;j<((bSplinesData[i])->placedSpanMatrices_).size();j++){
      BIAS::Matrix<double> newPlacedSpan(res->order_,res->numBasePolynoms_);
      newPlacedSpan.SetZero();
      for(k=0;k<res->order_;k++){
        for(l=numBasePolynomsTmp;
            l<((unsigned int)(((bSplinesData[i])->placedSpanMatrices_)[j]
                ).num_cols()+numBasePolynomsTmp);l++){
          newPlacedSpan[k][l]=((bSplinesData[i])->placedSpanMatrices_
            )[j][k][l-numBasePolynomsTmp];
        }
      }
      placedSpanMatrices.push_back(newPlacedSpan);
    }
    numBasePolynomsTmp+=(bSplinesData[i])->numBasePolynoms_;
  }

  res->placedSpanMatrices_=placedSpanMatrices;

  res->ComputeSplineMetricMatrix_();
  res->ComputeSplineMetricMatrixParametric_();
  res->ComputeAreaCoefficientMatrix_();

  //save pointer and init reference counter for pointer to 1
  (res->reference_)=1;
  pntr_.push_back(res);

  return res;
}

ContourBSplineData* ContourBSplineData::Register(const unsigned int order,
                  const ContourBSplineType::Type bType,
                  const std::vector<unsigned int>& mPnts,
                  const unsigned int numBasePolynoms,
                  const unsigned int numSpans,
                  const std::vector<unsigned int>& numSpansVec,
                  const std::vector<BIAS::Matrix<double> >&placedSpanMatrices){
  ContourBSplineData* res;
  res = new ContourBSplineData;

  res->order_ = order;
  res->bType_ = bType;
  res->mPnts_ = mPnts;
  res->numBasePolynoms_ = numBasePolynoms;
  res->numSpans_ = numSpans;
  res->numSpansVec_ = numSpansVec;
  res->placedSpanMatrices_ = placedSpanMatrices;
  res->ComputeSplineMetricMatrix_();
  res->ComputeSplineMetricMatrixParametric_();
  res->ComputeAreaCoefficientMatrix_();

  //save pointer and init reference counter for pointer to 1
  (res->reference_)=1;
  pntr_.push_back(res);

  return res;
}

//protected:

std::list<ContourBSplineData*> ContourBSplineData::pntr_;

ContourBSplineData* ContourBSplineData::LookUpData_(const unsigned int order,
                                      const ContourBSplineType::Type bType,
                                      const std::vector<unsigned int>& mPnts){
  ContourBSplineData* res=NULL; //default return value

  std::list<ContourBSplineData*>::iterator it;
  for(it=pntr_.begin(); it!=pntr_.end() ;it++){
    if( ((*it)->order_==order) && ((*it)->bType_==bType) &&
        ((*it)->mPnts_==mPnts) ){
      ((*it)->reference_)++; //increase reference counter
      return *it;
    }
  }
  return res;
}

void ContourBSplineData::RegisterPeriodicBSpline_(){
  unsigned int i; //counter variable

  //compute numSpans
  unsigned int multipleKnotCount=0;
  for(i=0;i<numBasePolynoms_;i++){
    multipleKnotCount+=mPnts_[i]-1;
  }
  numSpans_ = numBasePolynoms_ - multipleKnotCount;

  /* initialise periodic case
   * - is treated like aperiodic case; therefore virtual spans are used
   *   - 2x1 virtual termination spans at front and back
   *   - 2x1 virtual spans before and after real spans
   */

  /* build a temporary vector of knot multiplicities for bspline
   * -  vector has size numspans_+1
   *
   * - loop backwards! through mPnts_
   * - for first multiple control point look if it causes an offset to later
   *   bIndex (periodicOffset)
   * - for all other multiple control point compute index to vector
   */
  std::vector<unsigned int> kMultTmp(numSpans_+1,1);
  //initial offset var for building kMultTmp
  unsigned int j=multipleKnotCount-(order_-2);
  unsigned int periodicOffset=0;
  for(i=mPnts_.size();i>0;i--){
    if(mPnts_[i-1]>1){
      if( (mPnts_.size()-(i-1)) < (order_-1) ){
        /* nonrecurring case
         * if this case occurs, then only for the first multiple control point
         * (counted backward)
         */
        periodicOffset=(i-1)+(order_-1)-mPnts_.size();
      }
      kMultTmp[(i-1)-j-periodicOffset]=order_-1;//restricted to order-1
                                               //should be mPnts_[i]
      j-=order_-2;
    }
  }
  //multiplicity of last knot must be same as first knot
  kMultTmp.back()=kMultTmp.front();

  //correct placement matrix if first knot has simple multiplicity
  if(kMultTmp[0]==1){
    periodicOffset-=order_-2;
  }

  cout<<"In RegisterPeriodicBSpline 1"<<endl;

  /*
   * - build up vector of knot multiplicities used to build spans of the
   *   bspline
   * - includes virtual spans
   */
  std::vector<unsigned int> kMult;
  /* - first and last entry is order_, to share build code with aperiodic
   *   bsplines (knot multiplicity of order_ causes termination)
   * - actually the first and last spans are not computed
   * - no need to decrease periodicOffset here (ComputeSpan_(...) cares)
   */
  kMult.push_back(order_);
  /* - add one more knot to vector which causes one more virtual span
   * - this knot has the same multiplicity as the beginning knot of the last
   *   real span
   * - decrease periodicOffset
   */
  unsigned int multiplicity = kMultTmp[kMultTmp.size()-2];
  kMult.push_back(multiplicity);
  periodicOffset-=multiplicity;

  cout<<"In RegisterPeriodicBSpline 2"<<endl;
  /* - fill kMult_ with kMultTmp (causes the real spans)
   * - decrease periodicOffset by first entry of kMultTmp (2 share code with
   *   aperiodic case)
   */
  for(i=0;i<kMultTmp.size();i++){
    kMult.push_back(kMultTmp[i]);
  }
  periodicOffset-=kMultTmp[0];
  /* - add one more knot to vector which causes one more virtual span
   * - this knot hast the same multiplicity as the end knot of the first real
   *   span
   * - no need to decrease periodicOffset because this virtual span is behind
   *   the real spans
   */
  kMult.push_back(kMultTmp[1]);
  // termintion knot (virtual span)
  kMult.push_back(order_);

  //build knots (breakpoints) (including virtual spans)
  std::vector<int> k;
  int q=0;
  for(i=0;i<kMult.size();i++){
    for(j=0;j<kMult[i];j++){
      k.push_back(q);
    }
    q++;
  }
  cout<<"In RegisterPeriodicBSpline 3"<<endl;
  //compute the spans (ignore virtual spans)
  for(i=2;i<numSpans_+2;i++){
    ComputeSpan_(i,periodicOffset,kMult,k);
  }
  cout<<"In RegisterPeriodicBSpline 4"<<endl;
}

void ContourBSplineData::RegisterAperiodicBSpline_(){
  unsigned int i; //counter variable

  //compute numSpans
  unsigned int multipleKnotCount=0;
  for(i=1;i<numBasePolynoms_-1;i++){
    multipleKnotCount+=mPnts_[i]-1;
  }
  numSpans_ = numBasePolynoms_ - order_+1 - multipleKnotCount;

  /*
   * the construction of the aperiodic case is treated like in the book
   * "active contours"
   */

  /*
   * build vector of knot multiplicities from control point multiplicities
   * - loop over control points (first possible multiple control is the
   *   order_th controlpoint (index: order_-1)
   * - first and last knot are termination knots and thus have multiplicity
   *   order
   */
  std::vector<unsigned int> kMult(numSpans_+1,1);
  kMult.front()=order_;
  unsigned int j=order_-2;//initial offset need to build kMult
  for(i=order_-1;i<numBasePolynoms_-1;i++){
    if(mPnts_[i]>1){
      kMult[i-j]=mPnts_[i];//should possibly restricted to order_-1
      j+=order_-2;
    }
  }
  kMult.back()=order_;

  //build knots (breakpoints)
  std::vector<int> k;
  int q=0;
  for(i=0;i<=numSpans_;i++){
    for(j=0;j<kMult[i];j++){
      k.push_back(q);
    }
    q++;
  }


  //compute spans
  for(i=0;i<numSpans_;i++){
    ComputeSpan_(i, 0, kMult, k);
  }
}

BIAS::Polynom ContourBSplineData::ComputeBasePolynomRecursive_(const unsigned int span,
                                                              const unsigned int knotIndex,
                                                              const unsigned int order,
                                                              const std::vector<int>& breakpoints){
  //ground instance of the recursive algorithm
  if (order==1){
    BIAS::Polynom p(1);
    p[0]=(span>=(unsigned int)breakpoints[knotIndex] &&
          span<(unsigned int)breakpoints[knotIndex+1]) ? 1. : 0.;
    return p;
  }
  //recursive rule
  else{
    int denominator = breakpoints[knotIndex+order-1] - breakpoints[knotIndex];
    BIAS::Polynom p1(0);
    if(denominator){
      BIAS::Polynom pTmp(0), pTmp2(0);
      pTmp.SetOrder(2);
      pTmp[1]=1.; pTmp[0]=(double)span-(double)breakpoints[knotIndex];
      //BIAS::Polynom pTmp2(0);
      pTmp2 = ComputeBasePolynomRecursive_(span, knotIndex, order-1, breakpoints);
     // cout << pTmp2 << endl;
     // cout << pTmp << endl;
      p1 = (pTmp*pTmp2);// * (1./(double)denominator);
      //cout << p1 << endl;
      p1 *= (1./(double)denominator);
    }
    denominator = breakpoints[knotIndex+order] - breakpoints[knotIndex+1];
    BIAS::Polynom p2(1);p2[0]=0.;
    if(denominator){
      BIAS::Polynom pTmp(2);
      pTmp.SetOrder(2);
      pTmp[1]=-1.; pTmp[0]=(double)breakpoints[knotIndex+order]-(double)span;
      BIAS::Polynom pTmp2(0);
      pTmp2 = ComputeBasePolynomRecursive_(span, knotIndex+1, order-1, breakpoints);
      p2 = pTmp*pTmp2;
      p2 = p2 * (1./(double)denominator);
    }
    BIAS::Polynom result(0);
    result=p1+p2;
    return result;
  }
}


void ContourBSplineData::ComputeSpan_(const unsigned int span,
                                     const int placeOff,
                                     const std::vector<unsigned int>& kMult,
                                     const std::vector<int>& k){
  unsigned int i,j; //loop variables
  cout<<"In ComputeSpan 1"<<endl;
  // compute index of first basepolynom which is non-zero for the span
  int bIndex=0;
  for(i=0; i<=span;i++){
    bIndex+=(int) kMult[i];
  }
  bIndex-=(int)order_;
   cout<<"In ComputeSpan 2"<<endl;
  /*
   * compute span and placement matrix
   */
  BIAS::Matrix<double> spanMatrix(order_, order_);
  BIAS::Matrix<double> placementMatrix(order_, numBasePolynoms_);
  BIAS::Polynom basePolynom;
  for(i=0;i<order_;i++){
    cout<<"In ComputeSpan 2."<<i<<endl;
    basePolynom=ComputeBasePolynomRecursive_(span, (unsigned int)bIndex+i,
                                            order_, k);
    cout<<"In ComputeSpan 2."<<i<<endl;
    //fill jth columns on current span matrix with coefficients of basePolynom
    for(j=0;j<order_;j++){
      spanMatrix[j][i]=basePolynom[j];
    }

    //compute placementMatrix (different for aperiodic and periodic case)
    for(j=0;j<numBasePolynoms_;j++){
      //aperiodic case
      if(bType_==ContourBSplineType::Open){
        placementMatrix[i][j]= i+bIndex==j ? 1.:0.;
      }
      //periodic case
      else
      {
        //in periodic case the placementMatrix needs extra placement offset
        placementMatrix[i][j]=((int)i+bIndex+placeOff+(int)numBasePolynoms_)
                              % numBasePolynoms_ == j ? 1.:0.;
      }
    }
  }

  cout<<"In ComputeSpan 3"<<endl;
  placementMatrix.MultLeft(spanMatrix); //inplace mult!
  placedSpanMatrices_.push_back(placementMatrix);
}

void ContourBSplineData::ComputeSplineMetricMatrix_(){
  //see p.292 for details
  unsigned int i,j; //counter variables

  //compute Hilbert matrix
  BIAS::Matrix<double> hilbert(order_,order_);
  for(i=1;i<=order_;i++)
    for(j=1;j<=order_;j++)
      hilbert[i-1][j-1]=1./(double)(i+j-1);

  /* formula:
   * BSplineMetricMatrix = 1/numSpans_ * sum_(i=0)^(numSpans_-1)(
   * placementMatrices_[i]^T * spanMatrices_[i]^T * hilbert * spanMatrices_[i]
   * * placementMatrices_[i]
   * )
   */
  splineMetricMatrix_.newsize(numBasePolynoms_,numBasePolynoms_);
  splineMetricMatrix_.SetZero();
  BIAS::Matrix<double> tmpMat(numBasePolynoms_,numBasePolynoms_);
  for(i=0;i<numSpans_;i++){
    tmpMat.SetIdentity();
    tmpMat.MultLeft(placedSpanMatrices_[i]);
    tmpMat.MultLeft(hilbert);
    tmpMat.MultLeft(placedSpanMatrices_[i].Transpose());
    splineMetricMatrix_+=tmpMat;
  }
  splineMetricMatrix_.MultiplyIP(1./(double)numSpans_);
}

void ContourBSplineData::ComputeSplineMetricMatrixParametric_(){
  /*
   * converts the metric matrix of the bspline to a parametric version
   *  U = (B 0)
   *      (0 B)
   */
  unsigned int i,j; //loop variables
  splineMetricMatrixParametric_.newsize(2*numBasePolynoms_,2*numBasePolynoms_);
  splineMetricMatrixParametric_.SetZero();

  //fill upper left
  for(i=0;i<numBasePolynoms_;i++){
    for(j=0;j<numBasePolynoms_;j++){
      splineMetricMatrixParametric_[i][j]=splineMetricMatrix_[i][j];
    }
  }
  //fill lower right
  for(i=numBasePolynoms_;i<2*numBasePolynoms_;i++){
    for(j=numBasePolynoms_;j<2*numBasePolynoms_;j++){
      splineMetricMatrixParametric_[i][j]=
         splineMetricMatrix_[i-numBasePolynoms_][j-numBasePolynoms_];
    }
  }
}



void ContourBSplineData::ComputeAreaCoefficientMatrix_(){
  //see p.293 for details
  unsigned int i,j; //counter variables

  //compute Matrix P'
  BIAS::Matrix<double> pQuote(order_,order_);

  for(i=1;i<=order_;i++)
    for(j=1;j<=order_;j++)
      pQuote[i-1][j-1] = ((i==1) && (j==1))? 0. : (double)(j-1)/(double)(i+j-2);

  /*
   * formula bQuote:
   * bQuote=1/numSpans * sum_(i=0)^(numSpans-1)(
   * placedSpanMatrices[i]^T  * pQuote * placedSpanMatrices[i]
   * )
   */
  BIAS::Matrix<double>bQuote(numBasePolynoms_,numBasePolynoms_);
  bQuote.SetZero();
  BIAS::Matrix<double>tmpMat;
  for(i=0;i<numSpans_;i++){
    tmpMat=placedSpanMatrices_[i];
    tmpMat.MultLeft(pQuote);
    tmpMat.MultLeft(placedSpanMatrices_[i].Transpose());
    bQuote+=tmpMat;
  }
  bQuote.MultiplyIP(1./(double)numSpans_);

  //finaly compute area coefficient matrix
  areaCoefficientMatrix_.newsize(2*numBasePolynoms_,2*numBasePolynoms_);
  areaCoefficientMatrix_.SetZero();

  for(i=0;i<numBasePolynoms_;i++){
    for(j=0;j<numBasePolynoms_;j++){
      //fill upper left
      areaCoefficientMatrix_[i][j]=bQuote[i][j];
      //fill lower right
      areaCoefficientMatrix_[i+numBasePolynoms_][j+numBasePolynoms_]=-bQuote[i][j];
    }
  }
}