From: lgrezesb <lgrezesb>
Date: Wed, 10 Feb 2010 14:55:00 +0000 (+0000)
Subject: Good version of clitkSignalMeanPosition, from clitk2
X-Git-Tag: v1.2.0~813
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=253e4eea321f0a154e1c57834f36d4e4c4b34ef1;p=clitk.git

Good version of clitkSignalMeanPosition, from clitk2
---

diff --git a/tools/clitkSignalMeanPositionFilter.cxx b/tools/clitkSignalMeanPositionFilter.cxx
index 912d9ee..905af5e 100644
--- a/tools/clitkSignalMeanPositionFilter.cxx
+++ b/tools/clitkSignalMeanPositionFilter.cxx
@@ -1,7 +1,10 @@
 /*=========================================================================
                                                                                 
   Program:   clitk
+  Module:    $RCSfile: clitkSignalMeanPositionFilter.cxx,v $
   Language:  C++
+  Date:      $Date: 2010/02/10 14:55:00 $
+  Version:   $Revision: 1.4 $
                                                                                 
   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   l'Image). All rights reserved. See Doc/License.txt or
@@ -25,6 +28,8 @@ void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPos
   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;
@@ -48,6 +53,28 @@ void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPos
     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);
+  }
+
 }
 //---------------------------------------------------------------------
 
@@ -55,45 +82,19 @@ void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPos
 //---------------------------------------------------------------------
 void clitk::SignalMeanPositionFilter::Update() {
 
-  // DEBUG
-//   int e = 5;
-//   clitk::Signal temp;
-//   temp.resize(mInput.size()*e);
-//   for(unsigned int i=0; i<temp.size(); i++) {
-//     int index = lrint(floor((double)i/e));
-//     double res = (double)i/e - index;
-//     //    DD(index);
-//     //DD(res);
-//     temp[i] = mInput[index] + res*(mInput[index+1]-mInput[index])/e;
-//   }
-//   mInput.resize(temp.size());
-  for  (unsigned int i=0; i<mInput.size(); i++) mInput[i] = mInput[i]+(double)i/(double)mInput.size();
-  srand ( time(0) );
-  for  (unsigned int i=0; i<mInput.size(); i++) {
-    mInput[i] = (0.8+((double)rand()/ (double)RAND_MAX)*0.4) * mInput[i];
-  }
-
-   {
-     std::ofstream os;
-     openFileForWriting(os, "input.dat");
-     for(unsigned int i=0; i<mInput.size(); i++) os << mInput[i] << std::endl;
-     os.close();
-   }
-
-
   // Compute augmented space
   if (mVerbose) {
     std::cout << "Computing augmented space with delay = " << mAugmentationDelay
               << ", number of samples is " << mInput.size()-mAugmentationDelay 
               << std::endl;
   }
-  ComputeAugmentedSpace(mInput, mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
+  mInput.ComputeAugmentedSpace(mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
 
   // Output augmented state
   if (mOutputAugmentedFilenameIsSet) {
     std::ofstream os;
     openFileForWriting(os, mOutputAugmentedFilename);
-    for(unsigned int i=0; i<(unsigned int)mWindowLength; i++) {
+    for(unsigned int i=0; i<(unsigned int)mAugmentedInputX.size(); i++) {
       os << mAugmentedInputX[i] << " " << mAugmentedInputY[i] << std::endl;
     }
     os.close();
@@ -102,6 +103,10 @@ void clitk::SignalMeanPositionFilter::Update() {
   // 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) {
@@ -116,7 +121,6 @@ void clitk::SignalMeanPositionFilter::Update() {
   // Fitting method
   if (mIsAdaptiveMethod) {
     AdaptiveFitEllipse(An);
-    //FitEllipse(An);   
   }
   else FitEllipse(An);
 
@@ -130,25 +134,33 @@ void clitk::SignalMeanPositionFilter::Update() {
   if (mOutputFilenameIsSet) {
     std::ofstream os;
     openFileForWriting(os, mOutputFilename);
-    os << AnInitial.ComputeCenter()[0] << " " << AnInitial.ComputeCenter()[1] << " " 
-       << AnInitial.ComputeSemiAxeLengths()[0] << " " << AnInitial.ComputeSemiAxeLengths()[1] << " " 
-       << AnInitial.ComputeAngleInRad();
-    for(int i=0; i<6; i++) os << AnInitial[i] << " " ;
-    os << std::endl;
     
-    os << An.ComputeCenter()[0] << " " << An.ComputeCenter()[1] << " " 
-       << An.ComputeSemiAxeLengths()[0] << " " << An.ComputeSemiAxeLengths()[1] << " " 
-       << An.ComputeAngleInRad();
-    for(int i=0; i<6; i++) os << An[i] << " ";
-    os << std::endl;
-    os.close();
+    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] << " ";
+        }
+      */
 
