ContourDetectorBSpline.hh

/*
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
*/

#ifndef __CONTOURDETECTORBSPLINE_HH_
#define __CONTOURDETECTORBSPLINE_HH_


#include <cstring>
#include <cmath>
#include <vector>

#include <Base/Debug/Debug.hh>
#include <Base/Math/Vector.hh>
#include <Base/Math/Vector2.hh>
#include <Base/Image/Image.hh>

#include "ContourBSplineData.hh"
#include "ContourBSplineShapeMatrix.hh"
#include "ContourBSplineType.hh"
#include "ContourDetectorBase.hh"


#define D_CONTOURBSPLINE_DEBUG        1<<0
#define D_CONTOURBSPLINE_INIT         1<<1
#define D_CONTOURBSPLINE_INITFIT      1<<2
#define D_CONTOURBSPLINE_FIT          1<<3
#define D_CONTOURBSPLINE_BOUNDINGBOX  1<<4

namespace BIAS{
  /**
     @class ContourDetectorBSpline
     @brief represents a special B-Spline and holds functions to handle 
            fitting
   
      
     This class represents a special B-Spline curve for "Active Contours". 
     Based on book: "Active Contours" by Andrew Blake and Michael Isard.
     (http://www.robots.ox.ac.uk/~contours/)
    
     Class offers functions to manipulate Position, Orientation and Scale of 
     Curve, Drawing and Fitting Curve to gradient features in images.
   
   
    
     USAGE construction:
      ContourBSpline bspline; //call constructor
       bspline.Init(...); //initialise B-Spline
       
     The call of the Init(...) function looksup if the basis of the B-Spline
     is already computed. If not, a new ContourBSplineData Objects becomes
     generated, which holds data of this Object. If a B-Spline with same
     basis is available the pointer to the available ContourBSplineData gets
     copied. The purpose is to save/share memory. 
     This is hidden to users of this class.
     After the Init(..) call, you can draw and manipulate 
     the B-Spline. 
    
    
    
     For fitting you have to set (nonrecurring) Shape-Space, 
     invariant Shape-Space (SubShapeSpace) and compute the 
     Regularisation-Matrix. Look at "Active Contours" chapter 6 for details. 
     In common cases you want to set Shape-Space to identity 
     (SetShapeSpaceIdentity()) and invariant Shape-Space to planar affin
     transformations (SetSubShapeSpacePlanarAffin()).
    
     USAGE fitting:
       bspline.SetShapeSpace...(...);
       bspline.SetSubShapeSpace...(...);
      bspline.ComputeRegularisationMatrix();
       bspline.Fit(...);
    
     The functions Set*ShapeSpace...(...) and ComputeRegulatisationMatrix()
     register/generate a ContourBSplineShapeMatrix object.
     The purpose is to save/share memory. This should be 
     transparent to users of this class. 
            
     @ingroup g_feature     
     @author Marcel Lilienthal, 2005
   */
  template <class StorageType>
  class BIASFeatureDetector_EXPORT ContourDetectorBSpline : 
    ContourDetectorBase<StorageType>
  { 
  public:          

    /**
     * @brief standard constructor 
     * 
     *  Sets pointers to NULL and initialises some default values for Drawing
     *  and Fitting.
     */
    ContourDetectorBSpline(); 
      
    /**
     * @brief copy constructor
     * @param[in] toCopy object to copy
     * 
     * Uses operator= to copy referred B-Spline as new Object.
     * 
     */
    inline ContourDetectorBSpline(const ContourDetectorBSpline<StorageType> &toCopy){
      operator=(toCopy);
    }

    /** 
     * @brief destructor
     * 
     * Unregisters  ContourBSplineData object (with reference counting).
     * Unregisters  ContourBSplineShapeMatrix objects (with reference 
     * counting).
     * 
     */
    virtual ~ContourDetectorBSpline();

    /** 
     * @brief detect function, unfinished
     * 
     */
    int Detect(Image<StorageType>& image,
               std::vector<BIAS::BIASContour>& contour);

    /**
     * @brief builds a new ContourBspline object from given vector of 
     *        ContourBSpline objects
     * @param[in] bSplines vector of ContourBSpline objects to cluster
     * 
     * With this function you can cluster B-Spline curves in one B-Spline
     * curve. For example you can represent the contour of glasses with one
     * B-Spline curve constructed of two closed B-Spline curves.
     * 
     */
    static ContourDetectorBSpline 
      Cluster(const std::vector<ContourDetectorBSpline<StorageType> >& bSplines);


    /**
     * @brief copy operator
     * @param[in] toCopy object to copy
     */      
    inline ContourDetectorBSpline& operator=
      (const ContourDetectorBSpline<StorageType> &toCopy);
            
    /**
     * @brief initialisation function for simple curve without edges
     * @param[in] order order of B-Spline curve (implies degree)
     * @param[in] bType type of curve: Open or Closed
     * @param[in] cPnts control points of curve with x,y-coordinates
     * 
     * Generates/registers basis for B-Spline curve without edges 
     * transparently.
     * 
     * Example:
     * bspline.Init(order,BIAS::ContourBSplineType::Open,cPnts);
     */
    void Init(const unsigned int order, 
              const ContourBSplineType::Type bType,
              const std::vector<BIAS::Vector2<double> > &cPnts);

    /**
     * @brief initialisation function for curve with edges
     * @param[in] order order of B-Spline curve (implies degree)
     * @param[in] bType type of curve: Open or Closed
     * @param[in] mPnts multiplicity of control points
     * @param[in] cPnts control points of curve with x,y-coordinates
     * 
     * Generates/registers basis for B-Spline curve with edges 
     * transparently. 
     * 
     ** 
     * 
     * Param vector mPnts must have same length as vector of
     * control points. A multiplicity of "order-1" at index "i" implies an 
     * edge in the curve at control point with index "i". 
     *
     * If type of the curve is "Open" the first control point with an edge
     * is only possible at index "order-1" (begin counting with zero).
     * Last possible control point with an edge is "last index - (order-1)".
     * Edges in control points must have distance in indices of "order-1".
     * Example: If order==3 (quadratic B-Spline) and there is an edge in
     * control point with index 2 (mPnts[2]=2) then the next possible edge 
     * is in control point with index 4 (mPnts[4]=2).
     * 
     * If type of curve is "Closed" the indicies of control points with edges
     * must have periodic distance of "order-1".
     * Exampe: If order==3 (quadratic B-Spline) and there is an edge in
     * control point with index 0 (mPnts[0]=2) then the next possible edge
     * is in control point with index 2 (mPnts[2]=2). Periodic distance means
     * that the last possible edge is at index "last index - 1" in this case.
     * 
     */
    void Init(const unsigned int order, 
              const ContourBSplineType::Type bType,
              const std::vector<unsigned int> &mPnts,
              const std::vector<BIAS::Vector2<double> > &cPnts);

    /**
     * @brief initialisation function for simple curve without edges
     * @param[in] order order of B-Spline curve (implies degree)
     * @param[in] bType type of curve: Open or Closed
     * @param[in] Q control points of curve first all x-coordinates then all
     *              y-coordinates (adopted from book "Active Contours")
     *            Q= (x) - vector of x-coordinates
     *               (y) - vector of y-coordinates
     * 
     * Generates/registers basis for B-Spline curve without edges 
     * transparently.
     * 
     */
    void Init(const unsigned int order, 
              const ContourBSplineType::Type bType,
              const BIAS::Vector<double>& Q);

    /**
     * @brief initialisation function for curve with edges
     * @param[in] order order of B-Spline curve (implies degree)
     * @param[in] bType type of curve: Open or Closed
     * @param[in] mPnts multiplicity of control points
     * @param[in] Q control points of curve first all x-coordinates then all
     *              y-coordinates (adopted from book "Active Contours")
     *            Q= (x) - vector of x-coordinates
     *               (y) - vector of y-coordinates
     * 
     * Generates/registers basis for B-Spline curve with edges 
     * transparently. 
     * 
     ** 
     * 
     * Param vector mPnts must have same length as vector of
     * control points. A multiplicity of "order-1" at index "i" implies an 
     * edge in the curve at control point with index "i". 
     *
     * If type of the curve is "Open" the first control point with an edge
     * is only possible at index "order-1" (begin counting with zero).
     * Last possible control point with an edge is "last index - (order-1)".
     * Edges in control points must have distance in indices of "order-1".
     * Example: If order==3 (quadratic B-Spline) and there is an edge in
     * control point with index 2 (mPnts[2]=2) then the next possible edge 
     * is in control point with index 4 (mPnts[4]=2).
     * 
     * If type of curve is "Closed" the indicies of control points with edges
     * must have periodic distance of "order-1".
     * Exampe: If order==3 (quadratic B-Spline) and there is an edge in
     * control point with index 0 (mPnts[0]=2) then the next possible edge
     * is in control point with index 2 (mPnts[2]=2). Periodic distance means
     * that the last possible edge is at index "last index - 1" in this case.
     * 
     */
    void Init(const unsigned int order, 
              const ContourBSplineType::Type bType,
              const std::vector<unsigned int>& mPnts,
              const BIAS::Vector<double>& Q);
  
    /**
     * @brief fits B-Spline curve to gradient features found in the image
     * @param[in] greyImg input image
     * @param[in] alpha regularisation factor for fitting
     * @param[in] normalWidth width of normal for searching gradient features
     *
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for features
     * (normalWidth).
     * 
     * Along the curve normals get computed. Number of normals depends
     * on variable fitSampleWidth_ (SetFitSampleWidth(...)). 
     * Along the normals the feature points get assigned by the image 
     * points width highest gradient and distance less then normalWidth 
     * from curve (in both directions). Thus a feature curve is given.
     * 
     * The fitting algorithm is described in the book "Active Contours",
     * chapter 6. The parameter alpha determines the regularisation of the
     * fitting. For numeric stability alpha has to be greater than 0 even if
     * no regularisation is forced. the regularisation determines how the
     * curve should be fitted to the feature curve. high values of alpha 
     * cause the fitted curve to keep the ground shape and only adapt
     * position, orientation and scale. small values causes a possibly
     * distorted shape which could have nothing in common with the initial
     * curve.
     * 
     */
    void Fit(BIAS::Image<StorageType>& greyImg, const double& alpha, 
             const unsigned int normalWidth);

    /**
     * @brief fits B-Spline curve to gradient features found in the image
     * @param[in] greyImg input image
     * @param[in] alpha regularisation factor for fitting
     * @param[in] normalWidth width of normal for searching gradient features
     * @param[in] minClip value of image gradients which are discarded in
     *            the search after feature points
     *
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for features
     * (normalWidth).
     * 
     * Along the curve normals get computed. Number of normals depends
     * on variable fitSampleWidth_ (SetFitSampleWidth(...)). 
     * Along the normals the feature points get assigned by the image 
     * points width highest gradient and distance less then normalWidth 
     * from curve (in both directions). Image gradients with a value less 
     * than parameter minClip are discarded.
     * Thus a feature curve is given.
     * 
     * The fitting algorithm is described in the book "Active Contours",
     * chapter 6. The parameter alpha determines the regularisation of the
     * fitting. For numeric stability alpha has to be greater than 0 even if
     * no regularisation is forced. the regularisation determines how the
     * curve should be fitted to the feature curve. high values of alpha 
     * cause the fitted curve to keep the ground shape and only adapt
     * position, orientation and scale. small values causes a possibly
     * distorted shape which could have nothing in common with the initial
     * curve.
     * 
     */
    void Fit(BIAS::Image<StorageType>& greyImg, const double& alpha, 
             const unsigned int normalWidth, const double& minClip);
 
    /**
     * @brief fits B-Spline curve to gradient features found in the image
     * @param[in] featImg input gradient image
     * @param[in] alpha regularisation factor for fitting
     * @param[in] normalWidth width of normal for searching gradient features
     *
     * Along the curve normals get computed. Number of normals depends
     * on variable fitSampleWidth_ (SetFitSampleWidth(...)). 
     * Along the normals the feature points get assigned by the image 
     * points width highest gradient and distance less then normalWidth 
     * from curve (in both directions). Thus a feature curve is given.
     * 
     * The fitting algorithm is described in the book "Active Contours",
     * chapter 6. The parameter alpha determines the regularisation of the
     * fitting. For numeric stability alpha has to be greater than 0 even if
     * no regularisation is forced. the regularisation determines how the
     * curve should be fitted to the feature curve. high values of alpha 
     * cause the fitted curve to keep the ground shape and only adapt
     * position, orientation and scale. small values causes a possibly
     * distorted shape which could have nothing in common with the initial
     * curve.
     * 
     */              
    void Fit(const BIAS::Image<float>& featImg, const double& alpha, 
             const unsigned int normalWidth);

    /**
     * @brief fits B-Spline curve to gradient features found in the image
     * @param[in] featImg input gradient image
     * @param[in] alpha regularisation factor for fitting
     * @param[in] normalWidth width of normal for searching gradient features
     * @param[in] minClip value of image gradients which are discarded in
     *            the search after feature points
     *
     * Along the curve normals get computed. Number of normals depends
     * on variable fitSampleWidth_ (SetFitSampleWidth(...)). 
     * Along the normals the feature points get assigned by the image 
     * points width highest gradient and distance less then normalWidth 
     * from curve (in both directions). Image gradients with a value less 
     * than parameter minClip are discarded.
     * Thus a feature curve is given.
     * 
     * The fitting algorithm is described in the book "Active Contours",
     * chapter 6. The parameter alpha determines the regularisation of the
     * fitting. For numeric stability alpha has to be greater than 0 even if
     * no regularisation is forced. the regularisation determines how the
     * curve should be fitted to the feature curve. high values of alpha 
     * cause the fitted curve to keep the ground shape and only adapt
     * position, orientation and scale. small values causes a possibly
     * distorted shape which could have nothing in common with the initial
     * curve.
     * 
     */          
    void Fit(const BIAS::Image<float>& featImg, const double& alpha, 
             const unsigned int normalWidth, const double& minClip);

    /**
     * @brief saves all values which are needed to reinitialise the object
     *        to a binary file
     * @param[in] filename name of destination file 
     */      
    bool Save(const std::string& filename);

    /**
     * @brief load/initialise a ContourBSpline object from an given binary
     *        file
     * @param[in] filename name of source file
     */
    bool Load(const std::string& filename);

    //setter

    /**
     * @brief sets new control points
     * @param[in] controlPoints vector of control points with x,y-coordinates
     */      
    inline void SetControlPoints(
                                 const std::vector<BIAS::Vector2<double> >& controlPoints);

    /**
     * @brief sets new control points
     * @param[in] Q vector of control points consisting of vector of 
     *              x-coordinates and of vector of y-coordinates
     *              Q= (x) vector of x-coordinates
     *                 (y) vector of y-coordinates
     */
    inline void SetControlPoints(const BIAS::Vector<double>& Q);
  
    /**
     * @brief displace all control points
     * @param[in] x x-offset to x-coordinates of control points
     * @param[in] y y-offset to y-coordinates of control points
     */  
    inline void Displace(const double& x, const double&y);
      
    /**
     * @brief scales curve isotropic
     * @param[in] scale scaling factor
     * 
     * Scales the curve by factor "scale". Curve must have its origin at 
     * zero! There is no test!
     */
    inline void Scale(const double& scale);
      
    /**
     * @brief scales curve in x-direction (horizontally)
     * @param[in] scale scaling factor
     * 
     * Scales the curve by factor "scale" horizontally. Curve must have 
     * its origin at zero! There is no test!
     */
    inline void ScaleX(const double& scale);
      
    /**
     * @brief scales curve in y-direction (vertically)
     * @param[in] scale scaling factor
     * 
     * Scales the curve by factor "scale" vertically. Curve must have 
     * its origin at zero! There is no test!
     */
    inline void ScaleY(const double& scale);
    
    /**
     * @brief rotates curve by angle (rad)
     * @param[in] angle rotation angle
     * 
     * Rotates curve by angle. Curve must have its origin at zero! There is
     * no test!
     */  
    void Rotate(const double& angle);
      
    /**
     * @brief centers the curve by computing center over control points 
     *        (not accurate!)
     */
    inline void CenterOverControlPoints(); 

    /**
     * @brief sets spape-space for regularised fitting 
     * @param[in] shapeSpace matrix of shape-space 
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetShapeSpaceMatrix(BIAS::Matrix<double>& shapeSpace);
      
    /**
     * @brief sets spape-space for regularised fitting to identity
     * 
     * Read book "Active Contours" chapter 6 for details. 
     * If youre not sure which shape-space to set... set this.
     */
    inline void SetShapeSpaceIdentity();

    /**
     * @brief sets spape-space for regularised fitting to euclidian 
     *        similarities
     * 
     * Read book "Active Contours" chapter 6 for details. 
     */ 
    inline void SetShapeSpaceEuclidian();
      
    /**
     * @brief sets spape-space for regularised fitting to planar affine
     *        transformations
     * 
     * Read book "Active Contours" chapter 6 for details. 
     */ 
    inline void SetShapeSpacePlanarAffin();     

    /**
     * @brief set invariant shape-space for regularised fitting
     * @param[in] subShapeSpace matrix of invariant shape-space
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetSubShapeSpaceMatrix(BIAS::Matrix<double>& subShapeSpace);

    /**
     * @brief set invariant shape-space for regularised fitting to identity
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetSubShapeSpaceIdentity();
      
    /**
     * @brief set invariant shape-space for regularised fitting to euclidian
     *        similarities
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetSubShapeSpaceEuclidian();
      
    /**
     * @brief set invariant shape-space for regularised fitting to planar
     *        affine transformation
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetSubShapeSpacePlanarAffin();
      
    /**
     * @brief set invariant shape-space for regularised fitting to zero
     * 
     * If you set invariant shape-space matrix to zero, you have to set
     * the normal shape-space before (SetShapeSpace...(...))!
     * 
     * Read book "Active Contours" chapter 6 for details.
     */
    inline void SetSubShapeSpaceZero();

    /**
     * @brief computes the regularisation matrix needed for fitting 
     *        algorithm. 
     * 
     * The regularisation matrix for the fitting algorithm gets computed.
     * You have to set shape-space and invariant shape-space before!
     * Different to the book "Active Contours" the regularisation factor
     * alpha isnt included in the resulting matrix.
     * 
     **
     * 
     * Usage:
     * bspline.SetShapeSpace...(...);
     * bspline.SetSubShapeSpace...(...);
     * bspline.ComputeRegularisationMatrix(); 
     * 
     */
    inline void ComputeRegularisationMatrix();
      
    /**
     * @brief sets the amount of samples used to search
     *        features along the curve in fitting algorithm.
     * @param[in] width distance of sample points along the curve (<1.0)
     * 
     * Computes the amount of samples which are used to search features 
     * along the curve in the fitting algorithm. These sample points have 
     * distance "width". "width" should be considerably <1.0. 
     * the sample width is equidistant in sampling the basis of the curve,
     * but not in sampling the curve points!
     * 
     */
    inline void SetFitSampleWidth(double width);
      
    /**
     * @brief sets the amount of samples used to search features along the
     *        curve in fitting algorithm.
     * @param[in] samples amount of samples used to search features
     * 
     * Sets amount of samples which are used to search features along the
     * curve in the fitting algorithm. Distance between sample points 
     * gets computed (sample width).
     * the sample width is equidistant in sampling the basis of the curve,
     * but not in sampling the curve points!
     * 
     */
    inline void SetFitSamples(unsigned int samples);

    /**
     * @brief sets the amount of sample points used to draw the B-Spline
     *        curve.
     * @param[in] width distance of sample points along the curve (<1.0)
     * 
     * Computes the amount of samples which are used to draw the curve.
     * These sample points have distance "width". "width" should be 
     * considerably <1.0. 
     * the sample width is equidistant in sampling the basis of the curve,
     * but not in sampling the curve points!
     * 
     */      
    inline void SetDrawSampleWidth(double width);

    /**
     * @brief sets the amount of sample points used to draw the B-Spline
     *        curve.
     * @param[in] samples amount of samples used to draw curve
     * 
     * Sets amount of samples which are used to draw the curve. Distance
     * between sample points gets computed (sample width).
     * the sample width is equidistant in sampling the basis of the curve,
     * but not in sampling the curve points!
     * 
     */
    inline void SetDrawSamples(unsigned int samples);
      
    //getter

    /** 
     * @brief returns a vector of control points; first entries are 
     *        x-coordinates; last entries are y-coordinates
     * @param[out] Q vector of control points
     */
    inline void GetControlPoints(BIAS::Vector<double>& Q){Q=Q_;};
      
    /**
     * @brief returns a vector of control points in x,y-coordinates
     * @param[out] cPnts vector of control points in x,y-coordinates
     */
    inline void GetControlPoints(std::vector<BIAS::Vector2<double> >& cPnts){
      QToCPnts_(cPnts);
    };
      
    /**
     * @brief returns order of B-Spline curve
     */
    inline unsigned int GetOrder();

    /**
     * @brief returns type of B-Spline curve - possible values are Open,
     *        Closed or Cluster.
     */      
    inline ContourBSplineType::Type GetBType();       

    /**
     * @brief returns vector of multiple points; indicates which
     *        control points are modelled as edges
     */
    inline std::vector<unsigned int> GetMultiplePoints();       

    /**
     * @brief returns point on the curve
     * @param[in] t sample point in [0,1]
     * @param[out] res resulting point 
     * 
     * Parameter "t" is used to sample the basis of curve and compute
     * the curve point at "t".
     * 
     */
    void GetPoint(const double& t, BIAS::Vector2<double>& res);
      
    /**
     * @brief returns normal in curve point 
     * @param[in] t sample point in [0,1]
     * @param[out] res resulting normal
     * 
     * Parameter "t" is used to sample the basis of curve and compute the
     * normal in curve point at "t".
     * 
     */
    void GetNormal(const double& t, BIAS::Vector2<double>& res);
      
    /**
     * @brief returns the feature point with highest gradient along the 
     *        normal in curve point located at "t"
     * @param[in] greyImg input grey image
     * @param[in] t sample point in [0,1]
     * @param[in] normalWidth length of normal for searching gradient 
     *            feature along normal in curve point at "t"
     * @param[out] res resulting feature point
     * 
     * 
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for feature
     * (normalWidth).
     * 
     * Parameter "t" is used to sample the basis of the curve and computing
     * the curve point at "t". the normal at this curve point gets computed.
     * Along this normal (in both directions) the image point width the 
     * highest gradient gets determined. The resulting features point
     * has distance less equal than "normalWidth" to curve point.
     * 
     */
    void GetFeature(BIAS::Image<StorageType>& greyImg,
                    const double& t, const unsigned int normalWidth, 
                    BIAS::Vector2<double>& res);

    /**
     * @brief returns the feature point with highest gradient along the 
     *        normal in curve point located at "t"
     * @param[in] greyImg input grey image
     * @param[in] t sample point in [0,1]
     * @param[in] normalWidth length of normal for searching gradient 
     *            feature along normal in curve point at "t"
     * @param[in] minClip value of image gradients which are discarded in
     *            the search after the feature point
     * @param[out] res resulting feature point
     * 
     * 
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for feature
     * (normalWidth).
     * 
     * Parameter "t" is used to sample the basis of the curve and computing
     * the curve point at "t". the normal at this curve point gets computed.
     * Along this normal (in both directions) the image point width the 
     * highest gradient gets determined. image gradients less than "minClip"
     * get discarded. The resulting features point has distance less equal 
     * then "normalWidth" to curve point.
     * 
     */
    void GetFeature(BIAS::Image<StorageType>& greyImg,
                    const double& t, const unsigned int normalWidth,
                    const double& minClip, 
                    BIAS::Vector2<double>& res);

    /**
     * @brief returns the feature point with highest gradient along the 
     *        normal in curve point located at "t"
     * @param[in] featImg input gradient image
     * @param[in] t sample point in [0,1]
     * @param[in] normalWidth length of normal for searching gradient 
     *            feature along normal in curve point at "t"
     * @param[out] res resulting feature point
     * 
     * Parameter "t" is used to sample the basis of the curve and computing
     * the curve point at "t". the normal at this curve point gets computed.
     * Along this normal (in both directions) the image point width the 
     * highest gradient gets determined. The resulting feature point
     * has distance less equal than "normalWidth" to curve point.
     * 
     */
    void GetFeature(const BIAS::Image<float>& featImg, 
                    const double& t, const unsigned int normalWidth, 
                    BIAS::Vector2<double>& res);

    /**
     * @brief returns the feature point with highest gradient along the 
     *        normal in curve point located at "t"
     * @param[in] featImg input gradient image
     * @param[in] t sample point in [0,1]
     * @param[in] normalWidth length of normal for searching gradient 
     *            feature along normal at curve point at "t"
     * @param[in] minClip value of image gradients which are discarded in
     *            the search after the feature point
     * @param[out] res resulting feature point
     * 
     * 
     * Parameter "t" is used to sample the basis of the curve and computing
     * the curve point at "t". the normal at this curve point gets computed.
     * Along this normal (in both directions) the image point width the 
     * highest gradient gets determined. image gradients less than "minClip"
     * get discarded. The resulting feature point has distance less equal 
     * than "normalWidth" to curve point.
     * 
     */
    void GetFeature(const BIAS::Image<float>& featImg, 
                    const double& t, const unsigned int normalWidth,
                    const double& minClip, 
                    BIAS::Vector2<double>& res);

    /**
     * @brief returns the feature point with highest gradient along
     *        the normal in curve point "curvePnt"
     * @param[in] greyImg input grey image
     * @param[in] normalWidth length of normal for searching gradient
     *            feature along normal at curve point "curvePnt"
     * @param[in] curvePnt point on the curve which is start point for the
     *            search for feature point
     * @param[in] normal normal at curve point "curvePnt"
     * @param[out] res resulting feature point
     *
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for feature
     * (normalWidth).
     *  
     * Along the normal (in both directions) at curve point "curvePnt"
     * the image point width the highest gradient gets determined. The 
     * resulting feature point has distance less equal than "normalWidth"
     * to the curve point.
     * 
     */
    inline void GetFeature(BIAS::Image<StorageType>& greyImg, 
                           const unsigned int normalWidth,
                           const BIAS::Vector2<double>& curvePnt,
                           const BIAS::Vector2<double>& normal, 
                           BIAS::Vector2<double> & res);
                             
    /**
     * @brief returns the feature point with highest gradient along
     *        the normal in curve point "curvePnt"
     * @param[in] greyImg input grey image
     * @param[in] normalWidth length of normal for searching gradient
     *            feature along normal at curve point "curvePnt"
     * @param[in] curvePnt point on the curve which is start point for the
     *            search for feature point
     * @param[in] normal normal at curve point "curvePnt"
     * @param[in] minClip value of image gradients which are discarded in
     *            the search for the feature point
     * @param[out] res resulting feature point
     *
     * The input image gets converted to a single gradient image. A ROI gets
     * computed from controlpoints and the width to search for feature
     * (normalWidth).
     *  
     * Along the normal (in both directions) at curve point "curvePnt"
     * the image point width the highest gradient gets determined. 
     * Image gradients less than "minClip"get discarded.The 
     * resulting feature point has distance less equal than "normalWidth"
     * to the curve point.
     * 
     */
    inline void GetFeature(BIAS::Image<StorageType>& greyImg, 
                           const unsigned int normalWidth,
                           const BIAS::Vector2<double>& curvePnt,
                           const BIAS::Vector2<double>& normal,
                           const double& minClip,
                           BIAS::Vector2<double> & res);

    /**
     * @brief returns the feature point with highest gradient along
     *        the normal in curve point "curvePnt"
     * @param[in] featImg input gradient image
     * @param[in] normalWidth length of normal for searching gradient
     *            feature along normal at curve point "curvePnt"
     * @param[in] curvePnt point on the curve which is start point for the
     *            search for feature point
     * @param[in] normal normal at curve point "curvePnt"
     * @param[out] res resulting feature point
     *
     * Along the normal (in both directions) at curve point "curvePnt"
     * the image point width the highest gradient gets determined. The 
     * resulting feature point has distance less equal than "normalWidth"
     * to the curve point.
     * 
     */
    void GetFeature(const BIAS::Image<float>& featImg, 
                    const unsigned int normalWidth,
                    const BIAS::Vector2<double>& curvePnt,
                    const BIAS::Vector2<double>& normal, 
                    BIAS::Vector2<double> & res);

    /**
     * @brief returns the feature point with highest gradient along
     *        the normal in curve point "curvePnt"
     * @param[in] greyImg input gradient image
     * @param[in] normalWidth length of normal for searching gradient
     *            feature along normal at curve point "curvePnt"
     * @param[in] curvePnt point on the curve which is start point for the
     *            search for feature point
     * @param[in] normal normal at curve point "curvePnt"
     * @param[in] minClip value of image gradients which are discarded in
     *            the search for the feature point
     * @param[out] res resulting feature point
     *
     * Along the normal (in both directions) at curve point "curvePnt"
     * the image point width the highest gradient gets determined. 
     * Image gradients less than "minClip"get discarded.The 
     * resulting feature point has distance less equal than "normalWidth"
     * to the curve point.
     * 
     */
    void GetFeature(const BIAS::Image<float>& featImg, 
                    const unsigned int normalWidth,
                    const BIAS::Vector2<double>& curvePnt,
                    const BIAS::Vector2<double>& normal,
                    const double& minClip, 
                    BIAS::Vector2<double> & res);

    /**
     * @brief returns the parameter of the bounding box of the curve
     * @param[out] minX x-coordinate of top left corner of bounding box
     * @param[out] minY y-coordinate of top left corner of bounding box
     * @param[out] maxX x-coordinate of bottom right corner of bounding box
     * @param[out] maxX y-coordinate of bottom right corner of bounding box
     */
    void GetBoundingBox(int& minX, int& minY, int& maxX, int& maxY);
      
    //drawing
      
    /**
     * @brief draws the control points of the B-Spline curve into image
     * @param[out] img output image
     * @param[in] color
     */
    void DrawControlPoints(BIAS::Image<StorageType>& img,
                           const StorageType color[3]);

    /**
     * @brief draws B-Spline curve into image
     * @param[out] img output image
     * @param[in] color
     * 
     * Draws B-Spline curve into given image. Uses amount of samples given
     * by SetDrawSampleWidth(...) of SetDrawSamples(...). If amount of
     * samples isnt given by these function a default amount is used.
     * 
     */
    void DrawCurve(BIAS::Image<StorageType> &img,
                   const StorageType color[3]);

    /**
     * @brief draws B-Spline curve into image
     * @param[out] img output image
     * 
     * Draws B-Spline curve into given image. Uses amount of samples given
     * by SetDrawSampleWidth(...) or SetDrawSamples(...). If amount of
     * samples isnt given by these function a default amount is used.
     * 
     */
    void DrawCurve(BIAS::Image<float> &img);
      
    /**
     * @brief draws B-Spline normals into image
     * @param[out] img output image
     * 
     * Draws normals along the B-Spline curve into given image. Uses
     * amount of samples given by SetFitSampleWidth(...) or 
     * SetFitSamples(...). If amount of samples isnt given a default amount
     * is used.
     * 
     */
    void DrawNormals(const unsigned int normalWidth,
                     BIAS::Image<StorageType> &img,
                     const StorageType color[3]);


  protected:
    //methods

    /**
     * @brief initializes curve for drawing or for getting points on it
     * 
     * Computes curve vectors needed for determine the shape of the curve.
     * Depends on the basis and the current controlpoints.
     */
    void InitCurve_();
      
    /**
     * @brief converts a given control point vector in x,y-coordinates to
     *        the internal used flat control point vector Q whose first 
     *        entries are x-coordinates and last entries are y-coordinates
     * @param[in] cPnts vector of control points in x,y-coordinates
     */
    void CPntsToQ_(const std::vector<BIAS::Vector2<double> >& cPnts);

    /**
     * @brief converts a given control point vector in x,y-coordinates to
     *        the internal used flat control point vector Q whose first 
     *        entries are x-coordinates and last entries are y-coordinates
     * @param[in] cPnts vector of control points in x,y-coordinates
     * @param[out] Q vector of control points; first entries are 
     *               x-coordinates; last entries are y-coordinates
     */
    void CPntsToQ_(const std::vector<BIAS::Vector2<double> >& cPnts,
                   BIAS::Vector<double>& Q);
      
    /**
     * @brief converts a flat control point vector Q whose first entries 
     *        are x-coordinates and last entries are y-coordinates to a
     *        control point vector with x,y-coordinates
     * @param[out] cPnts control point vector with x,y-coordinates
     */
    void QToCPnts_(std::vector<BIAS::Vector2<double> >& cPnts);

    /**
     * @brief converts a flat control point vector Q whose first entries 
     *        are x-coordinates and last entries are y-coordinates to a
     *        control point vector with x,y-coordinates
     * @param[in] Q vector of control points; first entries are
     *              x-coordinates; last entries are y-coordinates
     * @param[out] cPnts control point vector with x,y-coordinates
     */
    void QToCPnts_(const BIAS::Vector<double>& Q, 
                   std::vector<BIAS::Vector2<double> >& cPnts);
                    
    /**
     * @brief computes a point on curve and a parameter vector (byproduct)
     *        which can be used for computing a special matrix (U) in 
     *        the fitting algorithm
     * @param[in] numSpan span of curve in which point is located
     * @param[in] s sample point of span in [0,1]
     * @param[out] resParamVec special parameter vector (byproduct) used in
     *             fitting algorithm
     * @param[out] res resulting point on curve
     * 
     * This function is internal used to compute a point on a curve. Each
     * B-Spline-Curve consist of a number of spans. Parameter "numSpan" and
     * "s" are used to sample the curve basis.
     * Parameter "resParamVec" is a byproduct of curve point computation. It
     * is build up on basis of "s" and can be used to compute the special 
     * matrix "U" in fitting algorihm.
     * More information about "U" can be found in the book "Active Contours".
     * 
     */
    void GetPointAndParamVec_(const unsigned int numSpan, const double& s, 
                              BIAS::Vector<double>& resParamVec,
                              BIAS::Vector2<double>& res);

    /**
     * @brief  computes a point on curve
     * @param[in] numSpan span of curve in which point is located
     * @param[in] s sample point of span in [0,1]
     * @param[out] res resulting point on curve
     * 
     * This function is internal used to compute a point on a curve. Each
     * B-Spline-Curve consist of a number of spans. Parameter "numSpan" and
     * "s" are used to sample the curve basis.
     */
    inline void GetPoint_(const unsigned int numSpan, const double& s,
                          BIAS::Vector2<double>& res);

    /**
     * @brief computes a normal at curve point which is parametrized by 
     *        "numSpan" and "s"
     * @param[in] numSpan span of curve in which normal is located
     * @param[in] s sample point of span in [0,1]
     * @param[out] res resulting normal
     * 
     * This function is internal used to compute a normal at a curve point.
     * Each B-Spline-Curve consists of a number of spans.
     * The curve point is parameterized by "numSpan" and "s".
     */
    void GetNormal_(const unsigned int numSpan, const double& s, 
                    BIAS::Vector2<double>& res);

    /**
     * @brief converts given grey image to gradient image with ROI
     * @param[in] greyImg input grey image
     * @param[out] featImg output gradient image
     * @param[in] normalWidth used to determine ROI
     * 
     * Converts a given grey image to a gradient image. gradient image is
     * needed for the fitting algorithm. In the fitting algorithm only a
     * small region of the image is of interest for searching feature points.
     * Thus a ROI gets computed from a bounding box of the curve and the
     * parameter "normalWidth". Parameter "normalWidth" determines how 
     * far along curve normal a feature point gets searched.
     * 
     */
    void GetFeatImage_(BIAS::Image<StorageType>& greyImg, 
                       BIAS::Image<float>& featImg, 
                       const unsigned int normalWidth);

    /**
     * @brief fits the curve to a features curve with regularized fitting
     * @param[in] featImage input gradient image
     * @param[in] alpha regularisation factor
     * @param[in] normalWidth length of normal used to search for feature
     *            curve
     * 
     * Internal used function to fit the B-Spline curve to a feature curve.
     * The feature curve is determined by searching for image points with
     * highest gradients along the normals at the curve. Parameter
     * "normalWidth" determines the maximal distance of feature point to
     * curve (in both directions).
     * 
     * More details in function Fit(..) and book "Active Contours" chapter 6.
     * 
     */
    void Fit_(const BIAS::Image<float>& featImg, const double& alpha, 
              const unsigned int normalWidth);

    /**
     * @brief fits the curve to a features curve with regularized fitting
     * @param[in] featImage input gradient image
     * @param[in] alpha regularisation factor
     * @param[in] normalWidth length of normal used to search for feature
     *            curve
     * @param[in] minClip gradients less than "minClip" are discarded in 
     *                    search for feature points
     * 
     * Internal used function to fit the B-Spline curve to a feature curve.
     * The feature curve is determined by searching for image points with
     * highest gradients along the normals at the curve. Parameter
     * "normalWidth" determines the maximal distance of feature point to
     * curve (in both directions).
     * 
     * More details in function Fit(..) and book "Active Contours" chapter 6.
     * 
     */
    void Fit_(const BIAS::Image<float>& featImg, const double& alpha, 
              const unsigned int normalWidth, const double& minClip);

    //DEBUG methods
      
#ifdef BIAS_DEBUG
    /**
     * @brief tests if B-Spline curve with given parameter can get 
     *        initialized.
     * @param[in] order order of B-Spline curve
     * @param[in] bType type of B-Spline curve (Open or Closed)
     * @param[in] mPnts vector of multiplicities - determines which control
     *                  points are modelled as edges
     * @param[in] Q vector of control points
     * 
     * Test if B-Spline curve with given parameter can get initialised.
     * Most important test if mPnts is correct. That means if edges could
     * be generated in these control points.
     */
    void TestInit_(const unsigned int order, 
                   const ContourBSplineType::Type bType,
                   const std::vector<unsigned int>& mPnts,
                   const BIAS::Vector<double>& Q);                    
#endif



    //class variables
      
    /**
     * vector of control points; first entries are x-coordinates; last 
     * entries are y-coordinates
     */
    BIAS::Vector<double> Q_; //control points 
      
    //! pointer to data of B-Spline curve initialised by call of Init(...)
    ContourBSplineData* pData_;
      
    //! vectors needed to draw curve (determined by InitCurve_)
    std::vector<std::vector<BIAS::Vector<double> > > curveVectors_;
    //! bool to indicate a change of control points 
    bool curveInitialised_;

    //! pointer to shape-space matrix object
    ContourBSplineShapeMatrix* pShapeSpaceMatrix_;
    //! pointer to invariant shape-space matrix object
    ContourBSplineShapeMatrix* pSubShapeSpaceMatrix_;
    //! pointer to regularisation matrix object
    ContourBSplineShapeMatrix* pRegMatrix_;//regularisation Matrix
      
    //! distance of sample points in basis used in fitting algorithm
    double* fitSampleWidth_;
    //! distance of sample points in basis used in drawing algorithm      
    double* drawSampleWidth_;
    //! amount of samples used in fitting algorithm
    unsigned int fitSamples_;
    //! amount of samples used in drawing algorithm
    unsigned int drawSamples_;
      
  };//class ContourBSpline

