common/clitkSignal.cxx

   1 #ifndef CLITKSIGNAL_CXX
   2 #define CLITKSIGNAL_CXX
   3
   4 #include "clitkSignal.h"
   5 #include <fstream>
   6
   7 namespace clitk {
   8
   9   //---------------------------------------------------------------------
  10   Signal::Signal(const Signal & s) { // Copy
  11     CopyFrom(s);
  12   }
  13   //---------------------------------------------------------------------
  14
  15
  16   //---------------------------------------------------------------------
  17   void Signal::CopyFrom(const Signal & s) {
  18     SetSamplingPeriod(s.GetSamplingPeriod());
  19     // DD(GetSamplingPeriod());
  20     resize(s.size());
  21     for(uint i=0; i<size(); i++) {
  22       m_Data[i] = s[i];
  23     }
  24   }
  25   //---------------------------------------------------------------------
  26
  27
  28
  29   //---------------------------------------------------------------------
  30   void Signal::Read(string fileName){
  31     std::ifstream signalStream(fileName.c_str());
  32     SignalValueType point;
  33     if(!signalStream.is_open()){
  34       std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
  35       return;
  36     }
  37     skipComment(signalStream);
  38     while(!signalStream.eof()) {
  39       skipComment(signalStream);
  40       signalStream >> point;
  41       skipComment(signalStream);
  42       m_Data.push_back(point);
  43       //      cout << point << endl;
  44     }
  45     signalStream.close();
  46     SetSamplingPeriod(1.);
  47   }
  48   //---------------------------------------------------------------------
  49
  50
  51   //---------------------------------------------------------------------
  52   bool Signal::Read(string fileName, int col){
  53     ifstream signalStream(fileName.c_str());
  54     SignalValueType point;
  55     if(!signalStream.is_open()){
  56       std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
  57       return false;
  58     }
  59     skipComment(signalStream);
  60     while(!signalStream.eof()) {
  61       skipComment(signalStream);
  62
  63       // Read one line
  64       std::string line;
  65       std::getline(signalStream, line);
  66
  67       // Get column nb col
  68       istringstream iss(line);
  69       for(int i=0; i<col; i++) { // skip col-1 first column
  70         string sub;
  71         iss >> sub;
  72         // DD(sub);
  73         // DD(sub.size());
  74         //         DD(iss.bad());
  75         if (!iss) {
  76           std::cerr << "ERROR: no col n" << col << " in the line '" << line << "' ?" << std::endl;
  77           return false; //exit(0);
  78         }
  79       }
  80       iss >> point;
  81       skipComment(signalStream);
  82       m_Data.push_back(point);
  83       //      cout << point << endl;
  84     }
  85     signalStream.close();
  86     SetSamplingPeriod(1.);
  87     return true;
  88   }
  89   //---------------------------------------------------------------------
  90
  91
  92
  93   //---------------------------------------------------------------------
  94   //Convert 1D image to signal
  95   Signal Signal::ConvertImageToSignal( Signal::ImageType::Pointer image)
  96   {
  97     //empty signal
  98     Signal signal;
  99
 100     //make an image iterator
 101     itk::ImageRegionConstIterator<ImageType> it(image,image->GetLargestPossibleRegion());
 102     it.Begin();
 103
 104     //copy
 105     while(!it.IsAtEnd())
 106       {
 107         signal.push_back(it.Get());
 108         ++it;
 109       }
 110
 111     //Spacing
 112     signal.SetSamplingPeriod(image->GetSpacing()[0]);
 113
 114     return signal;
 115   }
 116   //---------------------------------------------------------------------
 117
 118
 119   //---------------------------------------------------------------------
 120   //Convert 1D signal to image
 121   Signal::ImageType::Pointer Signal::ConvertSignalToImage(Signal signal)
 122   {
 123     //empty image
 124     ImageType::Pointer image =ImageType::New();
 125     ImageType::RegionType region;
 126     ImageType::RegionType::IndexType index;
 127     index[0]=0;
 128     ImageType::RegionType::SizeType size;
 129     size[0]=signal.size();
 130
 131     region.SetIndex(index);
 132     region.SetSize(size);
 133
 134     image->SetRegions(region);
 135     image->Allocate();
 136
 137     //make an image iterator
 138     itk::ImageRegionIterator<ImageType> mIt(image,image->GetLargestPossibleRegion());
 139     mIt.Begin();
 140
 141     //make a signal iterator
 142     Signal::const_iterator sIt=signal.begin();
 143
 144     //copy
 145     while(sIt!=signal.end())
 146       {
 147         mIt.Set(*sIt);
 148         sIt++;++mIt;
 149       }
 150
 151
 152     //spacing
 153     ImageType::SpacingType spacing;
 154     spacing[0]=signal.GetSamplingPeriod();
 155     image->SetSpacing(spacing);
 156
 157     return image;
 158
 159   }
 160   //---------------------------------------------------------------------
 161
 162
 163
 164   //---------------------------------------------------------------------
 165   void Signal::Write(const string fileName){
 166     ofstream signalStream(fileName.c_str());
 167     if(!signalStream.is_open()){
 168       cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Write" << endl;
 169       return;
 170     }
 171
 172     iterator it=begin();
 173     while(it!=end()) {
 174       signalStream << *it << endl;
 175       it++;
 176     }
 177     signalStream.close();
 178   }
 179   //---------------------------------------------------------------------
 180
 181
 182   //---------------------------------------------------------------------
 183   Signal & Signal::operator/=(Signal & div){
 184     if(size()!= div.size())
 185       {
 186         std::cerr << "Error: signal size must be the same!" << std::endl;
 187         return  *this;
 188       }
 189     iterator it=begin();
 190     iterator itD;
 191     itD=div.begin();
 192     while(it!=end())
 193       {
 194         if(*itD!=0)
 195           *it/=*itD;
 196         else
 197           cerr << "Division by 0 in operator/= skipped" << endl;
 198         it++;itD++;
 199       }
 200     return *this;
 201   }
 202   //---------------------------------------------------------------------
 203
 204
 205   //---------------------------------------------------------------------
 206   Signal & Signal::operator*=(Signal & mul){
 207     if(size()!= mul.size())
 208       {
 209         std::cerr << "Error: signal size must be the same!" << std::endl;
 210         return  *this;
 211       }
 212     iterator it=begin();
 213     iterator itD;
 214     itD=mul.begin();
 215     while(it!=end()){
 216       *it *= *itD;
 217       it++;itD++;
 218     }
 219     return *this;
 220   }
 221   //---------------------------------------------------------------------
 222
 223
 224   //---------------------------------------------------------------------
 225   Signal Signal::Normalize(double newMin,double newMax){
 226     Signal temp (m_SamplingPeriod);
 227     vector<double> extrema=GetGlobalMinMax();
 228     iterator itSig=begin();
 229     while(itSig!=end()){
 230       temp.push_back( ((*itSig)-extrema[0])*(newMax-newMin)/(extrema[1]-extrema[0]) + newMin );
 231       itSig++;
 232     }
 233     return temp;
 234   }
 235   //---------------------------------------------------------------------
 236
 237
 238   //---------------------------------------------------------------------
 239   vector<double> Signal::GetGlobalMinMax() const {
 240     vector<double> extrema(2);
 241     if(size()==0){
 242       cerr << "ERROR: GetExtrema / No signal" << endl;
 243       return extrema;
 244     }
 245     extrema[0]=m_Data[0];
 246     extrema[1]=m_Data[0];
 247     for(unsigned int i=1;i<m_Data.size();i++){
 248       if(extrema[0]>m_Data[i]) extrema[0]=m_Data[i];
 249       if(extrema[1]<m_Data[i]) extrema[1]=m_Data[i];
 250     }
 251     return extrema;
 252   }
 253   //---------------------------------------------------------------------
 254
 255
 256   //---------------------------------------------------------------------
 257   double Signal::GetGlobalMean() const {
 258     double m = 0.0;
 259     for(unsigned int i=0;i<m_Data.size();i++){
 260       m += m_Data[i];
 261     }
 262     return m/(double)m_Data.size();
 263   }
 264   //---------------------------------------------------------------------
 265
 266
 267   //---------------------------------------------------------------------
 268   Signal Signal::MovingAverageFilter ( unsigned int length)  {
 269
 270     Signal temp(m_SamplingPeriod);
 271
 272     for (unsigned int i=0; i <size(); i++)
 273       {
 274         double accumulator=0.;
 275         unsigned int scale=0;
 276         for (unsigned int j= max(0,static_cast<int>(i)-static_cast<int>(length));  j<=min(size(), i+length); j++)
 277           {
 278             accumulator+=m_Data[j];
 279             scale++;
 280           }
 281         temp.push_back(accumulator/scale);
 282       }
 283     return temp;
 284   }
 285   //---------------------------------------------------------------------
 286
 287
 288   //---------------------------------------------------------------------
 289   Signal Signal::GaussLikeFilter ( )  {
 290
 291     Signal temp(m_SamplingPeriod);
 292     if (size()<2)
 293       return *this;
 294     else
 295       {
 296         //first sample: mirrorring BC
 297         temp.push_back((2.*m_Data[0]+2*m_Data[1])/4.);
 298
 299         //middle samples
 300         for (unsigned int i=1; i <size()-1; i++)
 301           temp.push_back((2.*m_Data[i]+m_Data[i-1]+m_Data[i+1])/4.);
 302
 303         //end sample: mirrorring BC
 304         temp.push_back((2.*m_Data[size()-1]+2*m_Data[size()-2])/4.);
 305         return temp;
 306       }
 307   }
 308   //---------------------------------------------------------------------
 309
 310
 311   //---------------------------------------------------------------------
 312   Signal Signal::NormalizeMeanStdDev(double newMean,double newStdDev)
 313   {
 314     iterator itSig=begin();
 315     double sum=0, sum2=0;
 316     Signal temp(m_SamplingPeriod);
 317
 318     //Calculate the mean and the std dev
 319     while(itSig!=end())
 320       {
 321         sum += *itSig;
 322         sum2 += (*itSig) * (*itSig);
 323         itSig++;
 324       }
 325     double oldMean=sum/size();
 326     double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
 327
 328     //find a and b
 329     double a = newStdDev/oldStdDev;
 330     double b = newMean - a * oldMean;
 331     itSig=begin();
 332     while(itSig!=end())
 333       {
 334         temp.push_back( a *(*itSig) + b );
 335         itSig++;
 336       }
 337     return temp;
 338   }
 339   //---------------------------------------------------------------------
 340
 341
 342   //---------------------------------------------------------------------
 343 //   Signal Signal::HighPassFilter (double sampPeriod, double cutOffFrequency )
 344 //   {
 345 //     //output
 346 //     Signal temp(m_SamplingPeriod);
 347 //     temp.resize(size());
 348 //
 349 //     //the fft
 350 //     SIGNAL_FFT_TYPE fft;
 351 //
 352 //     //calculate the cut off frequency
 353 //     unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
 354 //
 355 //     //forward fft with empty fft
 356 //     if(fft.size()==0)  OneDForwardFourier(*this, fft);
 357 //
 358 //     //remove the frequencies
 359 //     for(unsigned int i=0;i<samp && i<fft.size();i++)
 360 //       fft[i]=complex<double>(0.,0.);
 361 //
 362 //     //backward with remaining frequencies
 363 //     OneDBackwardFourier(fft,temp);
 364 //     return temp;
 365 //   }
 366   //---------------------------------------------------------------------
 367
 368
 369   //---------------------------------------------------------------------
 370 //   Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency )
 371 //   {
 372 //     //output
 373 //     Signal temp(m_SamplingPeriod);
 374 //     temp.resize(size());
 375 //
 376 //     //the fft
 377 //     SIGNAL_FFT_TYPE fft;
 378 //
 379 //     //calculate the cut off frequency
 380 //     unsigned int samp=lrint(cutOffFrequency*static_cast<double>(size())*sampPeriod);
 381 //
 382 //     //forward fft with empty fft
 383 //     if(fft.size()==0)  OneDForwardFourier(*this, fft);
 384 //     unsigned int fsize=fft.size();
 385 //
 386 //     //remove the frequencies
 387 //     unsigned int limit=min (samp, fsize);
 388 //     for(unsigned int i=limit;i<fft.size();i++)
 389 //      fft[i]=complex<double>(0.,0.);
 390 //
 391 //     //backward with remaining frequencies
 392 //     OneDBackwardFourier(fft,temp);
 393 //     return temp;
 394 //   }
 395   //---------------------------------------------------------------------
 396
 397
 398   //---------------------------------------------------------------------
 399 //   void  Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft)
 400 //   {
 401 //     //Create output array
 402 //     fft.resize(input.size()/2+1);
 403 //     //Temp copy
 404 //     double *tempCopy=new double[size()];
 405 //     copy(begin(), end(), tempCopy);
 406 //
 407 //     //Forward Fourier Transform
 408 //     fftw_plan p;
 409 //     p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast<fftw_complex*>(&(fft[0])),FFTW_ESTIMATE);
 410 //     fftw_execute(p);
 411 //     fftw_destroy_plan(p);
 412 //     //delete tempCopy;
 413 //     return;
 414 //   }
 415   //---------------------------------------------------------------------
 416
 417
 418   //---------------------------------------------------------------------
 419 //   void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output)
 420 //   {
 421 //
 422 //     //Backward
 423 //     fftw_plan p;
 424 //     p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast<fftw_complex*>(&(fft[0])),&(output[0]),FFTW_ESTIMATE);
 425 //     fftw_execute(p);
 426 //     fftw_destroy_plan(p);
 427 //
 428 //     vector<double>::iterator it=output.begin();
 429 //     while(it!=output.end()){
 430 //       *it /= (double)output.size();
 431 //       it++;
 432 //     }
 433 //     return;
 434 //   }
 435   //---------------------------------------------------------------------
 436
 437
 438   //---------------------------------------------------------------------
 439 //   double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft)
 440 //   {
 441 //
 442 //     if(fft.size()==0) OneDForwardFourier(sig,fft);
 443 //     int posMax=1;
 444 //     double amplitude, amplitudeMax=abs(fft[1]);
 445 //     for(unsigned int i=1;i<fft.size();i++){
 446 //       amplitude=abs(fft[i]);
 447 //       if(amplitude>amplitudeMax){
 448 //      posMax=i;
 449 //      amplitudeMax=amplitude;
 450 //       }
 451 //     }
 452 //     return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod()));
 453 //   }
 454   //---------------------------------------------------------------------
 455
 456
 457   //---------------------------------------------------------------------
 458   Signal Signal::DetectLocalExtrema(unsigned int width)
 459   {
 460     Signal temp(m_SamplingPeriod);
 461     bool isMin, isMax;
 462     unsigned int upper, lower;
 463
 464     //has to be at least 1
 465     width=max(static_cast<int>(width), 1);
 466
 467     for(unsigned int i=0 ; i < size() ; i++){
 468       isMax=true;
 469       isMin=true;
 470
 471       for(unsigned int j=1; j< width+1; j++)
 472         {
 473           //set the boundaries
 474           upper = min( size(), i+j);
 475           lower = max( static_cast<int>(0), (int)i-(int)j);
 476
 477           //check if max
 478           if( ! (m_Data[i] >= m_Data[lower] && m_Data[i] >= m_Data[upper]))isMax=false;
 479
 480           //check if min
 481           if( ! (m_Data[i]<= m_Data[lower] && m_Data[i] <= m_Data[upper])) isMin=false;
 482
 483           //if neither, go to the next value
 484           if( (!isMax) && (!isMin)) break;
 485         }
 486
 487       if (isMax) temp.push_back(1);
 488       else if (isMin) temp.push_back(0);
 489       else temp.push_back(0.5);
 490     }
 491     return temp;
 492   }
 493   //---------------------------------------------------------------------
 494
 495
 496   //---------------------------------------------------------------------
 497   Signal Signal::LimPhase()
 498   {
 499
 500     //phase is defined as going from 0 to 1 linearly between end expiration
 501     Signal phase(m_SamplingPeriod);
 502     phase.resize(size());
 503
 504     unsigned int firstBeginCycle=0;
 505     unsigned int firstEndCycle=0;
 506     unsigned int beginCycle=0;
 507
 508     //=========================================
 509     //Search the first value in extrema not 0.5
 510     while(m_Data[beginCycle]!=1)
 511       {
 512         beginCycle++;
 513       }
 514
 515     //We search the corresponding end
 516     unsigned int endCycle=beginCycle+1;
 517     while(endCycle != size() && m_Data[endCycle]!=1){
 518       endCycle++;
 519
 520     }
 521
 522     //============================================================
 523     //Calculate phase at the beginning (before the first extremum)
 524     for(unsigned int i=0 ; i<beginCycle ; i++)
 525       {
 526
 527         //if before first cycle is longer than first cycle: let it go from 0 to 1
 528         if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
 529           {
 530             phase[i] = (double)(i-0)/(double)(beginCycle-0);
 531           }
 532         //copy the phase values later
 533         else
 534           {
 535             firstBeginCycle=beginCycle;
 536             firstEndCycle=endCycle;
 537           }
 538
 539       }
 540
 541     //===================================================================
 542     //Middle part
 543     while(endCycle != size()){
 544
 545       //fill between extrema
 546       for(unsigned int i=beginCycle ; i<endCycle ; i++){
 547         phase[i] = (double)(i-beginCycle)/(double)(endCycle-beginCycle);
 548       }
 549
 550       //Find next cycle
 551       beginCycle = endCycle++;
 552       while(endCycle != size() && m_Data[endCycle]!=1)
 553         endCycle++;
 554     }
 555
 556     //===================================================================
 557     //Calculate phase at the end (after the last extremum)
 558     endCycle = beginCycle--;
 559
 560     //count backwards till the previous
 561     while(m_Data[beginCycle]!=1) beginCycle--;
 562     for(unsigned int i=endCycle ; i<size() ; i++)
 563       {
 564
 565         //after last extrema is longer then last cycle
 566         if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
 567
 568           //make the last part go till 1
 569           phase[i] = (double)(i-endCycle)/(double)(size()-endCycle);
 570         }
 571
 572         //the last part is shorter, copy the last cycle values
 573         else{
 574           phase[i] = phase[i -(endCycle-beginCycle)];
 575         }
 576       }
 577
 578     //===================================================================
 579     //check it some remains to be copied in the beginning
 580     if (firstBeginCycle!=0)
 581       {
 582         for(unsigned int i=0 ; i<firstBeginCycle ; i++)
 583           phase[i] =phase[i+firstEndCycle-firstBeginCycle];
 584       }
 585
 586     return phase;
 587   }
 588   //---------------------------------------------------------------------
 589
 590
 591
 592   //---------------------------------------------------------------------
 593   Signal Signal::MonPhase()
 594   {//monPhase
 595
 596     //phase is defined as going from 0 to 1 linearly between end expiration
 597     Signal phase(m_SamplingPeriod);
 598     phase.resize(size());
 599
 600     unsigned int firstBeginCycle=0;
 601     unsigned int firstEndCycle=0;
 602     unsigned int cycleCounter=0;
 603     unsigned int beginCycle=0;
 604
 605     //===================================================================
 606     //Search the first value in extrema not 0.5
 607     while(m_Data[beginCycle]!=1)
 608       {
 609         beginCycle++;
 610
 611       }
 612     //We search the corresponding end
 613     unsigned int endCycle=beginCycle+1;
 614     while(endCycle != size() && m_Data[endCycle]!=1){
 615       endCycle++;
 616
 617     }
 618
 619     //===================================================================
 620     //Calculate phase at the beginning (before the first extremum)
 621     for(unsigned int i=0 ; i<beginCycle ; i++)
 622       {
 623
 624         //if before first cycle is longer than first cycle: let it go from 0 to 1
 625         if((double)beginCycle/(double)(endCycle-beginCycle) > 1)
 626           {
 627             phase[i] = (double)(i-0)/(double)(beginCycle-0);
 628           }
 629         //copy the phase values later
 630         else
 631           {
 632             firstBeginCycle=beginCycle;
 633             firstEndCycle=endCycle;
 634           }
 635
 636       }
 637
 638     //===================================================================
 639     //Middle part
 640     while(endCycle != size()){
 641
 642       cycleCounter++;
 643       //fill between extrema
 644       for(unsigned int i=beginCycle ; i<endCycle ; i++)
 645         phase[i] = (double)cycleCounter+(double)(i-beginCycle)/(double)(endCycle-beginCycle);
 646
 647
 648       //Find next cycle
 649       beginCycle = endCycle++;
 650       while(endCycle != size() && m_Data[endCycle]!=1)
 651         endCycle++;
 652     }
 653
 654     //===================================================================
 655     //Calculate phase at the end (after the last extremum)
 656     endCycle = beginCycle--;
 657
 658     //count backwards till the previous
 659     while(m_Data[beginCycle]!=1) beginCycle--;
 660     for(unsigned int i=endCycle ; i<size() ; i++)
 661       {
 662
 663         //after last extrema is longer then last cycle
 664         if((double)(size()-endCycle)/(double)(endCycle-beginCycle) > 1){
 665
 666           //make the last part go till 1
 667           phase[i] = (double)cycleCounter+(double)(i-endCycle)/(double)(size()-endCycle);
 668         }
 669
 670         //the last part is shorter, copy the last cycle values
 671         else{
 672           phase[i] = phase[i -(endCycle-beginCycle)]+1;
 673         }
 674       }
 675
 676     //===================================================================
 677     //check it some remains to be copied in the beginning
 678     if (firstBeginCycle!=0)
 679       {
 680         for(unsigned int i=0 ; i<firstBeginCycle ; i++)
 681           phase[i] =phase[i+firstEndCycle-firstBeginCycle]-1;
 682       }
 683
 684     return phase;
 685   }
 686   //---------------------------------------------------------------------
 687
 688
 689   //---------------------------------------------------------------------
 690   Signal Signal::MonPhaseDE(double eEPhaseValue, double eIPhaseValue)
 691   {
 692     //First calculate monPhase
 693     Signal monPhase=this->MonPhase();
 694
 695     //Create an empty signal
 696     Signal phase(size(), -1);
 697     phase.SetSamplingPeriod(m_SamplingPeriod);
 698
 699     //Fill in the values at the extrema position
 700     iterator phaseIt=phase.begin();
 701     iterator monPhaseIt=monPhase.begin();
 702     iterator extremaIt =begin();
 703     while (extremaIt!= end())
 704       {
 705         if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
 706         else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
 707         extremaIt++; phaseIt++;monPhaseIt++;
 708       }
 709
 710     return phase;
 711
 712   }
 713   //---------------------------------------------------------------------
 714
 715
 716   //---------------------------------------------------------------------
 717   Signal Signal::LinearlyInterpolateScatteredValues()
 718   {
 719     //Linearly interpolate the values in between
 720     unsigned int i=0;
 721     unsigned int beginCycle=0;
 722     unsigned int endCycle=0;
 723
 724     //Create a new signal
 725     Signal temp(size(),-1);
 726     temp.SetSamplingPeriod(m_SamplingPeriod);
 727
 728     //start from the first value
 729     while (m_Data[i]==-1)i++;
 730     beginCycle=i;
 731     i++;
 732
 733     //Go to the next
 734     while ( (m_Data[i]==-1) && (i<size()) )i++;
 735     while (i < size())
 736       {
 737         endCycle=i;
 738
 739         //fill in in between
 740         for (unsigned int k=beginCycle;k<endCycle+1; k++)
 741           temp[k]=m_Data[beginCycle]+(double)(k-beginCycle)*
 742             (m_Data[endCycle]-m_Data[beginCycle])/(double)(endCycle-beginCycle);
 743
 744         //swap and move on
 745         beginCycle=endCycle;
 746         i++;
 747         while( (i< size()) && (m_Data[i]==-1) ) i++;
 748       }
 749
 750     //For the last part
 751     if (beginCycle!= size()-1)
 752       {
 753         //continue with same slope
 754         for (unsigned int k=beginCycle+1; k<size();k++)
 755           temp[k]=temp[beginCycle]+(double)(k-beginCycle)*(temp[beginCycle]-temp[beginCycle-1]);
 756
 757       }
 758
 759     //For the first part
 760     if (temp[0]==-1)
 761       {
 762         i=0;
 763         while (temp[i]==-1)i++;
 764         beginCycle=i;
 765         i++;
 766         while (m_Data[i]==-1)i++;
 767         endCycle=i;
 768
 769         //if the first filled half cycle is longer, copy it
 770         if(beginCycle<(endCycle-beginCycle))
 771           {
 772             for (unsigned int k=0; k< beginCycle;k++)
 773               temp[k]=temp[k+endCycle-beginCycle];
 774           }
 775
 776         //if the first filled half cycle is longer, start from zero
 777         else
 778           {
 779             for (unsigned int k=0; k< beginCycle;k++)
 780               temp[k]=k*temp[beginCycle]/(beginCycle);
 781           }
 782       }
 783
 784     return temp;
 785   }
 786   //---------------------------------------------------------------------
 787
 788
 789
 790 //   //---------------------------------------------------------------------
 791 //   Signal Signal::ApproximateScatteredValuesWithBSplines(unsigned int splineOrder, unsigned int numberOfControlPoints)
 792 //   {
 793 //     //filter requires a vector pixelType
 794 //     typedef itk::PointSet<VectorType, 1> PointSetType;
 795 //     PointSetType::Pointer pointSet=PointSetType::New();
 796 //     typedef PointSetType::PointType PointType;
 797
 798 //     unsigned int i=0;
 799 //     unsigned int pointIndex=0;
 800 //     while (i< size())
 801 //       {
 802 //         if(m_Data[i]!=-1)
 803 //           {
 804 //             PointType p;
 805 //             p[0]= i;//JV spacing is allways 1
 806 //             pointSet->SetPoint( pointIndex, p );
 807 //             pointSet->SetPointData( pointIndex, m_Data[i] );
 808 //             pointIndex++;
 809 //           }
 810 //         i++;
 811 //       }
 812 // =======
 813 //  //---------------------------------------------------------------------
 814 //  Signal Signal::ApproximateScatteredValuesWithBSplines(unsigned int splineOrder, unsigned int numberOfControlPoints)
 815 //   {
 816 //     //filter requires a vector pixelType
 817 //     typedef itk::PointSet<VectorType, 1> PointSetType;
 818 //     PointSetType::Pointer pointSet=PointSetType::New();
 819 //     typedef PointSetType::PointType PointType;
 820
 821 //     unsigned int i=0;
 822 //     unsigned int pointIndex=0;
 823 //     while (i< size())
 824 //       {
 825 //       if(m_Data[i]!=-1)
 826 //      {
 827 //        PointType p;
 828 //        p[0]= i;//JV spacing is allways 1
 829 //        pointSet->SetPoint( pointIndex, p );
 830 //        pointSet->SetPointData( pointIndex, m_Data[i] );
 831 //        pointIndex++;
 832 //      }
 833 //       i++;
 834 //       }
 835
 836 //     //define the output signal properties
 837 //     ImageType::RegionType::SizeType outputSize;
 838 //     outputSize[0]= size();
 839 //     ImageType::PointType outputOrigin;
 840 //     outputOrigin[0]=0.0;//JV may need to be changed
 841 //     ImageType::SpacingType outputSpacing;
 842 //     outputSpacing[0]=1; //JV add relation to the original signal spacing
 843
 844 //     //Convert
 845 //     typedef itk::BSplineScatteredDataPointSetToImageFilter< PointSetType, VectorImageType > PointSetToImageFilterType;
 846 //     PointSetToImageFilterType::Pointer pointSetToImageFilter= PointSetToImageFilterType::New();
 847 //     pointSetToImageFilter->SetInput(pointSet);
 848 //     pointSetToImageFilter->SetSplineOrder(splineOrder);//JV
 849 //     pointSetToImageFilter->SetSize(outputSize);
 850 //     pointSetToImageFilter->SetOrigin(outputOrigin);
 851 //     pointSetToImageFilter->SetSpacing(outputSpacing);
 852
 853 //     //Convert to
 854 //     itk::FixedArray<unsigned int,1> num;
 855 //     num[0]=numberOfControlPoints;
 856 //     pointSetToImageFilter->SetNumberOfControlPoints(num);//JV
 857 //     pointSetToImageFilter->Update();
 858 //     VectorImageType::Pointer approximatedSignal=pointSetToImageFilter->GetOutput();
 859
 860 //     //Convert and return
 861 //     return ConvertVectorImageToSignal(approximatedSignal);
 862 //   }
 863 //   //---------------------------------------------------------------------
 864
 865
 866   //---------------------------------------------------------------------
 867   Signal Signal::ConvertVectorImageToSignal(VectorImageType::Pointer image)
 868   {
 869     //empty signal
 870     Signal signal;
 871
 872     //make an image iterator
 873     itk::ImageRegionConstIterator<VectorImageType> it(image,image->GetLargestPossibleRegion());
 874     it.Begin();
 875
 876     //copy
 877     while(!it.IsAtEnd())
 878       {
 879         signal.push_back(it.Get()[0]);
 880         ++it;
 881       }
 882
 883     //Spacing
 884     signal.SetSamplingPeriod(image->GetSpacing()[0]);
 885
 886     return signal;
 887   }
 888   //---------------------------------------------------------------------
 889
 890
 891   //---------------------------------------------------------------------
 892   Signal Signal::LimitSignalRange()
 893   {
 894     //empty signal
 895     Signal signal(m_SamplingPeriod);
 896     iterator it=begin();
 897     while(it != end())
 898       {
 899         signal.push_back(*it-floor(*it));
 900         it++;
 901       }
 902     return signal;
 903   }
 904
 905
 906   //---------------------------------------------------------------------
 907   double Signal::SSD(const Signal &sig2) const{
 908     if(sig2.size() != size()){
 909       cerr << "ERROR in Signal::SSD: signals don't have the same size" << endl;
 910       return -1;
 911     }
 912     double result=0.;
 913     for(unsigned int i=0;i<size();i++){
 914       result+=pow(sig2[i]-m_Data[i],2);
 915     }
 916     result/=size();
 917     result=sqrt(result);
 918     return result;
 919   }
 920   //---------------------------------------------------------------------
 921
 922
 923   //---------------------------------------------------------------------
 924   void Signal::AddValue(double v) {
 925     for(unsigned int i=0;i<size();i++) {
 926       m_Data[i] += v;
 927     }
 928   }
 929   //---------------------------------------------------------------------
 930
 931
 932   //---------------------------------------------------------------------
 933   void Signal::ComputeAugmentedSpace(clitk::Signal & outputX,
 934                                      clitk::Signal & outputY,
 935                                      unsigned int delay) const {
 936     if (size() <= delay) {
 937       std::cerr << "Error in signal length is " << size()
 938                 << " while delay is " << delay << " : too short. Abort." << std::endl;
 939       exit(0);
 940     }
 941     outputX.resize(size()-delay);
 942     outputY.resize(size()-delay);
 943     for(unsigned int i=0; i<outputX.size(); i++) {
 944       outputX[i] = m_Data[i+delay];
 945       outputY[i] = m_Data[i];
 946     }
 947   }
 948   //---------------------------------------------------------------------
 949
 950
 951   //---------------------------------------------------------------------
 952   void Signal::ResampleWithTimeWarpTo(clitk::Signal & output, bool lin) const {
 953     double sp = output.GetSamplingPeriod()/output.GetTotalTimeDuration()*GetTotalTimeDuration();
 954     // DD(GetTotalTimeDuration());
 955     //     DD(output.GetTotalTimeDuration());
 956     //     DD(sp);
 957     if (lin) {
 958       for(uint i=0; i<output.size(); i++) {
 959         output[i] = GetValueAtLin(i*sp);
 960       }
 961     }
 962     else {
 963       for(uint i=0; i<output.size(); i++) {
 964         output[i] = GetValueAtNN(i*sp);
 965       }
 966     }
 967   }
 968   //---------------------------------------------------------------------
 969
 970
 971   //---------------------------------------------------------------------
 972   double Signal::GetValueAtNN(double t) const {
 973     int index = (int)lrint(t/GetSamplingPeriod()); // closest point
 974     // DD(index);
 975     if ((index<0) || (index>=(int)size())) {
 976       FATAL("out of bound for time " << t << "in this signal having "
 977             << size() << " values with sampling " << GetSamplingPeriod() << std::endl);
 978     }
 979     return (*this)[index];
 980   }
 981   //---------------------------------------------------------------------
 982
 983
 984   //---------------------------------------------------------------------
 985   double Signal::GetValueAtLin(double t) const {
 986     double ti = t/GetSamplingPeriod();
 987        // DD(t);
 988 //          DD(ti);
 989     int index = (int)floor(ti); // floor point
 990     // DD(index);
 991     if (index<0) {
 992       FATAL("out of bound for time " << t << " in this signal having "
 993             << size() << " values with sampling " << GetSamplingPeriod() << std::endl);
 994     }
 995     if (index>=(int)size()-1) { // last value
 996       return (*this)[index];
 997     }
 998     // DD(index);
 999 //     DD(ti-index);
1000 //     DD((*this)[index]);
1001 //     DD((*this)[index+1]);
1002     return ((ti-index)*((*this)[index+1]) + (1.0-(ti-index))*((*this)[index]));
1003
1004   }
1005   //---------------------------------------------------------------------
1006
1007
1008   //---------------------------------------------------------------------
1009   void clitk::Signal::Shift(int d, int length) {
1010     std::cout << "Shift d=" << d << " length = " << length << std::endl;
1011     // int d = (int)lrint(s/GetSamplingPeriod()); // closest integer delta
1012     //     DD(d);
1013     clitk::Signal temp;
1014     temp.resize(length);//(uint)size());
1015     temp.SetSamplingPeriod(GetSamplingPeriod());
1016     for(uint i=0; i<temp.size(); i++) {
1017       int j = (i+d) % size();
1018       // DD(i);
1019       //       DD(j);
1020       temp[i] = (*this)[j];
1021     }
1022     CopyFrom(temp);
1023   }
1024   //---------------------------------------------------------------------
1025   //---------------------------------------------------------------------
1026   double Signal::UnormalizedXcov(const Signal &input2, Signal &output) const {
1027         unsigned int N;
1028         N=size();
1029         unsigned int max;
1030         max = 2*N;
1031         output.resize(max);
1032         assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1033
1034
1035         //calcul de la moyenne
1036         double mean1=0.,mean2=0.;
1037         mean1=this->GetGlobalMean();
1038         mean2=input2.GetGlobalMean();
1039
1040         for (unsigned int k= 0;k<N;k++){
1041           double value=0.;
1042           for (unsigned int n=0;n<=k;n++){
1043                 value=value+((*this)[n] - mean1)*(input2[n+N-k-1] - mean2);
1044           }
1045           output[k]=value;
1046         }
1047         for (unsigned int k=N;k<max-1;k++){
1048           double value=0.;
1049           for (unsigned int n=0;n<max-k-1;n++){
1050                 value=value+((*this)[n+k-N+1] - mean1)*(input2[n] - mean2);       //*input2??
1051           }
1052           output[k]=value;
1053         }
1054         double Coef;
1055         Coef = output[N-1];
1056         return Coef;
1057
1058   }
1059
1060  //---------------------------------------------------------------------
1061  //---------------------------------------------------------------------
1062   double Signal::UnormalizedCorrelation(const Signal &input2) const {
1063         unsigned int N;
1064         N=size();
1065         Signal output;
1066         output.resize(N);
1067         assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1068
1069         //calcul de la moyenne
1070         double mean1=0.,mean2=0.;
1071         mean1=this->GetGlobalMean();
1072         mean2=input2.GetGlobalMean();
1073         double value=0.;
1074         for (unsigned int n=0;n<N;n++){
1075                 value=value+((*this)[n] - mean1)*(input2[n] - mean2);
1076         }
1077         return value;
1078
1079   }
1080   //---------------------------------------------------------------------
1081  //---------------------------------------------------------------------
1082   double Signal::Correlation(const Signal &input2) const {
1083
1084         assert(size() == input2.size()); //les signaux à comparer doivent avoir la même taille
1085
1086         double cc11,cc22,cc;
1087         Signal input1(*this);
1088
1089         cc=this->UnormalizedCorrelation(input2);
1090         cc11=this->UnormalizedCorrelation(input1);
1091         cc22=input2.UnormalizedCorrelation (input2);
1092         cc=cc/sqrt(cc11*cc22);
1093         std::cout<<"coefficient de corrélation : "<<cc<<std::endl;
1094
1095         return cc;
1096
1097   }
1098
1099
1100
1101  //---------------------------------------------------------------------
1102  //---------------------------------------------------------------------
1103
1104   Signal Signal::Xcov(const Signal &input2 , int &shift, double &coef, double &cc) {
1105
1106         Signal temp_output,output;
1107         unsigned int N;
1108         N=size();
1109         unsigned int max;
1110         max = 2*N;
1111         output.resize(max);
1112         temp_output.resize(max);
1113         double cc11,cc22;
1114         Signal input1(*this);
1115
1116         cc=this->UnormalizedXcov(input2,output);
1117         cc11=this->UnormalizedCorrelation(input1);
1118         cc22=input2.UnormalizedCorrelation (input2);
1119
1120         unsigned int tailleOutput;
1121         tailleOutput=output.size();
1122         for (unsigned int k= 0;k<tailleOutput;k++){
1123           output[k]=output[k]/sqrt(cc11 * cc22 );
1124         }
1125
1126         coef=output[0];
1127         shift=0;
1128         for (unsigned int k= 0;k<tailleOutput;k++){
1129           if(output[k]>coef){
1130                 coef=output[k];
1131                 shift=k;
1132           }
1133         }
1134         cc=this->Correlation(input2);
1135         shift=shift-N+1;
1136         std::cout<<"shift: "<<shift<<std::endl;
1137         std::cout<<"XcovMax: "<<coef<<std::endl;
1138
1139
1140         return output;
1141   }
1142
1143   //---------------------------------------------------------------------
1144  //---------------------------------------------------------------------
1145
1146
1147   //  double Signal::Compare(Signal & sigRef) {
1148   //     double coeffCorrParam[5]={0.,0.,0.,0.,0.};
1149
1150   //     const_iterator itSig=begin();
1151   //     const_iterator itSigRef=sigRef.begin();//+offset;
1152   //     while(itSig!=end()){
1153   //       coeffCorrParam[0]+=*itSig;
1154   //       coeffCorrParam[1]+=*itSigRef;
1155   //       coeffCorrParam[2]+=(*itSig)*(*itSigRef);
1156   //       coeffCorrParam[3]+=(*itSig)*(*itSig);
1157   //       coeffCorrParam[4]+=(*itSigRef)*(*itSigRef);
1158   //       itSig++;itSigRef++;
1159   //     }
1160
1161   //     double coeffCorr = pow(size()*coeffCorrParam[2]-coeffCorrParam[0]*coeffCorrParam[1],2)/
1162   //       ((size()*coeffCorrParam[3]-pow(coeffCorrParam[0],2))*
1163   //        (size()*coeffCorrParam[4]-pow(coeffCorrParam[1],2)));
1164
1165   //     return coeffCorr;
1166   //   }
1167
1168   //   int Signal::DerivateSigne( const_iterator & it) const{
1169   //     int pos=-1;
1170   //     if(it==begin())
1171   //       pos=1;
1172   //     if((*it)==(*(it+pos)))
1173   //       return 0;
1174   //     else if((*it)<(*(it+pos)))
1175   //       return 1*pos;
1176   //     else // Case : ((*it)>(*(it+pos)))
1177   //       return -1*pos;
1178   //   }
1179
1180   //   void Signal::CenteredFiniteDifferences(Signal & result,int order,int* weights){
1181   //     const_iterator itSig=begin()+order;
1182   //     result.resize(size());
1183   //     iterator itDer=result.begin()+order;
1184   //     while(itSig!=end()-order){
1185   //       (*itDer)=0.;
1186   //       for(int i=-order;i<=order;i++){
1187   //    *itDer+=*(itSig-i)*weights[i+order];
1188   //       }
1189   //       itSig++;itDer++;
1190   //     }
1191   //   }
1192
1193   //   void Signal::FirstDerivate(Signal & result,int order){
1194   //     if(order==1){
1195   //       int weights[3]={-1,0,1};
1196   //       CenteredFiniteDifferences(result,order,weights);
1197   //     }
1198   //     else if(order==2){
1199   //       int weights[5]={1,-8,0,8,-1};
1200   //       CenteredFiniteDifferences(result,order,weights);
1201   //     }
1202   //   }
1203
1204   //   void Signal::SecondDerivate(Signal & result,int order){
1205   //     if(order==1){
1206   //       int weights[3]={1,-2,1};
1207   //       CenteredFiniteDifferences(result,order,weights);
1208   //     }
1209   //     else if(order==2){
1210   //       int weights[5]={-1,16,-30,16,-1};
1211   //       CenteredFiniteDifferences(result,order,weights);
1212   //     }
1213   //   }
1214
1215
1216
1217   //   void Signal::NormalizeMeanStdDev(double newMean,double newStdDev){
1218   //     iterator itSig=begin();
1219   //    double sum=0, sum2=0;
1220   //     while(itSig!=end()){
1221   //       sum += *itSig;
1222   //      sum2 += (*itSig) * (*itSig);
1223   //       itSig++;
1224   //     }
1225   //    double oldMean=sum/size();
1226   //    double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
1227
1228   //    double a = newStdDev/oldStdDev;
1229   //    double b = newMean - a * oldMean;
1230   //    itSig=begin();
1231   //    while(itSig!=end()){
1232   //      *itSig = a *(*itSig) + b;
1233   //       itSig++;
1234   //    }
1235   //   }
1236
1237
1238
1239   //   void Signal::print(ostream & os, const int level) const {
1240   //     os << "Size:" << m_Data.size() << endl;
1241   //     const_iterator it=m_Data.begin();
1242   //     while(it!=m_Data.end()){
1243   //       os << *it << endl;
1244   //       it++;
1245   //     }
1246   //   }
1247
1248
1249   //   //   }
1250
1251   //   //   istream& Signal::get(istream& is) {
1252   //   //     ERROR << "Signal::get NOT IMPLEMENTED";
1253   //   //     FATAL();
1254   //   //     return is;
1255   //   //   } ////
1256
1257   //   //   /** @b GridBase::put os
1258   //   //    * @param os
1259   //   //    * @return
1260   //   //    ***************************************************/
1261   //   //   ostream& Signal::put(ostream& os) const {
1262   //   //     print(os);
1263   //   //     return os;
1264   //   //   } ////
1265
1266
1267
1268   //   void Signal::Crop(unsigned int posmin, unsigned int posmax){
1269   //     if(posmin >= m_Data.size()) return;
1270   //     if(posmax >= m_Data.size()) posmax=m_Data.size();
1271   //     m_Data.erase(m_Data.begin()+posmax+1,m_Data.end());
1272   //     m_Data.erase(m_Data.begin(),m_Data.begin()+posmin);
1273   //   }
1274
1275   //   void Signal::LinearResample(const unsigned int newSize){
1276   //     SIGNAL newData;
1277   //     newData.push_back(front());
1278   //     double posInOld,leftWeight,rightWeight;
1279   //     int leftPos, rightPos;
1280   //     for(unsigned int i=1 ; i < newSize-1 ; i++){
1281   //       posInOld = (double)(i * (size()-1)) / (double)(newSize-1);
1282   //       leftPos = (int)floor(posInOld);
1283   //       rightPos = leftPos+1;
1284   //       leftWeight = (double)rightPos - posInOld;
1285   //       rightWeight =  posInOld - (double)leftPos;
1286   //       newData.push_back(m_Data[leftPos] * leftWeight + m_Data[rightPos] * rightWeight );
1287   //     }
1288
1289   //     newData.push_back(back());
1290   //     m_Data=newData;
1291   //   }
1292
1293
1294   //   int Signal::FreqToSamp(double freq){
1295   //     if(m_SamplingPeriod==-1.)
1296   //       cerr << "ERROR: you did not initialize the sampling period" << endl;
1297   //     return lrint(freq*(double)size()*m_SamplingPeriod);
1298   //   }
1299   //   double Signal::SampToFreq(int samp){
1300   //     if(m_SamplingPeriod==-1.)
1301   //       cerr << "ERROR: you did not initialize the sampling period" << endl;
1302   //     return ((double)(samp)/((double)size()*m_SamplingPeriod));
1303   //   }
1304
1305   //   //---------------------------------------------------------------------
1306   //   Signal Signal::limPhaseDE(eIPhaseValue, eEPhaseValue)
1307   //   {
1308
1309   //     //Create an empty signal
1310   //     phase.resize(size());
1311
1312   //     iterator phaseIt=initialPhaseValues.begin();
1313   //     iterator monPhaseIt=monPhase.begin();
1314   //     iterator extremaIt =begin();
1315
1316   //     while (extremaIt!= end())
1317   //     {
1318   //       if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
1319   //       else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
1320   //       extremaIt++; phaseIt++;monPhaseIt++;
1321   //     }
1322
1323   //   }
1324
1325
1326 }
1327
1328 #endif //#define CLITKSIGNAL