ADD_EXECUTABLE(clitkAffineTransform clitkAffineTransform.cxx clitkAffineTransform_ggo.c)
TARGET_LINK_LIBRARIES(clitkAffineTransform clitkCommon ITKIO clitkFilters)
-ADD_EXECUTABLE(clitkSignalMeanPositionTracking clitkSignalMeanPositionTracking.cxx clitkSignalMeanPositionFilter.cxx clitkEllipse.cxx clitkSignalMeanPositionTracking_ggo.c)
-TARGET_LINK_LIBRARIES(clitkSignalMeanPositionTracking clitkCommon ITKIO )
-
-ADD_EXECUTABLE(clitkSignalApparentMotionTracking clitkSignalApparentMotionTracking.cxx clitkTrajectory2D.cxx clitkSignalApparentMotionTrackingFilter.cxx clitkEllipse.cxx clitkSignalApparentMotionTracking_ggo.c)
-TARGET_LINK_LIBRARIES(clitkSignalApparentMotionTracking clitkCommon ITKIO)
-
-ADD_EXECUTABLE(clitkSignalFilter clitkSignalFilter.cxx clitkSignalFilter_ggo.c)
-TARGET_LINK_LIBRARIES(clitkSignalFilter clitkCommon ITKIO )
-
ADD_EXECUTABLE(clitkSetBackground clitkSetBackground.cxx
clitkSetBackgroundGenericFilter.cxx clitkSetBackground_ggo.c)
TARGET_LINK_LIBRARIES(clitkSetBackground clitkCommon ITKIO )
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkEllipse.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 10:47:48 $
- Version: $Revision: 1.3 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
- =========================================================================*/
-
-#include "clitkEllipse.h"
-
-typedef itk::Vector<double,2> Vector2d;
-typedef itk::Vector<double,6> Vector6d;
-typedef itk::Matrix<double,6,6> Matrix6x6d;
-typedef itk::Matrix<double,3,3> Matrix3x3d;
-
-//---------------------------------------------------------------------
-clitk::Ellipse::Ellipse():a((*this)[0]), b((*this)[1]),
- c((*this)[2]), d((*this)[3]),
- e((*this)[4]), f((*this)[5]) {
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-clitk::Ellipse::Ellipse(const Ellipse & e):a((*this)[0]), b((*this)[1]),
- c((*this)[2]), d((*this)[3]),
- e((*this)[4]), f((*this)[5]) {
- for(int i=0; i<6; i++) (*this)[i] = e[i];
-}
-
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-Vector2d clitk::Ellipse::ComputeCenter() {
- // See http://mathworld.wolfram.com/Ellipse.html
- // see Ruan 2008
- Vector2d center;
- center[0] = (2*c*d - b*e)/(b*b-4*a*c);
- center[1] = (2*a*e - b*d)/(b*b-4*a*c);
- return center;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-Vector2d clitk::Ellipse::ComputeSemiAxeLengths() {
- // See http://www.geometrictools.com/Documentation/InformationAboutEllipses.pdf
- Vector2d axis;
- Vector2d center = ComputeCenter();
- double & k1 = center[0];
- double & k2 = center[1];
- // DD(f);
- double mu = 1.0/(a*k1*k1 + b*k1*k2 + c*k2*k2 - f);
- // DD(a*k1*k1);
- // DD(b*k1*k2);
- // DD(c*k2*k2);
- // DD(a*k1*k1 + b*k1*k2 + c*k2*k2 - f);
- // DD(mu);
- double m11 = mu * a;
- double m12 = mu * 0.5 * b;
- double m22 = mu * c;
- double l1 = ( (m11+m22) + sqrt((m11-m22)*(m11-m22)+4*m12*m12) )/2.0;
- // DD(k1);
- // DD(k2);
- // DD(mu);
- // DD(m11);
- // DD(m12);
- // DD(m22);
- // DD(l1);
- assert(l1>=0.0);
- axis[1] = 1.0/sqrt(l1);
- double l2 = ((m11+m22)-sqrt((m11-m22)*(m11-m22)+4*m12*m12))/2.0;
- // DD(l2);
- assert(l2>0.0);
- axis[0] = 1.0/sqrt(l2);
- return axis;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-double clitk::Ellipse::ComputeAngleInRad() {
- // See http://www.geometrictools.com/Documentation/InformationAboutEllipses.pdf
- double theta;
- if (b==0) {
- if (a<c) { theta = 0; }
- else { theta = 0.5*M_PI; }
- }
- else {
- if (a<c) {
- theta = 0.5*atan(-b/(c-a));
- }
- else {
- theta = M_PI/2+0.5*atan(-b/(c-a));
- }
- }
- return theta;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Ellipse::InitialiseEllipseFitting(double eta, unsigned int n, clitk::Signal & inputX, clitk::Signal & inputY) {
- // Store data
- mEta = eta;
- mInputX = &inputX;
- mInputY = &inputY;
-
- // Initialise ellipse to global fit
- std::vector<double> extremaX(2);
- std::vector<double> extremaY(2);
- double x0 = 0.0;
- double y0 = 0.0;
- extremaX[0] = extremaY[0] = numeric_limits<double>::max();
- extremaX[1] = extremaY[1] = -numeric_limits<double>::max();
- for(unsigned int i=0; i<n; i++) {
- if (inputX[i] < extremaX[0]) extremaX[0] = inputX[i];
- if (inputX[i] > extremaX[1]) extremaX[1] = inputX[i];
- if (inputY[i] < extremaY[0]) extremaY[0] = inputY[i];
- if (inputY[i] > extremaY[1]) extremaY[1] = inputY[i];
- x0 += inputX[i];
- y0 += inputY[i];
- }
- x0 /= n;
- y0 /= n;
-
- // initialisation with an ellipse more small than real points extends
- double ax1 = (extremaX[1]-extremaX[0])/2.0;
- double ax2 = (extremaY[1]-extremaY[0])/2.0;
- if (ax2 >= ax1) ax2 = 0.99*ax1;
- SetCenterAndSemiAxes(x0, y0, ax1, ax2);
-
- // Initialisation of C
- C.Fill(0.0);
- C(0,0) = 0; C(0,1) = 0; C(0,2) = 2;
- C(1,0) = 0; C(1,1) = -1; C(1,2) = 0;
- C(2,0) = 2; C(2,1) = 0; C(2,2) = 0;
- Ct = C.GetVnlMatrix().transpose();
-
- // Compute initial S
- UpdateSMatrix(0,n, inputX, inputY);
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Ellipse::CopyBlock(Matrix6x6d & out, const Matrix6x6d & in,
- int ox, int oy, int l, double factor) {
- for(int x=ox; x<ox+l; x++)
- for(int y=oy; y<oy+l; y++) {
- out(x,y) = factor * in(x,y);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Ellipse::CopyBlock(Matrix6x6d & out, const Matrix3x3d & in,
- int ox, int oy, double factor) {
- for(int x=0; x<3; x++)
- for(int y=0; y<3; y++) {
- out(ox+x,oy+y) = factor * in(x,y);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-Matrix3x3d clitk::Ellipse::GetBlock3x3(const Matrix6x6d & M, int x, int y) {
- Matrix3x3d B;
- for(int i=0; i<3; i++) {
- for(int j=0; j<3; j++) {
- B(i,j) = M(x+i,y+j);
- }
- }
- return B;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-double clitk::Ellipse::EllipseFittingNextIteration() {
- Vector6d & current = (*this);
-
- // Normalize current vector. This allow to have a decreasing
- // residual r (no need to optimize)
- GetVnlVector().normalize();
- double r = (St*current)*current;
- // DD(r);
-
- // Temporary parameters
- Vector6d an;
-
- // Iterative update
- an = Sinv * C * current;
- double num = (Wt*current)*current; // anT W an = (WT an) an
- double denom = (Ct*current)*current;
- an = (mEta * num/denom) * an;
- an += (1.0-mEta)*current;
- SetVnlVector(an.GetVnlVector());
-
- // return residual
- return r;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Ellipse::SetCenterAndSemiAxes(double x0, double y0,
- double r1, double r2) {
- a = 1.0/(r1*r1);
- b = 0;
- c = 1.0/(r2*r2);
- if (a>c) {
- std::cerr << "Error major axis should be r1, not r2 (r1=" << r1
- << " r2=" << r2 << ")" << std::endl;
- exit(0);
- }
- d = (-2.0*x0)/(r1*r1);
- e = (-2.0*y0)/(r2*r2);
- f = (x0*x0)/(r1*r1) + (y0*y0)/(r2*r2) - 1.0;
- GetVnlVector().normalize();
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Ellipse::UpdateSMatrix(unsigned int begin, unsigned int n,
- clitk::Signal & inputX, clitk::Signal & inputY) {
- // Initialisation of z
- z.resize(n);
- int j = 0;
- for(unsigned int i=begin; i<begin+n; i++) {
- z[j][0] = inputX[i]*inputX[i];
- z[j][1] = inputX[i]*inputY[i];
- z[j][2] = inputY[i]*inputY[i];
- z[j][3] = inputX[i];
- z[j][4] = inputY[i];
- z[j][5] = 1;
- j++;
- }
-
- // Initialisation of S
- S.Fill(0.0);
- // j = 0;
- for(unsigned int i=0; i<n; i++) {
- for(unsigned int x=0; x<6; x++)
- for(unsigned int y=0; y<6; y++)
- S(x,y) += z[i][x]*z[i][y];
- // j++;
- }
- Sinv = S.GetInverse();
- St = S.GetVnlMatrix().transpose();
-
- // Initialisation of W
- W.Fill(0.0);
- CopyBlock(W, S, 0, 0, 3);
- CopyBlock(W, S, 3, 3, 3, -1);
- Wt = W.GetVnlMatrix().transpose();
-
- // Automated computation of mEta
- if (mEta<0) {
- Matrix3x3d E = GetBlock3x3(S, 0, 0);
- Matrix3x3d B = GetBlock3x3(S, 0, 3);
- Matrix3x3d Bt = GetBlock3x3(S, 3, 0);
- Matrix3x3d D = GetBlock3x3(S, 3, 3);
- Matrix3x3d Dinv(D.GetInverse());
- Matrix3x3d Stilde;
- Stilde = E - B*Dinv*Bt;
-
- Matrix3x3d Ctilde;
- Ctilde(0,0) = 0; Ctilde(0,1) = 0; Ctilde(0,2) = 2;
- Ctilde(1,0) = 0; Ctilde(1,1) = -1; Ctilde(1,2) = 0;
- Ctilde(2,0) = 2; Ctilde(2,1) = 0; Ctilde(2,2) = 0;
-
- // Ctilde is not inonneg-definite, so vnl print error on std cerr
- // the following "disableStdCerr" disable it ...
- disableStdCerr();
- vnl_generalized_eigensystem solver(Stilde.GetVnlMatrix(), Ctilde.GetVnlMatrix());
- //vnl_generalized_eigensystem solver(Ctilde.GetVnlMatrix(), Stilde.GetVnlMatrix());
- enableStdCerr();
- double dmax=0.0, dmin=9999999.9;
- dmax = solver.D(2);
- dmin = solver.D(1);
- mEta = 2.0/(fabs(dmax/dmin)+1);
- assert(mEta<1.0);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::ComputeIsoPhaseFromListOfEllipses(std::vector<clitk::Ellipse*> & l,
- clitk::Signal & sx,
- clitk::Signal & sy,
- int nbIsoPhase,
- int delay,
- int L,
- clitk::Signal & phase) {
- // Init
- DD(nbIsoPhase);
- std::vector<double> mIsoPhaseRefAngle(nbIsoPhase);
- phase.resize(l.size());
- double refphaseangle=0;
- double previousangle=0;
-
- // Loop on ellipses
- // DD(l.size());
- for (unsigned int i=0; i<l.size(); i++) {
- // DD("=================================");
- // DD(i);
- clitk::Ellipse & An = *l[i];
- // DD(An);
- // DD(An.ComputeCenter());
- // DD(An.ComputeSemiAxeLengths());
- // DD(rad2deg(An.ComputeAngleInRad()));
-
- // Compute current signed angle
- Vector2d x1(An.ComputeCenter());
- double a = l[0]->ComputeSemiAxeLengths()[0];
- double theta = l[0]->ComputeAngleInRad();
- Vector2d x2; x2[0] = x1[0]+a * cos(theta); x2[1] = x1[1]+a * sin(theta);
- Vector2d x3(x1);
- Vector2d x4; x4[0] = sx[i+L]; x4[1] = sy[i+L];
- Vector2d A(x2-x1);
- Vector2d B(x4-x3);
- // DD(sx[i+L]);
- // DD(sy[i+L]);
- // http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm
- double signed_angle = atan2(B[1], B[0]) - atan2(A[1], A[0]);
- if (signed_angle<0) signed_angle = 2*M_PI+signed_angle;
-
- // First time : set the angle
- if (i==0) {
- // DD(signed_angle);
- refphaseangle = signed_angle;
- for(int a=0; a<nbIsoPhase; a++) {
- if (a==0) mIsoPhaseRefAngle[a] = refphaseangle;//signed_angle;
- else mIsoPhaseRefAngle[a] = (refphaseangle //signed_angle
- + (2*M_PI)*a/nbIsoPhase);
- if (mIsoPhaseRefAngle[a] > 2*M_PI) mIsoPhaseRefAngle[a] -= 2*M_PI;
- if (mIsoPhaseRefAngle[a] < 0) mIsoPhaseRefAngle[a] = 2*M_PI-mIsoPhaseRefAngle[a];
- }
- int a=0;
- while ((a<nbIsoPhase) && (signed_angle>=mIsoPhaseRefAngle[a])) { a++; }
- phase[i] = a-1;
- if (nbIsoPhase == 1) phase[i] = 1;
- }
- else {
- phase[i] = phase[i-1];
-
- // Check if angle cross a ref angle
- for(int a=0; a<nbIsoPhase; a++) {
- // std::cout << "a=" << rad2deg(signed_angle) << " prev=" << rad2deg(previousangle)
- // << " ref=" << rad2deg(mIsoPhaseRefAngle[a]) << " " << phase[i] << std::endl;
- if (signed_angle > previousangle) {
- // DD("cas1");
- // (((signed_angle > mIsoPhaseRefAngle[a]) && (previousangle < mIsoPhaseRefAngle[a]))) ||
- // ((mIsoPhaseRefAngle[a]==0) && (signed_angle < previousangle)))
- if ((previousangle < mIsoPhaseRefAngle[a]) &&
- (signed_angle >= mIsoPhaseRefAngle[a]))
- {
- // DD(a);
- if (nbIsoPhase == 1) { // single phase, alternate 0 and 1
- phase[i] = -phase[i-1];
- }
- else phase[i] = a;
- }
- }
- else { // previousangle >= signed_angle (we turn around 0)
- // DD("cas2");
- if ((mIsoPhaseRefAngle[a] > previousangle) ||
- (mIsoPhaseRefAngle[a] < signed_angle)) {
- // DD(a);
- if (nbIsoPhase == 1) { // single phase, alternate 0 and 1
- phase[i] = -phase[i-1];
- }
- else phase[i] = a;
- }
- }
- }
- }
- previousangle = signed_angle;
- }
-}
-//---------------------------------------------------------------------
-
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkEllipse.h,v $
- Language: C++
- Date: $Date: 2010/03/03 10:47:48 $
- Version: $Revision: 1.3 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKELLIPSE_H
-#define CLITKELLIPSE_H
-
-#include "clitkCommon.h"
-#include "clitkSignal.h"
-#include "itkVector.h"
-#include <vnl/algo/vnl_generalized_eigensystem.h>
-#include <vnl/algo/vnl_symmetric_eigensystem.h>
-#include <vnl/algo/vnl_real_eigensystem.h>
-
-namespace clitk {
-
- //---------------------------------------------------------------------
- class Ellipse : public itk::Vector<double,6> {
- public:
- typedef itk::Vector<double,2> Vector2d;
- typedef itk::Vector<double,6> Vector6d;
- typedef itk::Matrix<double,6,6> Matrix6x6d;
- typedef itk::Matrix<double,3,3> Matrix3x3d;
-
- Ellipse();
- Ellipse(const Ellipse & e);
- void InitialiseEllipseFitting(double eta, unsigned int n, clitk::Signal & inputX, clitk::Signal & inputY);
- double EllipseFittingNextIteration();
- Vector2d ComputeCenter();
- Vector2d ComputeSemiAxeLengths();
- double ComputeAngleInRad();
- void SetCenterAndSemiAxes(double x0, double y0, double r1, double r2);
- void Copy(const Vector6d & a);
- double GetEta() { return mEta; }
- void UpdateSMatrix(unsigned int begin, unsigned int n,
- clitk::Signal & inputX, clitk::Signal & inputY);
-
- protected:
- double mEta;
- double & a;
- double & b;
- double & c;
- double & d;
- double & e;
- double & f;
- Matrix6x6d C;
- Matrix6x6d Ct;
- Matrix6x6d W;
- Matrix6x6d Wt;
- Matrix6x6d S;
- Matrix6x6d Sinv;
- Matrix6x6d St;
- std::vector<Vector6d> z;
- clitk::Signal * mInputX;
- clitk::Signal * mInputY;
-
- void CopyBlock(Matrix6x6d & out, const Matrix6x6d & in,
- int ox, int oy, int l, double factor=1.0);
- void CopyBlock(Matrix6x6d & out, const Matrix3x3d & in,
- int ox, int oy, double factor=1.0);
- Matrix3x3d GetBlock3x3(const Matrix6x6d & M, int x, int y);
-
- };
- //---------------------------------------------------------------------
-
-
- void ComputeIsoPhaseFromListOfEllipses(std::vector<clitk::Ellipse*> & l,
- clitk::Signal & sx,
- clitk::Signal & sy,
- int nbIsoPhase,
- int delay,
- int L,
- clitk::Signal & phase);
-
-
-} // end namespace
-
-#endif
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalApparentMotionTracking.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 12:47:31 $
- Version: $Revision: 1.1 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKSIGNALAPPARENTMOTIONTRACKING_CXX
-#define CLITKSIGNALAPPARENTMOTIONTRACKING_CXX
-/**
- =================================================
- * @file clitkSignalApparentMotionTracking.cxx
- * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
- * @date Sept 2009
- *
- =================================================*/
-
-// clitk include
-#include "clitkSignalApparentMotionTracking_ggo.h"
-#include "clitkSignalApparentMotionTrackingFilter.h"
-#include "clitkCommon.h"
-
-int main(int argc, char * argv[]) {
-
- // Init command line
- GGO(clitkSignalApparentMotionTracking,args_info);
- CLITK_INIT;
-
- // Init filter
- clitk::SignalApparentMotionTrackingFilter filter;
- filter.SetParameters(args_info);
-
- // Run
- filter.Update();
-
- // This is the end my friend
- return EXIT_SUCCESS;
-}
-
-#endif
+++ /dev/null
-#File clitkSignalApparentMotionTrackingFilter.ggo
-package "clitkSignalMeanPositionTracking"
-version "1.0"
-purpose "See Ruan 2008, compute mean position from a signal"
-
-option "config" - "Config file" string no
-option "input" i "Input trajectory filename" string yes
-option "ref" r "Input reference trajectory filename" string yes
-
-option "refsampling" - "Sampling between ref samples" double yes
-option "inputsampling" - "Sampling between input samples" double yes
-
-option "refphase" p "Input phase file" string yes
-option "delay" d "Delay for augmented space" int no default="5"
-option "windowLength" L "Sliding window length for ellipse fit" int no default="25"
-option "nbiso" n "Number of iso phase" int no default="4"
-
-option "output" o "Output phase (debug)" string no
-
-
-
-
-option "eta" e "Convergence param (auto if not given)" double no
-option "iter" - "Max # of iterations" int no default="100"
-option "t" t "DEBUG" int no default="1000"
-
-option "phase" - "Input phase file for validation" string no
-option "delta" - "Input learn delta phase file" string no
-
-
-section "Output"
-option "residual" - "Output optimisation residual filemane" string no
-option "augmented" a "Output augmented signal filemane" string no
-
-section "Verbose"
-option "verbose" v "Verbose" flag off
-option "verbose_iteration" - "Verbose each iteration" flag off
-
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalApparentMotionTrackingFilter.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 12:41:27 $
- Version: $Revision: 1.1 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#include "clitkSignalApparentMotionTrackingFilter.h"
-#include <limits>
-#include <fstream>
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::SetParameters(args_info_clitkSignalApparentMotionTracking & a) {
- args_info = a;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::ComputeMeanAndIsoPhase(clitk::Signal & s,
- int delay,
- int L,
- int nbIsoPhase,
- clitk::Signal & mean,
- clitk::Signal & isoPhase) {
- int mMaxIteration = 1;
-
- // Compute augmented signal
- clitk::Signal sx;
- clitk::Signal sy;
- s.ComputeAugmentedSpace(sx, sy, delay);
-
- // DEBUG augmented space
- static int tt=0;
- ofstream os;
- if (tt==0) os.open("augmentedU.sig");
- else os.open("augmentedV.sig");
- tt++;
- for(uint i=0; i<sx.size(); i++) {
- os << i << " " << sx[i] << " " << sy[i] << std::endl;
- }
- os.close();
-
- // Compute starting ellipse
- Ellipse An;
- An.InitialiseEllipseFitting(-1, L, sx, sy); // -1 is auto eta
- DD(An.ComputeSemiAxeLengths());
-
- // Fit successive ellipses
- std::vector<clitk::Ellipse*> e;
- int n = s.size()-L-delay;
- for(int t=0; t<n; t++) {
- // An.InitialiseEllipseFitting(-1, L, sx, sy); // -1 is auto eta
- An.UpdateSMatrix(t, L, sx, sy);
- // DD(An.GetEta());
-// DD(An);
-// DD(An.ComputeSemiAxeLengths());
- int nn = 0;
- while (nn<mMaxIteration) {
- An.EllipseFittingNextIteration();
- // DD(nn);
-// DD(An);
-// DD(An.ComputeSemiAxeLengths());
- nn++;
- }
- clitk::Ellipse * B = new clitk::Ellipse(An);
- e.push_back(B);
- DD(t);
- DD(B->ComputeCenter());
-// DD(B->ComputeSemiAxeLengths());
- }
-
- // Get mean and isoPhase
- //mean.resize(n);
- isoPhase.resize(n);
- ComputeIsoPhaseFromListOfEllipses(e, sx, sy, nbIsoPhase, delay, L, isoPhase);
- for(unsigned int i=0; i<e.size(); i++) {
- clitk::Ellipse & An = *e[i];
- mean[i+delay+L] = An.ComputeCenter()[0];
- // DD(An.ComputeCenter());
- }
- DD("------------------------");
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::ComputeMeanAndIsoPhase(clitk::Trajectory2D & s,
- int delay,
- int L,
- int nbIsoPhase,
- clitk::Trajectory2D & mean,
- std::vector<int> & isoPhase) {
- // Extract signal (slow I know)
- clitk::Signal U(s.GetU());
- clitk::Signal V(s.GetV());
-
- clitk::Signal & meanU = mean.GetU();
- clitk::Signal & meanV = mean.GetV();
-
- meanU.resize(s.size());
- meanV.resize(s.size());
-
- clitk::Signal isoPhaseU;
- clitk::Signal isoPhaseV;
-
- DD(U.size());
- DD(V.size());
-
- // Compute mean and iso independently
- ComputeMeanAndIsoPhase(U, delay, L, nbIsoPhase, meanU, isoPhaseU);
- ComputeMeanAndIsoPhase(V, delay, L, nbIsoPhase, meanV, isoPhaseV);
-
- // Check isphase ??? ???
- isoPhase.resize(s.size());
- for(unsigned int i=0; i<isoPhaseU.size(); i++) {
-
- // ********************* U OR V ??????????
-
- isoPhase[i+delay+L] = isoPhaseV[i];
- //isoPhase[i+delay+L] = isoPhaseU[i];
-
- // mean[i+delay+L] = mean[i];
- //DD(isoPhase[i]);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::ReadDB(std::string filename,
- std::vector<Trajectory2D> & bd) {
- // File format : 2 columns by point (U,V), time by line
- std::ifstream is(filename.c_str());
- skipComment(is);
- std::string line;
- std::getline(is, line);
-
- // Compute nb of column in this first line
- std::istringstream iss(line);
- std::string s;
- int nb = 0;
- while (iss>>s) {
- nb++;
- }
- DD(nb);
-
- if (nb%2 == 1) {
- std::cout << "ERROR in the file '" << filename << "', I read " << nb
- << " columns, which is not odd." << std::endl;
- }
- nb = nb/2;
- DD(nb);
-
- // Read trajectories
- bd.resize(nb);
- while (is) {
- std::istringstream iss(line);
- std::string s;
- int i=0;
- while (iss>>s) {
- double u = atof(s.c_str());
- iss >> s;
- double v = atof(s.c_str());
- bd[i].push_back(u,v);
- i++;
- }
- std::getline(is, line);
- }
- DD(bd[0].size());
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::Update() {
-
- // -----------------------------------------------------
- // Read and build initial DB
- ReadDB(args_info.ref_arg, mReferenceTrajectoriesBD);
- mReferencePhase.Read(args_info.refphase_arg);
- mReferencePhase.SetSamplingPeriod(args_info.refsampling_arg);
-
- DD(mReferencePhase.GetSamplingPeriod());
- DD(mReferencePhase.size());
- DD(mReferenceTrajectoriesBD.size());
-
- mReferenceMean.resize(mReferenceTrajectoriesBD.size());
- for(unsigned int i=0; i<mReferenceMean.size(); i++) {
- mReferenceTrajectoriesBD[i].GetMean(mReferenceMean[i]);
- // DD(mReferenceMean[i]);
- }
-
- for(uint i=0; i<mReferenceTrajectoriesBD.size(); i++)
- mReferenceTrajectoriesBD[i].SetSamplingPeriod(mReferencePhase.GetSamplingPeriod());
-
- // Relative reference motion
- for(uint i=0; i<mReferenceTrajectoriesBD.size(); i++) {
- // DD(i);
- // mReferenceTrajectoriesBD[i].Print();
- mReferenceTrajectoriesBD[i].Substract(mReferenceMean[i]);
- // mReferenceTrajectoriesBD[i].Print();
- }
-
- // -----------------------------------------------------
- // Read trajectory
- mInput.Read(args_info.input_arg);
- DD(mInput.size());
- mInput.SetSamplingPeriod(args_info.inputsampling_arg);
- mInputMean.SetSamplingPeriod(args_info.inputsampling_arg);
- mInputRelative.SetSamplingPeriod(args_info.inputsampling_arg);
- DD(mInput.GetSamplingPeriod());
-
- // -----------------------------------------------------
- // Compute input mean and isophase
- std::cout << "Compute Mean and IsoPhase" << std::endl;
- std::vector<int> isophase;
- double d = args_info.delay_arg;
- int L = args_info.windowLength_arg;
- int nbiso = args_info.nbiso_arg;
- isophase.resize(mInput.size());
- ComputeMeanAndIsoPhase(mInput, d, L, nbiso, mInputMean, isophase);
-
- // -----------------------------------------------------
- // Substract mean to trajectory
- std::cout << "Build relative traj" << std::endl;
- mInputRelative.resize(mInput.size());
- for(int i=0; i<mInput.size(); i++) {
- // DD(i);
-// DD(mInput.GetU(i));
-// DD(mInput.GetV(i));
-// DD(mInputMean.GetU(i));
-// DD(mInputMean.GetV(i));
- mInputRelative.SetPoint(i, mInput.GetU(i) - mInputMean.GetU(i),
- mInput.GetV(i) - mInputMean.GetV(i));
- // DD(mInputRelative.GetU(i));
-// DD(mInputRelative.GetV(i));
- }
-
- // DEBUG : output mean
- ofstream osm("mean-phase.sig");
- for(int i=0; i<mInput.size(); i++) {
- osm << mInputMean.GetU(i) << " " << mInputMean.GetV(i) << " " << isophase[i] << std::endl;
- }
- osm.close();
-
- // -----------------------------------------------------
- // Loop on each iso-phase segment
- std::cout << "Loop on each isophase" << std::endl;
- clitk::Signal mEstimatedPhase;
- mEstimatedPhase.resize(mInput.size());
-
- std::vector<int> isophaseIndex;
- int previous = isophase[0];
- for(uint i=0; i<isophase.size(); i++) {
- if (isophase[i] != previous) {
- isophaseIndex.push_back(i);
- previous = isophase[i];
- }
- }
- DDV(isophaseIndex, isophaseIndex.size());
-
- for(uint ips=0; ips<isophaseIndex.size()-1; ips++) {
- DD(ips);
- int begin = isophaseIndex[ips];
- int end = isophaseIndex[ips+1];
- DD(begin);
- DD(end);
-
-
- std::cout << "-----------------------------------------------" << std::endl;
- std::cout << "-----------------------------------------------" << std::endl;
- std::cout << "Find traj from " << begin << " to " << end << std::endl;
- mInput.Print("mInput", begin, end);
- mInputRelative.Print("mInputRelativ", begin, end);
-
-
- // Get current mean position
- Vector2d currentMeanPosition;
- mInputMean.GetPoint(begin, currentMeanPosition);
- DD(currentMeanPosition);
-
- // Find closest ref trajectory
- int index = FindClosestTrajectory(mReferenceMean, currentMeanPosition);
- DD(index);
- clitk::Trajectory2D & currentClosestReferenceTrajectory = mReferenceTrajectoriesBD[index];
-
- DD(mReferenceMean[index]);
- currentClosestReferenceTrajectory.Print("currentClosestReferenceTrajectory");
-
- // Estimate phase delta
- clitk::Signal tempPhase;
- EstimatePhase(mInputRelative, begin, end,
- currentClosestReferenceTrajectory,
- mReferencePhase, tempPhase);
-
- // Cat current phase
- int n=0;
- for(int i=begin; i<end; i++) {
- mEstimatedPhase[i] = tempPhase[n];
- n++;
- }
- }
-
- // -----------------------------------------------------
- // Output time - phase
- // DEBUG OUTPUT
- ofstream os(args_info.output_arg);
- for(int t=0; t<mInput.size(); t++) {
- os << t << " " << t*mInput.GetSamplingPeriod() << " " << mEstimatedPhase[t] << std::endl;
- }
- os.close();
-
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalApparentMotionTrackingFilter::EstimatePhase(const clitk::Trajectory2D & input,
- const int begin, const int end,
- const clitk::Trajectory2D & ref,
- const clitk::Signal & phaseref,
- clitk::Signal & phase) {
-
- // Create time-warped resampled trajectory
- clitk::Trajectory2D T;
- T.SetSamplingPeriod(ref.GetSamplingPeriod());
- DD(input.size());
- input.ResampleWithTimeWarpTo(begin, end, T);
- DD(input.GetSamplingPeriod());
- DD(ref.GetSamplingPeriod());
- ref.Print("ref");
- input.Print("input", begin, end);
- T.Print("T");
-
- // Find optimal delta
- int delta = FindPhaseDelta(T, ref);
- DD(delta);
-
- // Shift phase
- clitk::Signal temp(phaseref);
- DD(temp.GetSamplingPeriod());
- DDV(temp, temp.size());
- temp.Shift(delta, (int)lrint((end-begin)*input.GetSamplingPeriod()/phaseref.GetSamplingPeriod()));
- DD(temp.GetSamplingPeriod());
- DDV(temp, temp.size());
-
-
- // Output result with time-unwarp
- phase.resize(end-begin);
- phase.SetSamplingPeriod(input.GetSamplingPeriod());
- DD(phase.GetSamplingPeriod());
- DD(phaseref.GetSamplingPeriod());
- DD(temp.GetSamplingPeriod());
- DD(temp.size());
- DD(phase.size());
- temp.ResampleWithTimeWarpTo(phase, true); // linear
- DDV(phase, phase.size());
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-int clitk::SignalApparentMotionTrackingFilter::FindPhaseDelta(const clitk::Trajectory2D & A,
- const clitk::Trajectory2D & B) const {
- assert(A.GetSamplingPeriod() == B.GetSamplingPeriod());
-
- DD("A");
- A.Print("A");
- DD("B");
- B.Print("B");
-
- // exhaustive search
- double min=numeric_limits<double>::max();
- int imin =0;
- for(int i=0; i<B.size(); i++) {
- double d = A.DistanceTo(i, B);
- std::cout.precision(10);
- std::cout << "i=" << i << " d=" << d << " (min=" << min << " imin=" << imin << ")" << std::endl;
- if (d < min) {
- imin = i;
- min = d;
- }
- }
- DD(imin);
- return imin;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-int clitk::SignalApparentMotionTrackingFilter::FindClosestTrajectory(const std::vector<Vector2d> & meanlist,
- const Vector2d & mean) {
- // Brute force for the moment
- double dmin = numeric_limits<double>::max();
- double imin=0;
- for(uint i=0; i<meanlist.size(); i++) {
- double d = (meanlist[i]-mean).GetNorm();
- if (d< dmin) {
- dmin = d;
- imin = i;
- }
- }
- DD(dmin);
- DD(imin);
- return imin;
-}
-//---------------------------------------------------------------------
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalApparentMotionTrackingFilter.h,v $
- Language: C++
- Date: $Date: 2010/03/03 12:41:27 $
- Version: $Revision: 1.1 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef SIGNALAPPARENTMOTIONTRACKINGFILTER_H
-#define SIGNALAPPARENTMOTIONTRACKINGFILTER_H
-
-#include "clitkSignalApparentMotionTracking_ggo.h"
-#include "clitkSignal.h"
-#include "clitkEllipse.h"
-#include "clitkTrajectory2D.h"
-#include "itkVector.h"
-#include <math.h>
-
-namespace clitk {
-
- //---------------------------------------------------------------------
- class SignalApparentMotionTrackingFilter {
- public:
- typedef itk::Vector<double,2> Vector2d;
-
- void SetParameters(args_info_clitkSignalApparentMotionTracking & args_info);
- void Update();
-
- protected:
- args_info_clitkSignalApparentMotionTracking args_info;
- clitk::Signal mRefU;
- clitk::Signal mRefV;
- clitk::Signal mRefPhase;
- clitk::Signal mInputU;
- clitk::Signal mInputV;
- clitk::Signal mMeanU;
- clitk::Signal mMeanV;
- clitk::Signal mIsoPhaseU;
- clitk::Signal mIsoPhaseV;
- double mRefSampling;
-
- std::vector<Trajectory2D> mReferenceTrajectoriesBD;
- clitk::Signal mReferencePhase;
- clitk::Trajectory2D mInput;
- clitk::Trajectory2D mInputRelative;
- clitk::Trajectory2D mInputMean;
- std::vector<Vector2d> mReferenceMean;
-
- void ComputeMeanAndIsoPhase(clitk::Signal & s, int delay, int L, int nbIso,
- clitk::Signal & mean, clitk::Signal & isoPhase);
- void ComputeMeanAndIsoPhase(clitk::Trajectory2D & s, int delay, int L,
- int nbIsoPhase, clitk::Trajectory2D & mean,
- std::vector<int> & isoPhase);
- void EstimatePhase(const clitk::Trajectory2D & input,
- const int begin, const int end,
- const clitk::Trajectory2D & ref,
- const clitk::Signal & phaseref,
- clitk::Signal & phase);
- int FindPhaseDelta(const clitk::Trajectory2D & A, const clitk::Trajectory2D & B) const;
- int FindClosestTrajectory(const std::vector<Vector2d> & meanlist, const Vector2d & mean);
- void ReadDB(std::string filename, std::vector<Trajectory2D> & bd);
- };
- //---------------------------------------------------------------------
-
-} // end namespace
-
-#endif
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Language: C++
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKSIGNALFILTER_CXX
-#define CLITKSIGNALFILTER_CXX
-/**
- =================================================
- * @file clitkSignalFilter.cxx
- * @author Jef Vandemeulebroucke <jefvdmb@gmail.com>
- * @date 23 April 2008
- *
- * @brief "Apply a filter to the signal"
- =================================================*/
-
-// clitk include
-#include "clitkSignalFilter_ggo.h"
-#include "clitkSignal.h"
-#include "clitkIO.h"
-#include "clitkImageCommon.h"
-
-int main(int argc, char * argv[]) {
-
- // Init command line
- GGO(clitkSignalFilter, args_info);
- CLITK_INIT;
-
- //input
- clitk::Signal input;
- input.Read (args_info.input_arg);
- clitk::Signal input2;
- if (args_info.input2_given) input2.Read(args_info.input2_arg);
-
- double p1 = args_info.p1_arg;
- double p2 = args_info.p2_arg;
-
- //process
-
- if( args_info.multiply_flag) input*=input2;
- if( args_info.divide_flag) input/=input2;
- if( args_info.norm_flag) input=input.Normalize (p1, p2);
- //if( args_info.highPass_flag) input=input.HighPassFilter(p1,p2);
- //if( args_info.lowPass_flag) input=input.LowPassFilter(p1,p2);
- if( args_info.detect_flag) input=input.DetectLocalExtrema(static_cast<unsigned int >(p1));
- if( args_info.limPhase_flag) input=input.LimPhase();
- if( args_info.monPhase_flag) input=input.MonPhase();
- if( args_info.monPhaseDE_flag) input=input.MonPhaseDE(p1,p2);
- if( args_info.average_flag) input=input.MovingAverageFilter(static_cast<unsigned int> (p1));
- if( args_info.ssd_flag) std::cout<<"The sqrt of the mean SSD is "<< input.SSD(input2)<<std::endl;
- if (args_info.gauss_flag)input=input.GaussLikeFilter();
- if (args_info.rescale_flag)input=input.NormalizeMeanStdDev(p1, p2);
- if (args_info.interp_flag)input=input.LinearlyInterpolateScatteredValues();
- // if( args_info.approx_flag) input=input.ApproximateScatteredValuesWithBSplines(static_cast<unsigned int>(p1),static_cast<unsigned int>(p2));
- if( args_info.limit_flag) input=input.LimitSignalRange();
-// if( args_info._flag) input=input;
-// if( args_info._flag) input=input;
-// if( args_info._flag) input=input;
-// if( args_info._flag) input=input;
- // if( args_info._flag) input=input;
- // if( args_info._flag) input=input;
- // if( args_info._flag) input=input;
-
- if (args_info.output_given) input.Write(args_info.output_arg);
-
- return EXIT_SUCCESS;
-
-}
-
-#endif
+++ /dev/null
-#File clitkSignalFilter.ggo
-package "clitk"
-version "Pass an filter on an signal"
-
-option "config" - "Config file" string no
-option "verbose" v "Verbose" flag off
-option "input" i "Input grid filename" string yes
-option "input2" j "Second Input grid filename" string no
-option "p1" - "p1" double no
-option "p2" - "p2" double no
-option "output" o "Output grid filename" string no
-option "multiply" - "Mulitply values of input 1 with input 2 (same size)" flag off
-option "divide" - "Divide values of input 1 with input 2 (same size, zeros are skipped)" flag off
-option "norm" n "Normalize signal between p1=min and p2=max" flag off
-option "rescale" - "Rescale signal to p1=mean and p2=StdDev" flag off
-#option "highPass" - "High Pass Filter: p1=sampPeriod, p2= cutOffFrequecy" flag off
-#option "lowPass" - "Low Pass Filter: p1=sampPeriod, p2= cutOffFrequecy" flag off
-option "detect" - "Detect Local extrema: local window= 2*p1+1, p1>=1" flag off
-option "limPhase" - "Convert extrema signal into limited phase [0, 1[" flag off
-option "monPhase" - "Convert extrema signal into monotone phase [0, inf[" flag off
-option "monPhaseDE" - "Convert extrema signal into scattered monotone phase, taking into account both extrema (p1=eEPhaseValue,p2=eIPhaseValue)" flag off
-option "average" - "Moving Average filter: window= 2*p1+1 [0, inf[" flag off
-option "ssd" - "Calculate the sqrt of the mean SSD between i and j" flag off
-option "gauss" - "Pass a small gauss-like kernel (121-mask)" flag off
-option "interp" - "Linearly interpolate scattered values (unknow values=-1)" flag off
-option "approx" - "Bspline approximate scattered values(unknow values=-1, p1=Order, p2=number of control points" flag off
-option "limit" - "Limit the range to [0, 1[ (signal-floor (signal))" flag off
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalMeanPositionFilter.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 12:51:06 $
- Version: $Revision: 1.7 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#include "clitkSignalMeanPositionFilter.h"
-#include <fstream>
-
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPositionTracking & args) {
- args_info = args;
- mEtaIsSet = false;
- mOutputFilenameIsSet = false;
- mOutputResidualFilenameIsSet = false;
- mInput.Read(args_info.input_arg);
- mAugmentationDelay = args_info.delay_arg;
- mIsAdaptiveMethod = false;
- mWindowLength = -1;
- mValidationWithRealPhase = false;
- mUseLearnedDeltaPhase = false;
- if (args_info.eta_given) {
- mEta = args_info.eta_arg;
- mEtaIsSet = true;
- }
- mMaxIteration = args_info.iter_arg;
- if (args_info.output_given) {
- mOutputFilenameIsSet = true;
- mOutputFilename = args_info.output_arg;
- }
- if (args_info.residual_given) {
- mOutputResidualFilenameIsSet = true;
- mOutputResidualFilename = args_info.residual_arg;
- }
- if (args_info.augmented_given) {
- mOutputAugmentedFilenameIsSet = true;
- mOutputAugmentedFilename = args_info.augmented_arg;
- }
- mVerbose = args_info.verbose_flag;
- mVerboseIteration = args_info.verbose_iteration_flag;
- if (args_info.L_given) {
- mIsAdaptiveMethod = true;
- mWindowLength = args_info.L_arg;
- }
- if (args_info.phase_given) {
- mValidationWithRealPhase = true;
- mInputPhaseFilename = args_info.phase_arg;
- mInputPhase.Read(mInputPhaseFilename);
- }
- mNumberOfIsoPhase = args_info.nbiso_arg;
- if (args_info.delta_given) {
- mUseLearnedDeltaPhase = true;
- mLearnIsoPhaseDelta.Read(args_info.delta_arg);
- mNumberOfIsoPhase = mLearnIsoPhaseDelta.size();
- }
-
- // DEBUG
- if (args_info.phase_given) {
- std::cout << "PLEASE DO NO USE THIS OPTION --phase YET ..." << std::endl;
- exit(0);
- }
- if (args_info.delta_given) {
- std::cout << "PLEASE DO NO USE THIS OPTION --delta YET ..." << std::endl;
- exit(0);
- }
-
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::Update() {
-
- // Compute augmented space
- if (mVerbose) {
- std::cout << "Computing augmented space with delay = " << mAugmentationDelay
- << ", number of samples is " << mInput.size()-mAugmentationDelay
- << std::endl;
- }
- mInput.ComputeAugmentedSpace(mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
-
- // Output augmented state
- if (mOutputAugmentedFilenameIsSet) {
- std::ofstream os;
- openFileForWriting(os, mOutputAugmentedFilename);
- for(unsigned int i=0; i<(unsigned int)mAugmentedInputX.size(); i++) {
- os << mAugmentedInputX[i] << " " << mAugmentedInputY[i] << std::endl;
- }
- os.close();
- }
-
- // Initial starting ellipse
- Ellipse An; // starting point
- if (!mEtaIsSet) mEta = -1;
- DD(mWindowLength);
- if (mWindowLength == -1) {
- mWindowLength = mAugmentedInputY.size();
- }
- An.InitialiseEllipseFitting(mEta, mWindowLength, mAugmentedInputX, mAugmentedInputY);
- mEta = An.GetEta();
- if (mVerbose) {
- std::cout << "Eta is " << mEta << std::endl;
- }
- Ellipse AnInitial(An);
- if (mVerbose) {
- std::cout << "Set initial ellipse to c=" << An.ComputeCenter()
- << " and axes=" << An.ComputeSemiAxeLengths() << std::endl;
- }
-
- // Fitting method
- if (mIsAdaptiveMethod) {
- AdaptiveFitEllipse(An);
- }
- else FitEllipse(An);
-
- // Output
- if (mVerbose) {
- std::cout << "Result is " << An << std::endl;
- std::cout << "Center is " << An.ComputeCenter() << std::endl;
- std::cout << "SemiAxis are " << An.ComputeSemiAxeLengths() << std::endl;
- std::cout << "Angle is " << rad2deg(An.ComputeAngleInRad()) << " deg" << std::endl;
- }
- if (mOutputFilenameIsSet) {
- std::ofstream os;
- openFileForWriting(os, mOutputFilename);
-
- std::vector<double> phase;
- ComputeIsoPhase(mListOfEllipses, phase, mCycles);
-
- //int currentCycle = 0;
- for (unsigned int i=0; i<mListOfEllipses.size(); i++) {
- clitk::Ellipse & An = *mListOfEllipses[i];
- int time = i+mAugmentationDelay+mWindowLength;
- os << time << " " // current 'time' in input
- << mInput[time] << " " // current value in input
- << An.ComputeCenter()[0] << " " << An.ComputeCenter()[1] << " "
- << An.ComputeSemiAxeLengths()[0] << " " << An.ComputeSemiAxeLengths()[1] << " "
- << An.ComputeAngleInRad() << " "
- << rad2deg(An.ComputeAngleInRad()) << " _Phase_ "
- << mIsoPhaseIndex[i] << " " ;
-
- /*
- if (mUseLearnedDeltaPhase) {
- os << mIsoPhaseDelta[mIsoPhaseIndex[i]] << " " << phase[i] << " ";
- }
- */
-
- os << " _Ellipse_ ";
- for(int j=0; j<6; j++) os << An[j] << " ";
- os << std::endl;
- }
- os.close();
- }
-
- if (mOutputResidualFilenameIsSet) {
- std::ofstream os;
- openFileForWriting(os, mOutputResidualFilename);
- for(unsigned int i=0; i<mCurrentResidual.size(); i++) {
- os.precision(10);
- os << mCurrentResidual[i] << std::endl;
- }
- os.close();
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
- int n = 0;
- mCurrentResidual.clear();
- while (n<mMaxIteration) {
- double r = An.EllipseFittingNextIteration();
- mCurrentResidual.push_back(r);
- n++;
- if (mVerboseIteration) {
- std::cout.precision(3);
- std::cout << "A(" << n << ")=" << An
- << " c=" << An.ComputeCenter()
- << " ax=" << An.ComputeSemiAxeLengths()
- << " theta=" << rad2deg(An.ComputeAngleInRad());
- std::cout.precision(10);
- std::cout << " R=" << r << std::endl;
- }
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::AdaptiveFitEllipse(clitk::Ellipse & An) {
- // Compute the total number of iteration :
- // = nb of inputs minus the delay, minus the windowL, minus the first one
- int maxN = std::min((unsigned int)args_info.t_arg,
- mInput.size()-mWindowLength-mAugmentationDelay);
- for(int t=0; t<maxN; t++) {
- An.UpdateSMatrix(t, mWindowLength, mAugmentedInputX, mAugmentedInputY);
- FitEllipse(An);
- clitk::Ellipse * B = new clitk::Ellipse(An);
- mListOfEllipses.push_back(B);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::ComputeIsoPhase(std::vector<clitk::Ellipse*> & l,
- std::vector<double> & phase,
- std::vector<int> & cycles) {
- double refphaseangle=0;
- double previousangle=0;
- phase.resize(mListOfEllipses.size());
-
- // DD(mListOfEllipses.size());
- mIsoPhaseIndex.resize(mListOfEllipses.size());
- // mIsoPhaseDelta.resize(mNumberOfIsoPhase);
- // mIsoPhaseDeltaNb.resize(mNumberOfIsoPhase);
- mIsoPhaseRefAngle.resize(mNumberOfIsoPhase);
-
- for (unsigned int i=0; i<mListOfEllipses.size(); i++) {
- clitk::Ellipse & An = *mListOfEllipses[i];
-
- // Current time accordint to initial input signal
- // int time = i+mAugmentationDelay+mWindowLength; // not use yet
-
- // Compute current signed angle
- Vector2d x1(An.ComputeCenter());
- double a = mListOfEllipses[0]->ComputeSemiAxeLengths()[0];
- double theta = mListOfEllipses[0]->ComputeAngleInRad();
- Vector2d x2; x2[0] = x1[0]+a * cos(theta); x2[1] = x1[1]+a * sin(theta);
- Vector2d x3(x1);
- Vector2d x4; x4[0] = mAugmentedInputX[i+mWindowLength]; x4[1] = mAugmentedInputY[i+mWindowLength];
- Vector2d A(x2-x1);
- Vector2d B(x4-x3);
- double signed_angle = atan2(B[1], B[0]) - atan2(A[1], A[0]);
- // double signed_angle = atan2(B[1], B[0]) - atan2(0, 1);
- if (signed_angle<0) signed_angle = 2*M_PI+signed_angle;
- // http://www.euclideanspace.com/maths/algebra/vectors/angleBetween/index.htm
-
- // First time : set the angle
- if (i==0) {
- refphaseangle = signed_angle;
- for(int a=0; a<mNumberOfIsoPhase; a++) {
- if (a==0) mIsoPhaseRefAngle[a] = signed_angle;
- else mIsoPhaseRefAngle[a] = (signed_angle
- + (2*M_PI)*a/mNumberOfIsoPhase);
- if (mIsoPhaseRefAngle[a] > 2*M_PI) mIsoPhaseRefAngle[a] -= 2*M_PI;
- if (mIsoPhaseRefAngle[a] < 0) mIsoPhaseRefAngle[a] = 2*M_PI-mIsoPhaseRefAngle[a];
- // DD(rad2deg(mIsoPhaseRefAngle[a]));
- // mIsoPhaseDelta[a] = 0.0;
- // mIsoPhaseDeltaNb[a] = 0;
- }
- int a=0;
- // DD(rad2deg(signed_angle));
- while ((a<mNumberOfIsoPhase) && (signed_angle>=mIsoPhaseRefAngle[a])) { a++; }
- // DD(a);
- mIsoPhaseIndex[i] = a-1;
- // DD(mIsoPhaseIndex[0]);
- cycles.push_back(0);
- }
- else {
- mIsoPhaseIndex[i] = mIsoPhaseIndex[i-1];
-
- // Check if angle cross a ref angle
- for(int a=0; a<mNumberOfIsoPhase; a++) {
- std::cout << "a=" << rad2deg(signed_angle) << " p=" << rad2deg(previousangle)
- << " ref=" << rad2deg(mIsoPhaseRefAngle[a]) << std::endl;
- if (
- (((signed_angle > mIsoPhaseRefAngle[a]) && (previousangle < mIsoPhaseRefAngle[a]))) ||
- ((mIsoPhaseRefAngle[a]==0) && (signed_angle < previousangle)))
- {
- // if (mValidationWithRealPhase) {
- // mIsoPhaseDelta[a] += mInputPhase[time];
- // mIsoPhaseDeltaNb[a]++;
- // }
- mIsoPhaseIndex[i] = a;
- // DD(a);
- cycles.push_back(i);
- }
- }
- // DD(mIsoPhaseIndex[i]);
- }
-
- previousangle = signed_angle;
- }
-
- /*
- if (mValidationWithRealPhase) {
- // Mean of all deltas
- for(unsigned int a=0; a<mIsoPhaseDelta.size(); a++) {
- if (mIsoPhaseDelta[a] == 0) {
- std::cerr << "Error : less than one cyle ?" << std::endl;
- exit(0);
- }
- mIsoPhaseDelta[a] /= mIsoPhaseDeltaNb[a];
- DD(mIsoPhaseDelta[a]);
- }
- std::ofstream os;
- openFileForWriting(os, "delta.sig");
- for(unsigned int a=0; a<mIsoPhaseDelta.size(); a++) {
- os << mIsoPhaseDelta[a] << std::endl;
- }
- os.close();
- }
- */
-
- /*
- if (mUseLearnedDeltaPhase) {
- for(unsigned int a=0; a<mIsoPhaseDelta.size(); a++) {
- mIsoPhaseDelta[a] = mLearnIsoPhaseDelta[a];
- DD(mIsoPhaseDelta[a]);
- }
- }
- */
-
- // DEBUG UNUSED
- /*
- int j=0;
- for (unsigned int i=0; i<cycles.size(); i++) {
- DD(cycles[i]);
- int n;
- if (i!=0) {
- n=cycles[i]-cycles[i-1];
- for(int index=0; index<n; index++) {
-
- int indexOfPhase1 = mIsoPhaseIndex[cycles[i-1]];
- int indexOfPhase2 = mIsoPhaseIndex[cycles[i]];
- DD(indexOfPhase1);
- DD(indexOfPhase2);
- double ph1 = mIsoPhaseDelta[indexOfPhase1];
- double ph2 = mIsoPhaseDelta[indexOfPhase2];
- DD(ph1);
- DD(ph2);
-
- phase[j] = (ph1+(double)index*(ph2-ph1)/(double)n);
- if (phase[j]>1.0) phase[j] = 1.0-phase[j];
- DD(phase[j]);
- j++;
- }
- }
- else j+=cycles[0];
- }
- */
- for (unsigned int i=0; i<cycles.size(); i++) {
- DD(cycles[i]);
- }
-
-}
-//---------------------------------------------------------------------
-
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalMeanPositionFilter.h,v $
- Language: C++
- Date: $Date: 2010/03/03 10:47:48 $
- Version: $Revision: 1.6 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKSIGNALMEANPOSITIONFILTER_H
-#define CLITKSIGNALMEANPOSITIONFILTER_H
-
-#include "clitkSignalMeanPositionTracking_ggo.h"
-#include "clitkSignal.h"
-#include "clitkEllipse.h"
-#include "itkVector.h"
-#include <math.h>
-
-namespace clitk {
-
- //---------------------------------------------------------------------
- class SignalMeanPositionFilter {
- public:
- typedef itk::Vector<double,2> Vector2d;
-
- void SetParameters(args_info_clitkSignalMeanPositionTracking & args_info);
- void Update();
-
- protected:
- args_info_clitkSignalMeanPositionTracking args_info;
- clitk::Signal mInput;
- clitk::Signal mAugmentedInputX;
- clitk::Signal mAugmentedInputY;
- int mAugmentationDelay;
- int mMaxIteration;
- double mEta;
- bool mEtaIsSet;
- bool mOutputFilenameIsSet;
- bool mOutputResidualFilenameIsSet;
- bool mOutputAugmentedFilenameIsSet;
- std::string mOutputFilename;
- std::string mOutputResidualFilename;
- std::string mOutputAugmentedFilename;
- bool mVerbose;
- bool mVerboseIteration;
- bool mIsAdaptiveMethod;
- std::vector<double> mCurrentResidual;
- int mWindowLength;
- std::vector<clitk::Ellipse*> mListOfEllipses;
-
- bool mValidationWithRealPhase;
- std::string mInputPhaseFilename;
- clitk::Signal mInputPhase;
- std::vector<int> mCycles;
-
- std::vector<int> mIsoPhaseIndex;
- std::vector<double> mIsoPhaseDelta;
- std::vector<int> mIsoPhaseDeltaNb;
- std::vector<double> mIsoPhaseRefAngle;
-
- bool mUseLearnedDeltaPhase;
- clitk::Signal mLearnIsoPhaseDelta;
- int mNumberOfIsoPhase;
-
- void FitEllipse(clitk::Ellipse & An);
- void AdaptiveFitEllipse(clitk::Ellipse & An);
-
- void ComputeIsoPhase(std::vector<clitk::Ellipse*> & l,
- std::vector<double> & phase,
- std::vector<int> & cycles);
- };
- //---------------------------------------------------------------------
-
-} // end namespace
-
-#endif
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkSignalMeanPositionTracking.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 10:47:48 $
- Version: $Revision: 1.6 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKSIGNALMEANPOSITIONTRACKING_CXX
-#define CLITKSIGNALMEANPOSITIONTRACKING_CXX
-/**
- =================================================
- * @file clitkSignalMeanPositionTracking.cxx
- * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
- * @date Sept 2009
- *
- * @brief "See Ruan 2007, compute mean position from a signal"
- =================================================*/
-
-// clitk include
-#include "clitkSignalMeanPositionTracking_ggo.h"
-#include "clitkCommon.h"
-#include "clitkSignalMeanPositionFilter.h"
-
-int main(int argc, char * argv[]) {
-
- // Init command line
- GGO(clitkSignalMeanPositionTracking,args_info);
- CLITK_INIT;
-
- // Init filter
- clitk::SignalMeanPositionFilter filter;
- filter.SetParameters(args_info);
-
- // Run
- filter.Update();
-
- // This is the end my friend
- return EXIT_SUCCESS;
-}
-
-#endif
+++ /dev/null
-#File clitkSignalMeanPositionTracking.ggo
-package "clitkSignalMeanPositionTracking"
-version "1.0"
-purpose "See Ruan 2008, compute mean position from a signal"
-
-option "config" - "Config file" string no
-option "input" i "Input signal filename" string yes
-option "delay" d "Delay for augmented space" int no default="10"
-option "eta" e "Convergence param (auto if not given)" double no
-option "iter" - "Max # of iterations" int no default="100"
-option "L" L "Sliding window length" int no
-option "t" t "DEBUG" int no default="1000"
-
-option "phase" p "Input phase file for validation" string no
-option "delta" - "Input learn delta phase file" string no
-option "nbiso" n "Number of iso phase" int no default="1"
-
-
-section "Output"
-option "output" o "Output ellipse param filename (2 lines=initial/final ; 6 quadratic param + 5 parametric param=center/axis/theta)" string no
-option "residual" r "Output optimisation residual filemane" string no
-option "augmented" a "Output augmented signal filemane" string no
-
-section "Verbose"
-option "verbose" v "Verbose" flag off
-option "verbose_iteration" - "Verbose each iteration" flag off
-
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkTrajectory2D.cxx,v $
- Language: C++
- Date: $Date: 2010/03/03 12:41:27 $
- Version: $Revision: 1.1 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#include "clitkTrajectory2D.h"
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::ResampleWithTimeWarpTo(int begin, int end, Trajectory2D & output) const {
- DD(begin);
- DD(end);
- int n = (int)lrint((end-begin)*GetSamplingPeriod()/output.GetSamplingPeriod());
- output.resize(n);
- DD(n);
- double duration = (end-begin)*GetSamplingPeriod();
- DD(duration);
- double sp = output.GetSamplingPeriod();
- DD(output.GetTotalTimeDuration());
- // DD(sp);
- for(int i=0; i<n; i++) {
- Vector2d p;
- // DD(i);
- GetValueAtLin((i*sp)+(begin*GetSamplingPeriod()), p);
- // DD(p);
- output.SetPoint(i, p);
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-bool clitk::Trajectory2D::Read(std::string filename) {
- bool b = mU.Read(filename, 0);
- if (!b) return false;
- b = mV.Read(filename, 1);
- return b;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::GetValueAtLin(double t, Vector2d & p) const {
- // Can be faster, I know ...
- p[0] = mU.GetValueAtLin(t);
- p[1] = mV.GetValueAtLin(t);
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::ResampleWithTimeWarpTo(Trajectory2D & output) const {
- ResampleWithTimeWarpTo(0, size(), output);
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::Shift(double s) {
- std::cout << "Shift s=" << s << std::endl;
- int d = (int)lrint(s/GetSamplingPeriod()); // closest integer delta
- DD(d);
- clitk::Signal tempU((uint)size());
- clitk::Signal tempV((uint)size());
- for(int i=0; i<size(); i++) {
- tempU[0] = mU[(i+d)%size()];
- tempV[0] = mV[(i+d)%size()];
- }
- for(int i=0; i<size(); i++) {
- mU[i] = tempU[i];
- mV[i] = tempV[i];
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::GetPoint(const int index, Vector2d & p) {
- p[0] = mU[index];
- p[1] = mV[index];
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-clitk::Trajectory2D & clitk::Trajectory2D::operator+(Trajectory2D & d) {
- for(int i=0; i<size(); i++) {
- mU[i] += d.GetU(i);
- mV[i] += d.GetV(i);
- }
- return *this;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-clitk::Trajectory2D & clitk::Trajectory2D::operator-(Trajectory2D & d) {
- for(int i=0; i<size(); i++) {
- mU[i] -= d.GetU(i);
- mV[i] -= d.GetV(i);
- }
- return *this;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::GetMean(Vector2d & m) {
- m[0] = 0.0;
- m[1] = 0.0;
- for(int i=0; i<size(); i++) {
- m[0] += mU[i];
- m[1] += mV[i];
- }
- m[0] /= size();
- m[1] /= size();
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-double clitk::Trajectory2D::DistanceTo(int delta, const clitk::Trajectory2D & B) const {
- double d = 0.0;
- for(int n=0; n<size(); n++) {
- int i = n;//(n)%size();
- int j = (n+delta)%B.size();
- // DD(j);
-// DD(i);
-// DD(n);
- d += pow(mU[i] - B.GetU(j), 2) + pow(mV[i] - B.GetV(j), 2);
- }
- return d;
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::Print(const std::string & name) const {
- Print(name, 0, size());
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::Print(const std::string & name, int begin, int end) const {
- std::cout << "Traj " << name << " size = " << size()
- << " from " << begin << " to " << end << std ::endl;
- for(int i=begin; i<end; i++) {
- std::cout << mU[i] << " " << mV[i] << std::endl;
- }
-}
-//---------------------------------------------------------------------
-
-
-//---------------------------------------------------------------------
-void clitk::Trajectory2D::Substract(const Vector2d & m) {
- mU.AddValue(-m[0]);
- mV.AddValue(-m[1]);
-}
-//---------------------------------------------------------------------
+++ /dev/null
-/*=========================================================================
-
- Program: clitk
- Module: $RCSfile: clitkTrajectory2D.h,v $
- Language: C++
- Date: $Date: 2010/03/03 12:41:27 $
- Version: $Revision: 1.1 $
-
- Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
- l'Image). All rights reserved. See Doc/License.txt or
- http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-
-#ifndef CLITKTRAJECTORY2D_H
-#define CLITKTRAJECTORY2D_H
-
-#include "clitkCommon.h"
-#include "clitkSignal.h"
-#include "itkVector.h"
-#include <vnl/algo/vnl_generalized_eigensystem.h>
-#include <vnl/algo/vnl_symmetric_eigensystem.h>
-#include <vnl/algo/vnl_real_eigensystem.h>
-
-namespace clitk {
-
- //---------------------------------------------------------------------
- class Trajectory2D {
- public:
- typedef itk::Vector<double,2> Vector2d;
-
- void ResampleWithTimeWarpTo(int begin, int end, Trajectory2D & output) const;
- void ResampleWithTimeWarpTo(Trajectory2D & output) const;
- void Shift(double s);
-
- bool Read(std::string filename);
-
- void GetPoint(const int index, Vector2d & p);
- void GetMean(Vector2d & m);
- double GetU(const int index) const { return mU[index]; }
- double GetV(const int index) const { return mV[index]; }
- clitk::Signal & GetU() { return mU; }
- clitk::Signal & GetV() { return mV; }
-
- int size() const { return mU.size(); }
- void resize(int n) { mU.resize(n); mV.resize(n); }
-
- void SetPoint(const int i, Vector2d & p) { mU[i] = p[0]; mV[i] = p[1]; }
- void SetPoint(const int i, double u, double v) { mU[i] = u; mV[i] = v; }
- void GetValueAtLin(double t, Vector2d & p) const;
-
- double GetSamplingPeriod() const { return mU.GetSamplingPeriod(); }
- void SetSamplingPeriod(double sp){ mU.SetSamplingPeriod(sp); mV.SetSamplingPeriod(sp); }
- double DistanceTo(int delta, const clitk::Trajectory2D & B) const;
-
- Trajectory2D & operator+(Trajectory2D & d);
- Trajectory2D & operator-(Trajectory2D & d);
-
- void Print(const std::string & name) const;
- void Print(const std::string & name, int begin, int end) const;
-
- void push_back(double u, double v) { mU.push_back(u); mV.push_back(v); }
- void Substract(const Vector2d & m);
- double GetTotalTimeDuration() const { return mU.GetTotalTimeDuration(); }
-
- protected:
- clitk::Signal mU;
- clitk::Signal mV;
-
- };
- //---------------------------------------------------------------------
-
-
-} // end namespace
-
-#endif