#ifndef CLITKSIGNAL_CXX
#define CLITKSIGNAL_CXX
#include "clitkSignal.h"
namespace clitk {
//---------------------------------------------------------------------
void Signal::Read(string fileName){
ifstream signalStream(fileName.c_str());
SignalValueType point;
if(!signalStream.is_open()){
std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
return;
}
skipComment(signalStream);
while(!signalStream.eof()) {
skipComment(signalStream);
signalStream >> point;
skipComment(signalStream);
m_Data.push_back(point);
// cout << point << endl;
}
signalStream.close();
SetSamplingPeriod(1.);
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
void Signal::Read(string fileName, int col){
ifstream signalStream(fileName.c_str());
SignalValueType point;
if(!signalStream.is_open()){
std::cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Read" << std::endl;
return;
}
skipComment(signalStream);
while(!signalStream.eof()) {
skipComment(signalStream);
// Read one line
std::string line;
std::getline(signalStream, line);
// Get column nb col
istringstream iss(line);
for(int i=0; i
> sub;
}
iss >> point;
if (!iss) {
std::cerr << "ERROR: no col n" << col << " in the line '" << line << "' ?" << std::endl;
exit(0);
}
skipComment(signalStream);
m_Data.push_back(point);
// cout << point << endl;
}
signalStream.close();
SetSamplingPeriod(1.);
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
//Convert 1D image to signal
Signal Signal::ConvertImageToSignal( Signal::ImageType::Pointer image)
{
//empty signal
Signal signal;
//make an image iterator
itk::ImageRegionConstIterator it(image,image->GetLargestPossibleRegion());
it.Begin();
//copy
while(!it.IsAtEnd())
{
signal.push_back(it.Get());
++it;
}
//Spacing
signal.SetSamplingPeriod(image->GetSpacing()[0]);
return signal;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
//Convert 1D signal to image
Signal::ImageType::Pointer Signal::ConvertSignalToImage(Signal signal)
{
//empty image
ImageType::Pointer image =ImageType::New();
ImageType::RegionType region;
ImageType::RegionType::IndexType index;
index[0]=0;
ImageType::RegionType::SizeType size;
size[0]=signal.size();
region.SetIndex(index);
region.SetSize(size);
image->SetRegions(region);
image->Allocate();
//make an image iterator
itk::ImageRegionIterator mIt(image,image->GetLargestPossibleRegion());
mIt.Begin();
//make a signal iterator
Signal::const_iterator sIt=signal.begin();
//copy
while(sIt!=signal.end())
{
mIt.Set(*sIt);
sIt++;++mIt;
}
//spacing
ImageType::SpacingType spacing;
spacing[0]=signal.GetSamplingPeriod();
image->SetSpacing(spacing);
return image;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
void Signal::Write(const string fileName){
ofstream signalStream(fileName.c_str());
if(!signalStream.is_open()){
cerr << "ERROR: Couldn't open file " << fileName << " in Signal::Write" << endl;
return;
}
iterator it=begin();
while(it!=end()) {
signalStream << *it << endl;
it++;
}
signalStream.close();
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal & Signal::operator/=(Signal & div){
if(size()!= div.size())
{
std::cerr << "Error: signal size must be the same!" << std::endl;
return *this;
}
iterator it=begin();
iterator itD;
itD=div.begin();
while(it!=end())
{
if(*itD!=0)
*it/=*itD;
else
cerr << "Division by 0 in operator/= skipped" << endl;
it++;itD++;
}
return *this;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal & Signal::operator*=(Signal & mul){
if(size()!= mul.size())
{
std::cerr << "Error: signal size must be the same!" << std::endl;
return *this;
}
iterator it=begin();
iterator itD;
itD=mul.begin();
while(it!=end()){
*it *= *itD;
it++;itD++;
}
return *this;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::Normalize(double newMin,double newMax){
Signal temp (m_SamplingPeriod);
vector extrema=GetGlobalMinMax();
iterator itSig=begin();
while(itSig!=end()){
temp.push_back( ((*itSig)-extrema[0])*(newMax-newMin)/(extrema[1]-extrema[0]) + newMin );
itSig++;
}
return temp;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
vector Signal::GetGlobalMinMax() const {
vector extrema(2);
if(size()==0){
cerr << "ERROR: GetExtrema / No signal" << endl;
return extrema;
}
extrema[0]=m_Data[0];
extrema[1]=m_Data[0];
for(unsigned int i=1;im_Data[i]) extrema[0]=m_Data[i];
if(extrema[1](i)-static_cast(length)); j<=min(size(), i+length); j++)
{
accumulator+=m_Data[j];
scale++;
}
temp.push_back(accumulator/scale);
}
return temp;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::GaussLikeFilter ( ) {
Signal temp(m_SamplingPeriod);
if (size()<2)
return *this;
else
{
//first sample: mirrorring BC
temp.push_back((2.*m_Data[0]+2*m_Data[1])/4.);
//middle samples
for (unsigned int i=1; i (size())*sampPeriod);
//forward fft with empty fft
if(fft.size()==0) OneDForwardFourier(*this, fft);
//remove the frequencies
for(unsigned int i=0;i(0.,0.);
//backward with remaining frequencies
OneDBackwardFourier(fft,temp);
return temp;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::LowPassFilter (double sampPeriod, double cutOffFrequency )
{
//output
Signal temp(m_SamplingPeriod);
temp.resize(size());
//the fft
SIGNAL_FFT_TYPE fft;
//calculate the cut off frequency
unsigned int samp=lrint(cutOffFrequency*static_cast(size())*sampPeriod);
//forward fft with empty fft
if(fft.size()==0) OneDForwardFourier(*this, fft);
unsigned int fsize=fft.size();
//remove the frequencies
unsigned int limit=min (samp, fsize);
for(unsigned int i=limit;i(0.,0.);
//backward with remaining frequencies
OneDBackwardFourier(fft,temp);
return temp;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
void Signal::OneDForwardFourier(const Signal& input, SIGNAL_FFT_TYPE & fft)
{
//Create output array
fft.resize(input.size()/2+1);
//Temp copy
double *tempCopy=new double[size()];
copy(begin(), end(), tempCopy);
//Forward Fourier Transform
fftw_plan p;
p=fftw_plan_dft_r2c_1d(size(),tempCopy,reinterpret_cast(&(fft[0])),FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
//delete tempCopy;
return;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
void Signal::OneDBackwardFourier(SIGNAL_FFT_TYPE & fft, Signal &output)
{
//Backward
fftw_plan p;
p=fftw_plan_dft_c2r_1d(output.size(),reinterpret_cast(&(fft[0])),&(output[0]),FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
vector::iterator it=output.begin();
while(it!=output.end()){
*it /= (double)output.size();
it++;
}
return;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
double Signal::MaxFreq(const Signal &sig,SIGNAL_FFT_TYPE & fft)
{
if(fft.size()==0) OneDForwardFourier(sig,fft);
int posMax=1;
double amplitude, amplitudeMax=abs(fft[1]);
for(unsigned int i=1;iamplitudeMax){
posMax=i;
amplitudeMax=amplitude;
}
}
return ((double)(posMax)/((double)sig.size()*sig.GetSamplingPeriod()));
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::DetectLocalExtrema(unsigned int width)
{
Signal temp(m_SamplingPeriod);
bool isMin, isMax;
unsigned int upper, lower;
//has to be at least 1
width=max(static_cast(width), 1);
for(unsigned int i=0 ; i < size() ; i++){
isMax=true;
isMin=true;
for(unsigned int j=1; j< width+1; j++)
{
//set the boundaries
upper = min( size(), i+j);
lower = max( static_cast(0), (int)i-(int)j);
//check if max
if( ! (m_Data[i] >= m_Data[lower] && m_Data[i] >= m_Data[upper]))isMax=false;
//check if min
if( ! (m_Data[i]<= m_Data[lower] && m_Data[i] <= m_Data[upper])) isMin=false;
//if neither, go to the next value
if( (!isMax) && (!isMin)) break;
}
if (isMax) temp.push_back(1);
else if (isMin) temp.push_back(0);
else temp.push_back(0.5);
}
return temp;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::LimPhase()
{
//phase is defined as going from 0 to 1 linearly between end expiration
Signal phase(m_SamplingPeriod);
phase.resize(size());
unsigned int firstBeginCycle=0;
unsigned int firstEndCycle=0;
unsigned int beginCycle=0;
//=========================================
//Search the first value in extrema not 0.5
while(m_Data[beginCycle]!=1)
{
beginCycle++;
}
//We search the corresponding end
unsigned int endCycle=beginCycle+1;
while(endCycle != size() && m_Data[endCycle]!=1){
endCycle++;
}
//============================================================
//Calculate phase at the beginning (before the first extremum)
for(unsigned int i=0 ; i 1)
{
phase[i] = (double)(i-0)/(double)(beginCycle-0);
}
//copy the phase values later
else
{
firstBeginCycle=beginCycle;
firstEndCycle=endCycle;
}
}
//===================================================================
//Middle part
while(endCycle != size()){
//fill between extrema
for(unsigned int i=beginCycle ; i 1){
//make the last part go till 1
phase[i] = (double)(i-endCycle)/(double)(size()-endCycle);
}
//the last part is shorter, copy the last cycle values
else{
phase[i] = phase[i -(endCycle-beginCycle)];
}
}
//===================================================================
//check it some remains to be copied in the beginning
if (firstBeginCycle!=0)
{
for(unsigned int i=0 ; i 1)
{
phase[i] = (double)(i-0)/(double)(beginCycle-0);
}
//copy the phase values later
else
{
firstBeginCycle=beginCycle;
firstEndCycle=endCycle;
}
}
//===================================================================
//Middle part
while(endCycle != size()){
cycleCounter++;
//fill between extrema
for(unsigned int i=beginCycle ; i 1){
//make the last part go till 1
phase[i] = (double)cycleCounter+(double)(i-endCycle)/(double)(size()-endCycle);
}
//the last part is shorter, copy the last cycle values
else{
phase[i] = phase[i -(endCycle-beginCycle)]+1;
}
}
//===================================================================
//check it some remains to be copied in the beginning
if (firstBeginCycle!=0)
{
for(unsigned int i=0 ; iMonPhase();
//Create an empty signal
Signal phase(size(), -1);
phase.SetSamplingPeriod(m_SamplingPeriod);
//Fill in the values at the extrema position
iterator phaseIt=phase.begin();
iterator monPhaseIt=monPhase.begin();
iterator extremaIt =begin();
while (extremaIt!= end())
{
if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
extremaIt++; phaseIt++;monPhaseIt++;
}
return phase;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::LinearlyInterpolateScatteredValues()
{
//Linearly interpolate the values in between
unsigned int i=0;
unsigned int beginCycle=0;
unsigned int endCycle=0;
//Create a new signal
Signal temp(size(),-1);
temp.SetSamplingPeriod(m_SamplingPeriod);
//start from the first value
while (m_Data[i]==-1)i++;
beginCycle=i;
i++;
//Go to the next
while ( (m_Data[i]==-1) && (i PointSetType;
// PointSetType::Pointer pointSet=PointSetType::New();
// typedef PointSetType::PointType PointType;
// unsigned int i=0;
// unsigned int pointIndex=0;
// while (i< size())
// {
// if(m_Data[i]!=-1)
// {
// PointType p;
// p[0]= i;//JV spacing is allways 1
// pointSet->SetPoint( pointIndex, p );
// pointSet->SetPointData( pointIndex, m_Data[i] );
// pointIndex++;
// }
// i++;
// }
// //define the output signal properties
// ImageType::RegionType::SizeType outputSize;
// outputSize[0]= size();
// ImageType::PointType outputOrigin;
// outputOrigin[0]=0.0;//JV may need to be changed
// ImageType::SpacingType outputSpacing;
// outputSpacing[0]=1; //JV add relation to the original signal spacing
// //Convert
// typedef itk::BSplineScatteredDataPointSetToImageFilter< PointSetType, VectorImageType > PointSetToImageFilterType;
// PointSetToImageFilterType::Pointer pointSetToImageFilter= PointSetToImageFilterType::New();
// pointSetToImageFilter->SetInput(pointSet);
// pointSetToImageFilter->SetSplineOrder(splineOrder);//JV
// pointSetToImageFilter->SetSize(outputSize);
// pointSetToImageFilter->SetOrigin(outputOrigin);
// pointSetToImageFilter->SetSpacing(outputSpacing);
// //Convert to
// itk::FixedArray num;
// num[0]=numberOfControlPoints;
// pointSetToImageFilter->SetNumberOfControlPoints(num);//JV
// pointSetToImageFilter->Update();
// VectorImageType::Pointer approximatedSignal=pointSetToImageFilter->GetOutput();
// //Convert and return
// return ConvertVectorImageToSignal(approximatedSignal);
// }
// //---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::ConvertVectorImageToSignal(VectorImageType::Pointer image)
{
//empty signal
Signal signal;
//make an image iterator
itk::ImageRegionConstIterator it(image,image->GetLargestPossibleRegion());
it.Begin();
//copy
while(!it.IsAtEnd())
{
signal.push_back(it.Get()[0]);
++it;
}
//Spacing
signal.SetSamplingPeriod(image->GetSpacing()[0]);
return signal;
}
//---------------------------------------------------------------------
//---------------------------------------------------------------------
Signal Signal::LimitSignalRange()
{
//empty signal
Signal signal(m_SamplingPeriod);
iterator it=begin();
while(it != end())
{
signal.push_back(*it-floor(*it));
it++;
}
return signal;
}
//---------------------------------------------------------------------
double Signal::SSD(const Signal &sig2) const{
if(sig2.size() != size()){
cerr << "ERROR in Signal::SSD: signals don't have the same size" << endl;
return -1;
}
double result=0.;
for(unsigned int i=0;i(*(it+pos)))
// return -1*pos;
// }
// void Signal::CenteredFiniteDifferences(Signal & result,int order,int* weights){
// const_iterator itSig=begin()+order;
// result.resize(size());
// iterator itDer=result.begin()+order;
// while(itSig!=end()-order){
// (*itDer)=0.;
// for(int i=-order;i<=order;i++){
// *itDer+=*(itSig-i)*weights[i+order];
// }
// itSig++;itDer++;
// }
// }
// void Signal::FirstDerivate(Signal & result,int order){
// if(order==1){
// int weights[3]={-1,0,1};
// CenteredFiniteDifferences(result,order,weights);
// }
// else if(order==2){
// int weights[5]={1,-8,0,8,-1};
// CenteredFiniteDifferences(result,order,weights);
// }
// }
// void Signal::SecondDerivate(Signal & result,int order){
// if(order==1){
// int weights[3]={1,-2,1};
// CenteredFiniteDifferences(result,order,weights);
// }
// else if(order==2){
// int weights[5]={-1,16,-30,16,-1};
// CenteredFiniteDifferences(result,order,weights);
// }
// }
// void Signal::NormalizeMeanStdDev(double newMean,double newStdDev){
// iterator itSig=begin();
// double sum=0, sum2=0;
// while(itSig!=end()){
// sum += *itSig;
// sum2 += (*itSig) * (*itSig);
// itSig++;
// }
// double oldMean=sum/size();
// double oldStdDev=sqrt(sum2/size()-oldMean*oldMean);
// double a = newStdDev/oldStdDev;
// double b = newMean - a * oldMean;
// itSig=begin();
// while(itSig!=end()){
// *itSig = a *(*itSig) + b;
// itSig++;
// }
// }
// void Signal::print(ostream & os, const int level) const {
// os << "Size:" << m_Data.size() << endl;
// const_iterator it=m_Data.begin();
// while(it!=m_Data.end()){
// os << *it << endl;
// it++;
// }
// }
// // }
// // istream& Signal::get(istream& is) {
// // ERROR << "Signal::get NOT IMPLEMENTED";
// // FATAL();
// // return is;
// // } ////
// // /** @b GridBase::put os
// // * @param os
// // * @return
// // ***************************************************/
// // ostream& Signal::put(ostream& os) const {
// // print(os);
// // return os;
// // } ////
// void Signal::Crop(unsigned int posmin, unsigned int posmax){
// if(posmin >= m_Data.size()) return;
// if(posmax >= m_Data.size()) posmax=m_Data.size();
// m_Data.erase(m_Data.begin()+posmax+1,m_Data.end());
// m_Data.erase(m_Data.begin(),m_Data.begin()+posmin);
// }
// void Signal::LinearResample(const unsigned int newSize){
// SIGNAL newData;
// newData.push_back(front());
// double posInOld,leftWeight,rightWeight;
// int leftPos, rightPos;
// for(unsigned int i=1 ; i < newSize-1 ; i++){
// posInOld = (double)(i * (size()-1)) / (double)(newSize-1);
// leftPos = (int)floor(posInOld);
// rightPos = leftPos+1;
// leftWeight = (double)rightPos - posInOld;
// rightWeight = posInOld - (double)leftPos;
// newData.push_back(m_Data[leftPos] * leftWeight + m_Data[rightPos] * rightWeight );
// }
// newData.push_back(back());
// m_Data=newData;
// }
// int Signal::FreqToSamp(double freq){
// if(m_SamplingPeriod==-1.)
// cerr << "ERROR: you did not initialize the sampling period" << endl;
// return lrint(freq*(double)size()*m_SamplingPeriod);
// }
// double Signal::SampToFreq(int samp){
// if(m_SamplingPeriod==-1.)
// cerr << "ERROR: you did not initialize the sampling period" << endl;
// return ((double)(samp)/((double)size()*m_SamplingPeriod));
// }
// //---------------------------------------------------------------------
// Signal Signal::limPhaseDE(eIPhaseValue, eEPhaseValue)
// {
// //Create an empty signal
// phase.resize(size());
// iterator phaseIt=initialPhaseValues.begin();
// iterator monPhaseIt=monPhase.begin();
// iterator extremaIt =begin();
// while (extremaIt!= end())
// {
// if (*extremaIt==0.) *phaseIt=eIPhaseValue+floor(*monPhaseIt);
// else if (*extremaIt==1.) *phaseIt=eEPhaseValue+floor(*monPhaseIt);
// extremaIt++; phaseIt++;monPhaseIt++;
// }
// }
}
#endif //#define CLITKSIGNAL