-    openFileForWriting(os, "centers.dat");
-    for(unsigned int i=0; i<mCenters.size(); i++) {
-      os << mCenters[i][0] << " " << mCenters[i][1] << std::endl;
+      os << " _Ellipse_ ";
+      for(int j=0; j<6; j++) os << An[j] << " ";
+      os << std::endl;
     }
     os.close();
-
   }
   
   if (mOutputResidualFilenameIsSet) {
@@ -164,27 +176,6 @@ void clitk::SignalMeanPositionFilter::Update() {
 //---------------------------------------------------------------------
 
 
-//---------------------------------------------------------------------
-void clitk::SignalMeanPositionFilter::ComputeAugmentedSpace(const clitk::Signal & input, 
-                                                            clitk::Signal & outputX, 
-                                                            clitk::Signal & outputY, 
-                                                            unsigned int delay) {
-  if (input.size() <= delay) {
-    std::cerr << "Error in signal length is " << input.size()
-              << " while delay is " << delay << " : too short. Abort." << std::endl;
-    exit(0);
-  }
-  outputX.resize(input.size()-delay);
-  outputY.resize(input.size()-delay);
-  for(unsigned int i=0; i<outputX.size(); i++) {
-    outputX[i] = input[i+delay];
-    outputY[i] = input[i];
-  }
-}
-
-//---------------------------------------------------------------------
-
-
 //---------------------------------------------------------------------
 void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
   int n = 0;
@@ -209,11 +200,162 @@ void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
 
 //---------------------------------------------------------------------
 void clitk::SignalMeanPositionFilter::AdaptiveFitEllipse(clitk::Ellipse & An) {  
-  for(unsigned int t=0; t<(unsigned int)args_info.t_arg; t++) {
-    DD(t);
+  // 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);
-    mCenters.push_back(An.ComputeCenter());
     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]);
+  }
+
+}
+//---------------------------------------------------------------------
+
diff --git a/tools/clitkSignalMeanPositionFilter.h b/tools/clitkSignalMeanPositionFilter.h
index 4cbc4f0..980c0a5 100644
--- a/tools/clitkSignalMeanPositionFilter.h
+++ b/tools/clitkSignalMeanPositionFilter.h
@@ -1,7 +1,10 @@
 /*=========================================================================
                                                                                 
   Program:   clitk
+  Module:    $RCSfile: clitkSignalMeanPositionFilter.h,v $
   Language:  C++
+  Date:      $Date: 2010/02/10 14:55:00 $
+  Version:   $Revision: 1.4 $
                                                                                 
   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   l'Image). All rights reserved. See Doc/License.txt or
@@ -20,6 +23,7 @@
 #include "clitkSignal.h"
 #include "clitkEllipse.h"
 #include "itkVector.h"
+#include <math.h>
 
 namespace clitk {
 
@@ -51,15 +55,28 @@ namespace clitk {
     bool mIsAdaptiveMethod;
     std::vector<double> mCurrentResidual;
     int mWindowLength;
-    std::vector<Vector2d> mCenters;
+    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 ComputeAugmentedSpace(const clitk::Signal & input, 
-                               clitk::Signal & outputX, 
-                               clitk::Signal & outputY, 
-                               unsigned int delay);
+    void ComputeIsoPhase(std::vector<clitk::Ellipse*> & l, 
+                         std::vector<double> & phase,
+                         std::vector<int> & cycles);
   };
   //---------------------------------------------------------------------
   
diff --git a/tools/clitkSignalMeanPositionTracking.cxx b/tools/clitkSignalMeanPositionTracking.cxx
index e63e929..939d7b9 100644
--- a/tools/clitkSignalMeanPositionTracking.cxx
+++ b/tools/clitkSignalMeanPositionTracking.cxx
@@ -1,7 +1,10 @@
 /*=========================================================================
                                                                                 
   Program:   clitk
+  Module:    $RCSfile: clitkSignalMeanPositionTracking.cxx,v $
   Language:  C++
+  Date:      $Date: 2010/02/10 14:55:00 $
+  Version:   $Revision: 1.4 $
                                                                                 
   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   l'Image). All rights reserved. See Doc/License.txt or
diff --git a/tools/clitkSignalMeanPositionTracking.ggo b/tools/clitkSignalMeanPositionTracking.ggo
index 46baa88..a1a7939 100644
--- a/tools/clitkSignalMeanPositionTracking.ggo
+++ b/tools/clitkSignalMeanPositionTracking.ggo
@@ -9,7 +9,12 @@ 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
+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