X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkSignalMeanPositionFilter.cxx;h=912d9ee63e052f50de8c1c76c6411800b082581b;hb=96a65c0dafae2081d6da4975bddf7567ab6c05d0;hp=905af5ece48b8886b93db2e98db76d773dbe55b8;hpb=253e4eea321f0a154e1c57834f36d4e4c4b34ef1;p=clitk.git

diff --git a/tools/clitkSignalMeanPositionFilter.cxx b/tools/clitkSignalMeanPositionFilter.cxx
index 905af5e..912d9ee 100644
--- a/tools/clitkSignalMeanPositionFilter.cxx
+++ b/tools/clitkSignalMeanPositionFilter.cxx
@@ -1,10 +1,7 @@
 /*=========================================================================
                                                                                 
   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
@@ -28,8 +25,6 @@ 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;
@@ -53,28 +48,6 @@ 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);
-  }
-
 }
 //---------------------------------------------------------------------
 
@@ -82,19 +55,45 @@ 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;
   }
-  mInput.ComputeAugmentedSpace(mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
+  ComputeAugmentedSpace(mInput, 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++) {
+    for(unsigned int i=0; i<(unsigned int)mWindowLength; i++) {
       os << mAugmentedInputX[i] << " " << mAugmentedInputY[i] << std::endl;
     }
     os.close();
@@ -103,10 +102,6 @@ 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) {
@@ -121,6 +116,7 @@ void clitk::SignalMeanPositionFilter::Update() {
   // Fitting method
   if (mIsAdaptiveMethod) {
     AdaptiveFitEllipse(An);
+    //FitEllipse(An);   
   }
   else FitEllipse(An);
 
@@ -134,33 +130,25 @@ 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;
     
-    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 << 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();
 
-      os << " _Ellipse_ ";
-      for(int j=0; j<6; j++) os << An[j] << " ";
-      os << std::endl;
+    openFileForWriting(os, "centers.dat");
+    for(unsigned int i=0; i<mCenters.size(); i++) {
+      os << mCenters[i][0] << " " << mCenters[i][1] << std::endl;
     }
     os.close();
+
   }
   
   if (mOutputResidualFilenameIsSet) {
@@ -176,6 +164,27 @@ 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;
@@ -200,162 +209,11 @@ void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
 
 //---------------------------------------------------------------------
 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++) {
+  for(unsigned int t=0; t<(unsigned int)args_info.t_arg; t++) {
+    DD(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]);
-  }
-
-}
-//---------------------------------------------------------------------
-