DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
MeshFromDiscreteFiberODF.cxx
Go to the documentation of this file.
1 
18 #include "vtkPolyDataWriter.h"
19 
20 #include "MeshFromDiscreteFiberODFCLP.h"
23 
26 
28 #include "utl.h"
29 #include "utlVTK.h"
30 #include "vtkPolyDataViewer.h"
31 
38 int
39 main (int argc, char const* argv[])
40 {
41  PARSE_ARGS;
42 
43  // Time Probe
44  itk::TimeProbe clock;
45 
46  utlGlobalException(_BoxView.size()!=6, "need 6 parameters in --box");
47  utlGlobalException(_SliceView.size()!=3, "need 3 parameters in --box");
48  utlGlobalException(_Flip.size()!=3, "need 3 parameters in --flip");
49 
50  // Define Variables
51  typedef double PixelType;
54 
55  InputImageType::Pointer inputImage = InputImageType::New();
56  itk::ReadImage<InputImageType, ReaderType>(_InputFile, inputImage);
57 
58  // Create mesh
60  MeshCreatorType::Pointer filter = MeshCreatorType::New();
61  filter->SetScale( _Scale );
62  if (_MaskFileArg.isSet())
63  filter->SetMaskImage(_MaskFile);
64 
65  if (_Normalization=="NONE") filter->SetNormalization( MeshCreatorType::NONE );
66  if (_Normalization=="MIN_MAX") filter->SetNormalization( MeshCreatorType::MIN_MAX );
67  if (_Normalization=="UNIT_MAX") filter->SetNormalization( MeshCreatorType::UNIT_MAX );
68  if (_Normalization=="UNIT_INTEGRAL") filter->SetNormalization( MeshCreatorType::UNIT_INTEGRAL );
69 
70  filter->SetPow(_Pow);
71  filter->SetRemoveNegativeValues(_RemoveNegativeValuesArg.isSet());
72  filter->SetStretch(!_NoStretchArg.isSet());
73  filter->SetInput( inputImage );
74 
75  filter->SetBoxView(_BoxView[0], _BoxView[1], _BoxView[2], _BoxView[3], _BoxView[4], _BoxView[5]);
76  filter->SetSliceView(_SliceView[0], _SliceView[1], _SliceView[2]);
77  filter->SetFlip(_Flip[0], _Flip[1], _Flip[2]);
78 
79  if (_BasicShape == "TETRAHEDRON") { filter->SetBasicShape(MeshCreatorType::SphereTessellatorType::TETRAHEDRON); }
80  if (_BasicShape == "OCTAHEDRON") { filter->SetBasicShape(MeshCreatorType::SphereTessellatorType::OCTAHEDRON); }
81  if (_BasicShape == "ICOSAHEDRON") { filter->SetBasicShape(MeshCreatorType::SphereTessellatorType::ICOSAHEDRON); }
82  filter->SetTessellationOrder(_TessellationOrder);
83 
84  if (_ColorScheme == "DIRECTION") { filter->SetColorScheme(MeshCreatorType::DIRECTION); }
85  if (_ColorScheme == "MAGNITUDE") { filter->SetColorScheme(MeshCreatorType::MAGNITUDE); }
86 
88  if (_ShowProgressArg.isSet())
89  filter->AddObserver( itk::ProgressEvent(), observer );
90  if (_SingleThreadArg.isSet())
91  filter->SetNumberOfThreads(1);
92  filter->SetDebug(_DebugArg.isSet());
93 
95  if (_OrientationsFileArg.isSet())
96  gradBasis = utl::ReadGrad<double>(_OrientationsFile, DIRECTION_NODUPLICATE, CARTESIAN_TO_CARTESIAN);
97  else
98  gradBasis = utl::ReadGrad<double>(4, DIRECTION_NODUPLICATE, CARTESIAN_TO_CARTESIAN);
99  utlGlobalException((_OutputType=="DWI" || _OutputType=="EAP") && _Radius<0, "need to set _Radius ");
100  if (_Radius>0)
101  {
102  filter->SetRadius(_Radius);
103  }
104 
105  utlGlobalException(_TensorEigenValues.size()!=2 && _TensorEigenValues.size()!=3, "wrong size of --tensor");
106 
107  typedef itk::TensorBasisMatrixGenerator<double> TensorBasisGeneratorType;
108  TensorBasisGeneratorType::Pointer tensorBasisGenerator = NULL;
109  if (_BasisType=="TENSOR")
110  {
111  utlGlobalException(_TensorEigenValues[0]<_TensorEigenValues[1]-1e-10, "The first eigenvalue should be the largest");
112  tensorBasisGenerator = TensorBasisGeneratorType::New();
113  tensorBasisGenerator->SetEigenValue1(_TensorEigenValues[0]);
114  tensorBasisGenerator->SetEigenValue2(_TensorEigenValues[1]);
115  if (_TensorEigenValues.size()==2)
116  tensorBasisGenerator->SetEigenValue3(_TensorEigenValues[1]);
117  else
118  {
119  utlGlobalException(_TensorEigenValues[2]<_TensorEigenValues[3]-1e-10, "The second eigenvalue should be more than the third one");
120  tensorBasisGenerator->SetEigenValue3(_TensorEigenValues[2]);
121  }
122  if (_IsotroipicEigenValue>0)
123  tensorBasisGenerator->SetEigenValueISO(_IsotroipicEigenValue);
124 
125  tensorBasisGenerator->SetBasisOrientations(gradBasis);
126  tensorBasisGenerator->SetUseIsotropicTerm( _UseIsotropicComponentArg.isSet() );
127  // tensorBasisGenerator->Print(std::cout<<"basis");
128  filter->SetBasisMatrixGenerator(tensorBasisGenerator);
129  }
130  else if (_BasisType=="CYLINDER")
131  {
132  utlGlobalException(true, "TODO");
133  }
134 
135  filter->GetBasisMatrixGenerator()->SetODFOrder(_ODFOrder);
136  filter->GetBasisMatrixGenerator()->GetSamplingSchemeQSpace()->SetTau(_Tau);
137  filter->GetBasisMatrixGenerator()->GetSamplingSchemeRSpace()->SetTau(_Tau);
138  if (_OutputType=="DWI")
139  filter->GetBasisMatrixGenerator()->SetOutputType(TensorBasisGeneratorType::DWI);
140  else if (_OutputType=="EAP")
141  filter->GetBasisMatrixGenerator()->SetOutputType(TensorBasisGeneratorType::EAP);
142  else if (_OutputType=="ODF")
143  filter->GetBasisMatrixGenerator()->SetOutputType(TensorBasisGeneratorType::ODF);
144 
145  std::cout << "Generating mesh ... " << std::flush;
146  clock.Start();
147  filter->Update();
148  clock.Stop();
149  std::cout << clock.GetMean() << "s elapsed" << std::endl;
150 
151  MeshCreatorType::OutputMeshPolyDataType* mesh = filter->GetOutput();
152 
153  if (_OutputFileArg.isSet())
154  {
155  utl::WriteVtkPolyData(mesh, _OutputFile);
156  }
157  else
158  {
159  utlGlobalException(_WindowSize.size()!=2, "wrong window size");
160  utlGlobalException(_Angle.size()!=2, "wrong angle size");
161  utlGlobalException(_BackgroundColor.size()!=3, "wrong size of background color");
162  vtk::VisualizePolyData(mesh, _Angle, _WindowSize, !_NoNormalArg.isSet(), !_NoLightingArg.isSet(), _Zoom, _PNGFile, _BackgroundColor);
163  }
164 
165  return EXIT_SUCCESS;
166 }
Compute mesh from discrete fiber ODF represented by a set of basis functions.
helper functions specifically used in dmritool
Created "06-08-2017.
Created "06-06-2017.
utl_shared_ptr< MatrixType > MatrixPointer
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
void VisualizePolyData(vtkPolyData *mesh, const std::vector< double > &angle={0.0, 0.0}, const std::vector< int > &windowSize={600, 600}, const bool useNormal=true, const bool lighting=true, const double zoom=1.0, const std::string &pngfile="", const std::vector< double > &bgColor={0, 0, 0})
int main(int argc, char const *argv[])
Mesh From Discrete Fiber ODF Create mesh (a spherical profile of DWI/EAP/ODF) from discrete fiber ODF...
An n-dimensional vector image with a sparse memory model.
void WriteVtkPolyData(vtkPolyData *mesh, const std::string &filename)
Definition: utlVTK.h:45