common/clitkSignal.cxx

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