tools/clitkSignalMeanPositionFilter.cxx

   1 /*=========================================================================
   2
   3   Program:   clitk
   4   Language:  C++
   5
   6   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   7   l'Image). All rights reserved. See Doc/License.txt or
   8   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
   9
  10      This software is distributed WITHOUT ANY WARRANTY; without even
  11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12      PURPOSE.  See the above copyright notices for more information.
  13
  14 =========================================================================*/
  15
  16 #include "clitkSignalMeanPositionFilter.h"
  17
  18 //---------------------------------------------------------------------
  19 void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPositionTracking & args) {
  20   args_info = args;
  21   mEtaIsSet = false;
  22   mOutputFilenameIsSet = false;
  23   mOutputResidualFilenameIsSet = false;
  24   mInput.Read(args_info.input_arg);
  25   mAugmentationDelay = args_info.delay_arg;
  26   mIsAdaptiveMethod = false;
  27   mWindowLength = -1;
  28   if (args_info.eta_given) {
  29     mEta = args_info.eta_arg;
  30     mEtaIsSet = true;
  31   }
  32   mMaxIteration = args_info.iter_arg;
  33   if (args_info.output_given) {
  34     mOutputFilenameIsSet = true;
  35     mOutputFilename = args_info.output_arg;
  36   }
  37   if (args_info.residual_given) {
  38     mOutputResidualFilenameIsSet = true;
  39     mOutputResidualFilename = args_info.residual_arg;
  40   }
  41   if (args_info.augmented_given) {
  42     mOutputAugmentedFilenameIsSet = true;
  43     mOutputAugmentedFilename = args_info.augmented_arg;
  44   }
  45   mVerbose = args_info.verbose_flag;
  46   mVerboseIteration = args_info.verbose_iteration_flag;
  47   if (args_info.L_given) {
  48     mIsAdaptiveMethod = true;
  49     mWindowLength = args_info.L_arg;
  50   }
  51 }
  52 //---------------------------------------------------------------------
  53
  54
  55 //---------------------------------------------------------------------
  56 void clitk::SignalMeanPositionFilter::Update() {
  57
  58   // DEBUG
  59 //   int e = 5;
  60 //   clitk::Signal temp;
  61 //   temp.resize(mInput.size()*e);
  62 //   for(unsigned int i=0; i<temp.size(); i++) {
  63 //     int index = lrint(floor((double)i/e));
  64 //     double res = (double)i/e - index;
  65 //     //    DD(index);
  66 //     //DD(res);
  67 //     temp[i] = mInput[index] + res*(mInput[index+1]-mInput[index])/e;
  68 //   }
  69 //   mInput.resize(temp.size());
  70   for  (unsigned int i=0; i<mInput.size(); i++) mInput[i] = mInput[i]+(double)i/(double)mInput.size();
  71   srand ( time(0) );
  72   for  (unsigned int i=0; i<mInput.size(); i++) {
  73     mInput[i] = (0.8+((double)rand()/ (double)RAND_MAX)*0.4) * mInput[i];
  74   }
  75
  76    {
  77      std::ofstream os;
  78      openFileForWriting(os, "input.dat");
  79      for(unsigned int i=0; i<mInput.size(); i++) os << mInput[i] << std::endl;
  80      os.close();
  81    }
  82
  83
  84   // Compute augmented space
  85   if (mVerbose) {
  86     std::cout << "Computing augmented space with delay = " << mAugmentationDelay
  87               << ", number of samples is " << mInput.size()-mAugmentationDelay
  88               << std::endl;
  89   }
  90   ComputeAugmentedSpace(mInput, mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
  91
  92   // Output augmented state
  93   if (mOutputAugmentedFilenameIsSet) {
  94     std::ofstream os;
  95     openFileForWriting(os, mOutputAugmentedFilename);
  96     for(unsigned int i=0; i<(unsigned int)mWindowLength; i++) {
  97       os << mAugmentedInputX[i] << " " << mAugmentedInputY[i] << std::endl;
  98     }
  99     os.close();
 100   }
 101
 102   // Initial starting ellipse
 103   Ellipse An; // starting point
 104   if (!mEtaIsSet) mEta = -1;
 105   An.InitialiseEllipseFitting(mEta, mWindowLength, mAugmentedInputX, mAugmentedInputY);
 106   mEta = An.GetEta();
 107   if (mVerbose) {
 108     std::cout << "Eta is " << mEta << std::endl;
 109   }
 110   Ellipse AnInitial(An);
 111   if (mVerbose) {
 112     std::cout << "Set initial ellipse to c=" << An.ComputeCenter()
 113               << " and axes=" << An.ComputeSemiAxeLengths() << std::endl;
 114   }
 115
 116   // Fitting method
 117   if (mIsAdaptiveMethod) {
 118     AdaptiveFitEllipse(An);
 119     //FitEllipse(An);
 120   }
 121   else FitEllipse(An);
 122
 123   // Output
 124   if (mVerbose) {
 125     std::cout << "Result    is "  << An << std::endl;
 126     std::cout << "Center    is "  << An.ComputeCenter() << std::endl;
 127     std::cout << "SemiAxis are " << An.ComputeSemiAxeLengths() << std::endl;
 128     std::cout << "Angle     is " << rad2deg(An.ComputeAngleInRad()) << " deg" << std::endl;
 129   }
 130   if (mOutputFilenameIsSet) {
 131     std::ofstream os;
 132     openFileForWriting(os, mOutputFilename);
 133     os << AnInitial.ComputeCenter()[0] << " " << AnInitial.ComputeCenter()[1] << " "
 134        << AnInitial.ComputeSemiAxeLengths()[0] << " " << AnInitial.ComputeSemiAxeLengths()[1] << " "
 135        << AnInitial.ComputeAngleInRad();
 136     for(int i=0; i<6; i++) os << AnInitial[i] << " " ;
 137     os << std::endl;
 138
 139     os << An.ComputeCenter()[0] << " " << An.ComputeCenter()[1] << " "
 140        << An.ComputeSemiAxeLengths()[0] << " " << An.ComputeSemiAxeLengths()[1] << " "
 141        << An.ComputeAngleInRad();
 142     for(int i=0; i<6; i++) os << An[i] << " ";
 143     os << std::endl;
 144     os.close();
 145
 146     openFileForWriting(os, "centers.dat");
 147     for(unsigned int i=0; i<mCenters.size(); i++) {
 148       os << mCenters[i][0] << " " << mCenters[i][1] << std::endl;
 149     }
 150     os.close();
 151
 152   }
 153
 154   if (mOutputResidualFilenameIsSet) {
 155     std::ofstream os;
 156     openFileForWriting(os, mOutputResidualFilename);
 157     for(unsigned int i=0; i<mCurrentResidual.size(); i++) {
 158       os.precision(10);
 159       os << mCurrentResidual[i] << std::endl;
 160     }
 161   os.close();
 162   }
 163 }
 164 //---------------------------------------------------------------------
 165
 166
 167 //---------------------------------------------------------------------
 168 void clitk::SignalMeanPositionFilter::ComputeAugmentedSpace(const clitk::Signal & input,
 169                                                             clitk::Signal & outputX,
 170                                                             clitk::Signal & outputY,
 171                                                             unsigned int delay) {
 172   if (input.size() <= delay) {
 173     std::cerr << "Error in signal length is " << input.size()
 174               << " while delay is " << delay << " : too short. Abort." << std::endl;
 175     exit(0);
 176   }
 177   outputX.resize(input.size()-delay);
 178   outputY.resize(input.size()-delay);
 179   for(unsigned int i=0; i<outputX.size(); i++) {
 180     outputX[i] = input[i+delay];
 181     outputY[i] = input[i];
 182   }
 183 }
 184
 185 //---------------------------------------------------------------------
 186
 187
 188 //---------------------------------------------------------------------
 189 void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
 190   int n = 0;
 191   mCurrentResidual.clear();
 192   while (n<mMaxIteration) {
 193     double r = An.EllipseFittingNextIteration();
 194     mCurrentResidual.push_back(r);
 195     n++;
 196     if (mVerboseIteration) {
 197       std::cout.precision(3);
 198       std::cout << "A(" << n << ")=" << An
 199                 << " c="  << An.ComputeCenter()
 200                 << " ax=" << An.ComputeSemiAxeLengths()
 201                 << " theta=" << rad2deg(An.ComputeAngleInRad());
 202       std::cout.precision(10);
 203       std::cout << " R=" << r << std::endl;
 204     }
 205   }
 206 }
 207 //---------------------------------------------------------------------
 208
 209
 210 //---------------------------------------------------------------------
 211 void clitk::SignalMeanPositionFilter::AdaptiveFitEllipse(clitk::Ellipse & An) {
 212   for(unsigned int t=0; t<(unsigned int)args_info.t_arg; t++) {
 213     DD(t);
 214     An.UpdateSMatrix(t, mWindowLength, mAugmentedInputX, mAugmentedInputY);
 215     mCenters.push_back(An.ComputeCenter());
 216     FitEllipse(An);
 217   }
 218 }
 219 //---------------------------------------------------------------------