ContourBSplineShapeMatrix.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 <Base/Debug/Debug.hh>
#include <MathAlgo/SVD.hh>

#include "ContourBSplineShapeMatrix.hh"

using namespace BIAS;

//public:

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::SetShapeSpaceMatrix(
                                           const BIAS::Matrix<double>& matrix){
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  //copy matrix
  mat->newsize(matrix.num_rows(),matrix.num_cols());
  int i,j;
  for(i=0;i<matrix.num_rows();i++){
    for(j=0;j<matrix.num_cols();j++){
      (*mat)[i][j] = matrix[i][j];
    }
  }
  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::SetShapeSpaceIdentity(
                                           const unsigned int numBasePolynoms){
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  mat->newsize(2*numBasePolynoms,2*numBasePolynoms);
  mat->SetIdentity();
  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::SetShapeSpaceEuclidian(
                                            const unsigned int numBasePolynoms,
                                            const BIAS::Vector<double>& Q){
  unsigned int i; //counter variable;
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  mat->newsize(2*numBasePolynoms,4);
  mat->SetZero();
  
  for(i=0;i<numBasePolynoms;i++){
    (*mat)[i][0]=1.;
    (*mat)[i+numBasePolynoms][1]=1.;
    (*mat)[i][2]=Q[i];
    (*mat)[i+numBasePolynoms][2]=Q[i+numBasePolynoms];
    (*mat)[i][3]=-Q[i+numBasePolynoms];
    (*mat)[i+numBasePolynoms][3]=Q[i];
  }  
  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::SetShapeSpacePlanarAffin(
                                            const unsigned int numBasePolynoms,
                                            const BIAS::Vector<double>& Q){
  unsigned int i; //counter variable
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  mat->newsize(2*numBasePolynoms,6);
  mat->SetZero();

  for(i=0;i<numBasePolynoms;i++){
    (*mat)[i][0]=1.;
    (*mat)[i+numBasePolynoms][1]=1.;
    (*mat)[i][2]=Q[i];
    (*mat)[i+numBasePolynoms][3]=Q[i+numBasePolynoms];
    (*mat)[i+numBasePolynoms][4]=Q[i];
    (*mat)[i][5]=Q[i+numBasePolynoms];
  }
  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::SetSubShapeSpaceZero(
                                        ContourBSplineShapeMatrix* shapeSpace){
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  mat->newsize(shapeSpace->num_rows(),shapeSpace->num_cols());
  mat->SetZero();
  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

ContourBSplineShapeMatrix* ContourBSplineShapeMatrix::
  ComputeRegularisationMatrix(ContourBSplineData* data,
                              ContourBSplineShapeMatrix* shapeSpace,
                              ContourBSplineShapeMatrix* subShapeSpace){
  //build Hscript and HscriptSub
  BIAS::Matrix<double> Hscript(shapeSpace->num_cols(),shapeSpace->num_cols());
  BIAS::Matrix<double> HscriptSub(subShapeSpace->num_cols(),
                                  subShapeSpace->num_cols());
                            
  BIAS::Matrix<double> Uscript=data->splineMetricMatrixParametric_;
  
  Hscript=*shapeSpace;
  HscriptSub=*subShapeSpace;
  Hscript.MultLeft(Uscript);
  HscriptSub.MultLeft(Uscript);
  Hscript.MultLeft(shapeSpace->Transpose());
  HscriptSub.MultLeft(subShapeSpace->Transpose());

  //build (pseudo) inverses of Hscript and HscriptSub
  BIAS::SVD HscriptSVD(Hscript);
  BIAS::SVD HscriptSubSVD(HscriptSub);
  BIAS::Matrix<double> HscriptInverse = HscriptSVD.Invert();
  BIAS::Matrix<double> HscriptSubInverse = HscriptSubSVD.Invert();
  
  //build pseudo inverses for shapeSpaceMatrix and subShapeSpaceMatrix
  BIAS::SVD shapeSpaceSVD(*shapeSpace);
  BIAS::SVD subShapeSpaceSVD(*subShapeSpace);
  BIAS::Matrix<double> shapeSpaceInverse = shapeSpaceSVD.Invert();
  BIAS::Matrix<double> subShapeSpaceInverse = subShapeSpaceSVD.Invert();

// below uncommented version of inverse computing after book "active contours"
//  BIAS::Matrix<double> shapeSpaceMatrixPInverse(2*numBasePolynoms_,
//                                                2*numBasePolynoms_);
//  BIAS::Matrix<double> subShapeSpaceMatrixPInverse(
//                                               2*numBasePolynoms_,
//                                               2*numBasePolynoms_);
//
//  shapeSpaceMatrixPInverse.SetIdentity();
//  shapeSpaceMatrixPInverse.MultLeft(Uscript);
//  shapeSpaceMatrixPInverse.MultLeft(shapeSpaceMatrix_.Transpose());
//  shapeSpaceMatrixPInverse.MultLeft(HscriptInverse);
//      
//  subShapeSpaceMatrixPInverse.SetIdentity();
//  subShapeSpaceMatrixPInverse.MultLeft(Uscript);
//  subShapeSpaceMatrixPInverse.MultLeft(subShapeSpaceMatrix_.Transpose());
//  subShapeSpaceMatrixPInverse.MultLeft(HscriptSubInverse);
    
  BIAS::Matrix<double> Es(shapeSpace->num_cols(),
                          shapeSpace->num_cols());
  Es=*shapeSpace;
  Es.MultLeft(subShapeSpaceInverse);
  Es.MultLeft(*subShapeSpace);
  Es.MultLeft(shapeSpaceInverse);
  
  BIAS::Matrix<double> Ed(shapeSpace->num_cols(),
                          shapeSpace->num_cols());
  Ed.SetIdentity();
  Ed.SubIP(Es);

  //final regularisation matrix                                                            
  ContourBSplineShapeMatrix* mat = new ContourBSplineShapeMatrix;
  mat->newsize(shapeSpace->num_cols(),shapeSpace->num_cols());
  int i,j;//conter variables 
  for(i=0;i<mat->num_rows();i++){
    for(j=0;j<mat->num_cols();j++){
      (*mat)[i][j]=Ed[i][j];
    }
  }
  mat->MultLeft(Hscript);
  mat->MultLeft(Ed.Transpose());

  (mat->reference_)++;
  pntr_.push_back(mat);
  return mat;
}

//protected:

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