tools/clitkSignalMeanPositionFilter.cxx

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