00001
00006 #ifndef TRACTOGRAPHY_H_
00007 #define TRACTOGRAPHY_H_
00008
00009 #include <string>
00010 #include <vector>
00011 #include "fiber.h"
00012 #include "seed.h"
00013
00014 #include <vnl/vnl_matrix.h>
00015 #include <vnl/algo/vnl_determinant.h>
00016 #include <vnl/vnl_vector_ref.h>
00017 #include <vnl/algo/vnl_svd.h>
00018 #include <vnl/algo/vnl_qr.h>
00019
00020 class ISignalData;
00021
00026 class Tractography
00027 {
00028
00029 public:
00030
00032 enum model_type { _1T,
00033 _1T_FW,
00034 _1T_FULL,
00035 _1T_FW_FULL,
00036 _2T,
00037 _2T_FW,
00038 _2T_FULL,
00039 _2T_FW_FULL,
00040 _3T,
00041 _3T_FULL
00042 };
00043
00045 Tractography(FilterModel *model, model_type filter_model_type,
00046
00047 const std::string& output_file, const std::string &output_file_with_second_tensor,
00048 const bool record_fa, const bool record_nmse, const bool record_trace,
00049 const bool record_state, const bool record_cov, const bool record_free_water, const bool record_tensors,
00050 const bool transform_position, const bool store_glyphs, const bool branchesOnly,
00051
00052 const double fa_min, const double ga_min, const double seedFALimit,
00053 const int num_tensors, const int seeds_per_voxel,
00054 const double minBranchingAngle, const double maxBranchingAngle,
00055 const bool is_full_model, const bool free_water,
00056 const double stepLength, const double maxHalfFiberLength,
00057 const std::vector<int>& labels,
00058
00059 double p0, double sigma_signal, double sigma_mask,
00060 double min_radius, double full_brain_ga_min,
00061
00062 const int num_threads
00063 );
00064
00066 ~Tractography();
00067
00072 bool LoadFiles(const std::string& data_file,
00073 const std::string& seed_file,
00074 const std::string& mask_file,
00075 const bool normalized_DWI_data,
00076 const bool output_normalized_DWI_data
00077 ) ;
00078
00084 void Init(std::vector<SeedPointInfo>& seed_infos) ;
00085
00087 void Run() ;
00088
00092 void Follow3T(const int thread_id, const size_t seed_index, const SeedPointInfo& seed, Fiber& fiber, bool is_branching,
00093 std::vector<SeedPointInfo>& branching_seeds, std::vector<BranchingSeedAffiliation>& branching_seed_affiliation) ;
00097 void Follow2T(const int thread_id, const size_t seed_index, const SeedPointInfo& seed, Fiber& fiber, bool is_branching,
00098 std::vector<SeedPointInfo>& branching_seeds, std::vector<BranchingSeedAffiliation>& branching_seed_affiliation) ;
00102 void Follow1T(const int thread_id, const SeedPointInfo& seed, Fiber& fiber);
00103
00104 private:
00109 void UnpackTensor(const std::vector<double>& b,
00110 const std::vector<vec_t>& u,
00111 std::vector<std::vector<double> >& s,
00112 std::vector<std::vector<double> >& ret);
00113
00115 void Step3T(const int thread_id, vec_t& x, vec_t& m1, vec_t& l1, vec_t& m2, vec_t& l2, vec_t& m3,
00116 vec_t& l3, double& fa, double& fa2, State& state,
00117 vnl_matrix<double>& covariance, double& dNormMSE, double& trace, double& trace2);
00119 void Step2T(const int thread_id, vec_t& x, vec_t& m1, vec_t& l1, vec_t& m2, vec_t& l2, double& fa, double& fa2,
00120 State& state, vnl_matrix<double>& covariance, double& dNormMSE, double& trace, double& trace2);
00122 void Step1T(const int thread_id, vec_t& x, double& fa, State& state,
00123 vnl_matrix<double>& covariance, double& dNormMSE, double& trace);
00124
00125
00130 void SwapState3T(State& state, vnl_matrix<double>& covariance, int i);
00135 void SwapState2T(State& state, vnl_matrix<double>& covariance);
00136
00141 void Record(const vec_t& x, double fa, double fa2, const State& state,
00142 const vnl_matrix<double> p, Fiber& fiber, double dNormMSE, double trace, double trace2);
00143
00145 std::vector<UnscentedKalmanFilter *> _ukf ;
00146
00148 ISignalData *_signal_data;
00150 FilterModel *_model;
00151
00153 model_type _filter_model_type ;
00154
00156 const std::string _output_file ;
00158 const std::string _output_file_with_second_tensor ;
00160 const bool _record_fa ;
00165 const bool _record_nmse ;
00167 const bool _record_trace ;
00169 const bool _record_state ;
00171 const bool _record_cov ;
00173 const bool _record_free_water ;
00178 const bool _record_tensors;
00182 const bool _transform_position ;
00184 const bool _store_glyphs ;
00186 bool _branches_only;
00187
00188
00189 bool _is_branching ;
00190 const double _p0 ;
00191 const double _sigma_signal;
00192 const double _sigma_mask ;
00193 const double _min_radius ;
00194 const double _full_brain_ga_min ;
00196 const int _max_length ;
00197 bool _full_brain ;
00199 int _nPosFreeWater;
00200
00201
00202 const double _fa_min;
00203 const double _ga_min;
00204 const double _seedFALimit ;
00205 const int _num_tensors ;
00206 const int _seeds_per_voxel;
00207 double _cos_theta_min ;
00208 double _cos_theta_max ;
00209 const bool _is_full_model ;
00210 const bool _free_water;
00211 const double _stepLength ;
00212 const std::vector<int> _labels ;
00213
00214
00215 const int _num_threads ;
00216 };
00217
00218 #endif // TRACTOGRAPHY_H_