BALL
1.4.1
|
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