BALL: spectrum.h Source File

00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 // $Id: spectrum.h,v 1.14.18.6 2007-04-12 13:53:57 anne Exp $
00005 //
00006 
00007 #ifndef BALL_NMR_SPECTRUM_H
00008 #define BALL_NMR_SPECTRUM_H
00009 
00010 #ifndef BALL_NMR_PEAKLIST_H
00011 # include<BALL/NMR/peakList.h>
00012 #endif
00013 
00014 #ifndef BALL_DATATYPE_REGULARDATA1D_H
00015 # include<BALL/DATATYPE/regularData1D.h>
00016 #endif
00017 
00018 #ifndef BALL_DATATYPE_REGULARDATA2D_H
00019 # include<BALL/DATATYPE/regularData2D.h>
00020 #endif
00021 
00022 #ifndef BALL_DATATYPE_REGULARDATA3D_H
00023 # include<BALL/DATATYPE/regularData3D.h>
00024 #endif
00025 
00026 #ifdef BALL_HAS_FFTW
00027 # include <BALL/MATHS/complex.h>
00028 #ifndef BALL_MATHS_FFT1D_H
00029 # include <BALL/MATHS/FFT1D.h>
00030 #endif
00031 #ifndef BALL_MATHS_FFT2D_H
00032 # include <BALL/MATHS/FFT2D.h>
00033 #endif
00034 #endif
00035 
00036 
00037 using namespace std;
00038 
00039 namespace BALL 
00040 {
00047   template <typename DataT, typename PeakT, typename PositionT = typename PeakT::Position>
00048   class Spectrum
00049   {
00050     public:
00051     
00055 
00056     typedef DataT DataType;
00058     typedef PositionT PositionType;
00060     typedef PeakT PeakType;
00062     typedef typename DataT::Iterator Iterator;
00064     typedef typename DataT::ConstIterator ConstIterator;
00066 
00070       // ?????
00071     Spectrum()
00072       : data_(),
00073         sticks_(),
00074         spacing_(),
00075         min_(),
00076         max_()
00077     {}
00078 
00079     Spectrum(const DataType& data)
00080       : data_(data),
00081         sticks_(),
00082         spacing_(),
00083         min_(),
00084         max_()
00085     {}
00086 
00087     Spectrum(const std::vector<PeakType>& peaks, const PositionType& origin, const PositionType& dimension, const PositionType& spacing)
00088     {}
00089 
00092     virtual ~Spectrum() {}
00094 
00098 
00099     const DataType& getData() const;
00101     DataType& getData();
00102     
00103     const vector<float>& getHuInvariants() const;
00104     
00105     vector<float>& getHuInvariants();
00106 
00108 
00109     virtual void clear();
00110     virtual void clearSticks();
00111     virtual double difference(const Spectrum<DataT, PeakT, PositionT>& spectrum) const;
00112     virtual  Spectrum<DataT, PeakT, PositionT> differenceSpectrum(const Spectrum<DataT, PeakT, PositionT>& spectrum);
00113     virtual double earthMoversDistance(const Spectrum<DataT,PeakT, PositionT>& spectrum) const; 
00114     
00115     virtual void convertToGaussian(); 
00116     virtual void convertToLorentzian(); 
00117     virtual void computeAllMoments(int  moment_number);
00118 
00119     virtual void setSpacing(const PositionType& spacing);
00120     virtual PositionType getSpacing() const;
00121 
00122     virtual void setSticks(std::vector<PeakType> sticks) {sticks_ = sticks;};
00123     virtual std::vector<PeakType> getSticks() const {return sticks_;};
00124     
00125   
00126     // computes the integral over the fabs() of the spectrum
00127     virtual double getAbsIntegral() const;
00128     virtual void computeHuInvariants();
00129     virtual vector<float> computeHuInvariantsDifferences(vector<Spectrum<DataT, PeakT, PositionT> >& spectra);
00130 
00137     virtual double getFourierDifference(const Spectrum<DataT,PeakT, PositionT>& spectrum, float min_freq = 1e6, float max_freq = -1e6);
00138 
00142     virtual double getNormalMomentsDifference(Spectrum<DataT,PeakT, PositionT>& spectrum, int moment_number);
00143     virtual double getCentralMomentsDifference(Spectrum<DataT,PeakT, PositionT>& spectrum, int moment_number);
00144     virtual double getStandardizedMomentsDifference(Spectrum<DataT,PeakT, PositionT>& spectrum, int moment_number);
00145     
00147     std::vector<float> normal_moments;
00149     std::vector<float> central_moments;
00151     std::vector<float> standardized_moments;
00152   
00154     void binaryWrite(const String& filename);
00155 
00157     void binaryRead(const String& filename);
00158 
00159     protected:
00160       DataType              data_;
00161       std::vector<PeakType> sticks_;
00162       PositionType          spacing_; // rausschmeissen , wie auch immer 
00163       PositionType          min_;     // rausschmeissen , wie auch immer 
00164       PositionType          max_;     // rausschmeissen , wie auch immer 
00165       std::vector<float>    Hu_invariants_;
00166   };
00167 
00171   template <typename DataT, typename PeakT, typename PositionT>
00172   void Spectrum<DataT, PeakT, PositionT>::clear()
00173   {
00174     data_.clear();
00175     sticks_.clear();
00176   }
00177   
00178   template <typename DataT, typename PeakT, typename PositionT>
00179   void Spectrum<DataT, PeakT, PositionT>::clearSticks()
00180   {
00181     sticks_.clear();
00182   }
00183 
00184 
00185   // TODO: muss die hier stehen??? 
00188   template <typename DataT, typename PeakT, typename PositionT>
00189   double Spectrum<DataT, PeakT, PositionT>::difference(const Spectrum<DataT, PeakT, PositionT>& /* spectrum */) const
00190   {
00191     // ?????
00192     return 0.0;
00193   }
00194 
00197   template <typename DataT, typename PeakT, typename PositionT>
00198   typename Spectrum<DataT, PeakT, PositionT>::PositionType Spectrum<DataT, PeakT, PositionT>::getSpacing() const
00199   {
00200     return spacing_;
00201   }
00202 
00205   template <typename DataT, typename PeakT, typename PositionT>
00206   void Spectrum<DataT, PeakT, PositionT>::setSpacing(const typename Spectrum<DataT, PeakT, PositionT>::PositionType& spacing)
00207   {
00208     spacing_ = spacing;
00209   }
00210 
00214   template <typename DataT, typename PeakT, typename PositionT>
00215   double operator - (const Spectrum<DataT, PeakT, PositionT>& s1, const Spectrum<DataT, PeakT, PositionT>& s2)
00216   {
00217     return s1.difference(s2);
00218   }
00219   
00220   template <typename DataT, typename PeakT, typename PositionT>
00221   double Spectrum<DataT, PeakT, PositionT>::getNormalMomentsDifference(Spectrum<DataT, PeakT, PositionT>& spectrum, int moment_number)
00222   {
00223     if (normal_moments.size() != (Size)moment_number)
00224       computeAllMoments(moment_number);
00225     if (spectrum.normal_moments.size() != (Size)moment_number)
00226       spectrum.computeAllMoments(moment_number);
00227 
00228     double diff = 0.;
00229     for (int current_moment=0; current_moment<moment_number; current_moment++)
00230       diff += fabs(normal_moments[current_moment] - spectrum.normal_moments[current_moment]);
00231 
00232     return diff;
00233   }
00234 
00235   template <typename DataT, typename PeakT, typename PositionT>
00236   double Spectrum<DataT, PeakT, PositionT>::getCentralMomentsDifference(Spectrum<DataT, PeakT, PositionT>& spectrum, int moment_number)
00237   {
00238     if (central_moments.size() != (Size)moment_number)
00239       computeAllMoments(moment_number);
00240     if (spectrum.central_moments.size() != (Size)moment_number)
00241       spectrum.computeAllMoments(moment_number);
00242 
00243     double diff = 0.;
00244     for (int current_moment=0; current_moment<moment_number; current_moment++)
00245       diff += fabs(central_moments[current_moment] - spectrum.central_moments[current_moment]);
00246 
00247     return diff;
00248   }
00249 
00250   template <typename DataT, typename PeakT, typename PositionT>
00251   double Spectrum<DataT, PeakT, PositionT>::getStandardizedMomentsDifference(Spectrum<DataT, PeakT, PositionT>& spectrum, int moment_number)
00252   {
00253     if (standardized_moments.size() != (Size)moment_number)
00254       computeAllMoments(moment_number);
00255     if (spectrum.standardized_moments.size() != (Size)moment_number)
00256       spectrum.computeAllMoments(moment_number);
00257 
00258     double diff = 0.;
00259     for (int current_moment=0; current_moment<moment_number; current_moment++)
00260       diff += fabs(standardized_moments[current_moment] - spectrum.standardized_moments[current_moment]);
00261 
00262     return diff;
00263   } 
00264   
00265   template <typename DataT, typename PeakT, typename PositionT>
00266   const DataT& Spectrum<DataT, PeakT, PositionT>::getData() const
00267   {
00268     return data_;
00269   }
00270   
00271   template <typename DataT, typename PeakT, typename PositionT>
00272   DataT& Spectrum<DataT, PeakT, PositionT>::getData()
00273   {
00274     return data_;
00275   }
00276 /*  
00277   const vector<float>& getHuInvariants();
00278 */
00279 
00280   template <typename DataT, typename PeakT, typename PositionT>
00281   const vector<float>& Spectrum<DataT, PeakT, PositionT>::getHuInvariants() const
00282   {
00283     return Hu_invariants_;
00284   }
00285 
00286   template <typename DataT, typename PeakT, typename PositionT>
00287   vector<float>& Spectrum<DataT, PeakT, PositionT>::getHuInvariants() 
00288   {
00289     return Hu_invariants_;
00290   }
00291 
00292 
00293   
00294   template <typename DataT, typename PeakT, typename PositionT>
00295   void Spectrum<DataT, PeakT, PositionT>::computeHuInvariants()
00296   {
00297     Log.error()<< "computeHuInvariants() only implemented in 2D" << std::endl;
00298     return;
00299   }
00300 
00301   template <typename DataT, typename PeakT, typename PositionT>
00302   vector<float> Spectrum<DataT, PeakT, PositionT>::computeHuInvariantsDifferences(vector<Spectrum<DataT, PeakT, PositionT> >& /*spectra*/)
00303   {
00304     Log.error()<< "computeHuInvariantsDifferences() only implemented in 2D" << std::endl;
00305     std::vector<float> result;
00306     return result;
00307   }
00308 
00309   template <typename DataT, typename PeakT, typename PositionT>
00310   double Spectrum<DataT, PeakT, PositionT>::getFourierDifference(const Spectrum<DataT, PeakT, PositionT>& spectrum, float min_freq, float max_freq)
00311   {
00312     Log.error() << "getFourierDifference only implemented in 1D and 2D" << std::endl;
00313     return 0.;
00314   }
00315 
00320 
00321   typedef Spectrum<RegularData1D, Peak1D> Spectrum1D;
00322 
00324   typedef Spectrum<RegularData2D, Peak2D> Spectrum2D;
00325 
00327   typedef Spectrum<RegularData3D, Peak3D> Spectrum3D;
00329   
00330   template <typename DataT, typename PeakT, typename PositionT>
00331   std::ostream& operator << (std::ostream& os, const Spectrum<DataT, PeakT, PositionT>& spectrum);
00332   
00333   template <typename DataT, typename PeakT, typename PositionT>
00334   std::istream& operator >> (std::istream& is, Spectrum<DataT, PeakT, PositionT>& spectrum);
00335 
00336 # ifndef BALL_NO_INLINE_FUNCTIONS
00337 #   include <BALL/NMR/spectrum.iC>
00338 # endif
00339 } // namespace BALL
00340 
00341 
00342 #endif // BALL_NMR_SPECTRUM_H