   template <class StorageType> ContourDetectorBSpline<StorageType>& 
   ContourDetectorBSpline<StorageType>::
   operator=(const ContourDetectorBSpline<StorageType>& toCopy){
    //copies all important class member variables from object toCopy
      Q_=toCopy.Q_;
      pData_=toCopy.pData_;
      if(pData_!=NULL){
        (pData_->reference_)++; //increase data reference counter
      }        
      curveInitialised_=false; //curveVectors_ get not copied
      pShapeSpaceMatrix_=toCopy.pShapeSpaceMatrix_;
      if(pShapeSpaceMatrix_!=NULL){
        (pShapeSpaceMatrix_->reference_)++;//increase ref counter
      }
      pSubShapeSpaceMatrix_=toCopy.pSubShapeSpaceMatrix_;
      if(pSubShapeSpaceMatrix_!=NULL){
        (pSubShapeSpaceMatrix_->reference_)++;//increase ref counter
      }
      pRegMatrix_=toCopy.pRegMatrix_;
      if(pRegMatrix_!=NULL){
        (pRegMatrix_->reference_)++;//increase ref counter
      }
      fitSampleWidth_=toCopy.fitSampleWidth_;
      drawSampleWidth_=toCopy.drawSampleWidth_;
      fitSamples_=toCopy.fitSamples_;
      drawSamples_=toCopy.drawSamples_;
  
      return *this;
    }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetControlPoints(const std::vector<BIAS::Vector2<double> >& controlPoints){
    //convert vector of control points into internal used format
    CPntsToQ_(controlPoints);
    curveInitialised_=false;
  }

  
  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetControlPoints(const BIAS::Vector<double>& Q){
    Q_=Q;
    curveInitialised_=false;
  }

  
  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  Displace(const double& x, const double&y){
    unsigned int i; //counter variable
    for(i=0;i<(pData_->numBasePolynoms_);i++){
      Q_[i]+=x;
      Q_[i+(pData_->numBasePolynoms_)]+=y;
    }
    curveInitialised_=false;
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  Scale(const double& scale){
    //curve must have origin at zero!
    int i; // counter variable;
    for(i=0;i<Q_.size();i++){
      Q_[i]*=scale;
    }
  }

   template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  ScaleX(const double& scale){
    //curve must have origin at zero!
    int i; // counter variable;
    for(i=0;i<(Q_.size()/2);i++){
      Q_[i]*=scale;
    }
  }

   template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  ScaleY(const double& scale){
    //curve must have origin at zero!
    int i; // counter variable;
    for(i=(Q_.size()/2);i<Q_.size();i++){
      Q_[i]*=scale;
    }
  }

  template <class StorageType> void 
  ContourDetectorBSpline<StorageType>::
  CenterOverControlPoints(){
    unsigned int i; //counter variable
    double x,y;
    x=0;y=0;
    for(i=0; i<pData_->numBasePolynoms_; i++){
      x+=Q_[i];
      y+=Q_[i+pData_->numBasePolynoms_];
    }
    x/=(double) pData_->numBasePolynoms_;
    y/=(double) pData_->numBasePolynoms_;
    Displace(-x,-y);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  SetShapeSpaceMatrix(
                      BIAS::Matrix<double>& shapeSpace){
    pShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpaceMatrix(
                                                                        shapeSpace);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "shapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void 
  ContourDetectorBSpline<StorageType>:: 
  SetShapeSpaceIdentity(){
    pShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpaceIdentity(
                                                                          pData_->numBasePolynoms_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "shapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void 
  ContourDetectorBSpline<StorageType>::
  SetShapeSpaceEuclidian(){
    pShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpaceEuclidian(
                                                                           pData_->numBasePolynoms_,Q_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "shapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void 
  ContourDetectorBSpline<StorageType>:: 
  SetShapeSpacePlanarAffin(){
    pShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpacePlanarAffin
      (pData_->numBasePolynoms_,Q_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "shapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  SetSubShapeSpaceMatrix(
                         BIAS::Matrix<double>& subShapeSpace){
    pSubShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpaceMatrix(
                                                                           subShapeSpace);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "subShapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pSubShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetSubShapeSpaceIdentity(){
    pSubShapeSpaceMatrix_ = ContourBSplineShapeMatrix::SetShapeSpaceIdentity
      (pData_->numBasePolynoms_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "subShapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pSubShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetSubShapeSpaceEuclidian(){
    pSubShapeSpaceMatrix_ = ContourBSplineShapeMatrix::
      SetShapeSpaceEuclidian(pData_->numBasePolynoms_,Q_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "subShapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pSubShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetSubShapeSpacePlanarAffin(){
    pSubShapeSpaceMatrix_ = ContourBSplineShapeMatrix::
      SetShapeSpacePlanarAffin(pData_->numBasePolynoms_,Q_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "subShapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pSubShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetSubShapeSpaceZero(){
    pSubShapeSpaceMatrix_ = 
      ContourBSplineShapeMatrix::SetSubShapeSpaceZero(pShapeSpaceMatrix_);
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "subShapeSpaceMatrix: " << 
              (BIAS::Matrix<double>)*pSubShapeSpaceMatrix_ <<std::endl<<
              std::flush);
  }

 template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  ComputeRegularisationMatrix(){
    if(pData_==NULL)
      BIASERR("ComputeRegularisationMatrix: pData_ NULL pointer\r\n");
    if(pShapeSpaceMatrix_==NULL)
      BIASERR("ComputeRegularisationMatrix: "<<
              "pShapeSpaceMatrix_ NULL pointer\r\n");
    if(pSubShapeSpaceMatrix_==NULL)
      BIASERR("ComputeRegularisationMatrix: "<<
              "pSubShapeSpaceMatrix_ NULL pointer\r\n");
    pRegMatrix_= ContourBSplineShapeMatrix::ComputeRegularisationMatrix(
                           pData_, pShapeSpaceMatrix_, pSubShapeSpaceMatrix_ );
    BIASCDOUT(D_CONTOURBSPLINE_INITFIT, "regMatrix: " << 
              (BIAS::Matrix<double>)*pRegMatrix_ <<std::endl<<std::flush);
  }

 template <class StorageType> void  
  ContourDetectorBSpline<StorageType>:: 
  SetFitSampleWidth(double width){
    if(fitSampleWidth_==NULL){
      fitSampleWidth_= new double;
    }
    fitSamples_=(unsigned int)(1./width);
    *fitSampleWidth_=1./(double)fitSamples_;
  }

  template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  SetFitSamples(unsigned int samples){
    if(fitSampleWidth_==NULL){
      fitSampleWidth_= new double;
    }
    fitSamples_=samples;
    *fitSampleWidth_=1./(double)fitSamples_;
  }

  template <class StorageType> void  
  ContourDetectorBSpline<StorageType>::
  SetDrawSampleWidth(double width){
    if(drawSampleWidth_==NULL){
      drawSampleWidth_= new double;
    }
    drawSamples_=(unsigned int)(1./width);
    *drawSampleWidth_=1./(double)drawSamples_;
  }

 template <class StorageType> void 
  ContourDetectorBSpline<StorageType>::
  SetDrawSamples(unsigned int samples){
    if(drawSampleWidth_==NULL){
      drawSampleWidth_= new double;
    }
    drawSamples_=samples;
    *drawSampleWidth_=1./(double)drawSamples_;
  }

  template <class StorageType> unsigned int
  ContourDetectorBSpline<StorageType>::GetOrder(){
    BIASASSERT(pData_!=NULL);
    return pData_->order_; 
  }

   template <class StorageType> ContourBSplineType::Type 
  ContourDetectorBSpline<StorageType>::GetBType(){
    BIASASSERT(pData_!=NULL);
    return pData_->bType_;
  }

   template <class StorageType> std::vector<unsigned int> 
  ContourDetectorBSpline<StorageType>::GetMultiplePoints(){
    BIASASSERT(pData_!=NULL);
    return pData_->mPnts_;
  }

   template <class StorageType> void  
  ContourDetectorBSpline<StorageType>::
  GetFeature(BIAS::Image<StorageType>& greyImg, 
             const unsigned int normalWidth,
             const BIAS::Vector2<double>& curvePnt,
             const BIAS::Vector2<double>& normal, 
             BIAS::Vector2<double> & res){
    BIAS::Image<float> featImg;
    GetFeatImage_(greyImg, featImg, normalWidth);
    GetFeature(featImg,normalWidth,curvePnt,normal,res);
  }

   template <class StorageType> void
  ContourDetectorBSpline<StorageType>:: 
  GetFeature(BIAS::Image<StorageType>& greyImg, 
             const unsigned int normalWidth,
             const BIAS::Vector2<double>& curvePnt,
             const BIAS::Vector2<double>& normal,
             const double& minClip,
             BIAS::Vector2<double> & res){
    BIAS::Image<float> featImg;
    GetFeatImage_(greyImg, featImg, normalWidth);
    GetFeature(featImg,normalWidth,curvePnt,normal,minClip,res);
  }

   template <class StorageType> void
  ContourDetectorBSpline<StorageType>::
  GetPoint_(const unsigned int numSpan, 
            const double& s,
            BIAS::Vector2<double>& res){
    BIAS::Vector<double> param;
    GetPointAndParamVec_(numSpan, s, param, res);
  }

}//namespace BIAS
#endif//__CONTOURDETECTORBSPLINE_HH_