DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
OrientationsViewer.cxx
Go to the documentation of this file.
1 
18 #include <vtkCellArray.h>
19 #include <vtkProperty.h>
20 #include <vtkDataSetMapper.h>
21 #include <vtkLODActor.h>
22 #include <vtkPoints.h>
23 #include <vtkPolyData.h>
24 #include <vtkPolygon.h>
25 #include <vtkSmartPointer.h>
26 #include <vtkRenderer.h>
27 #include <vtkRenderWindow.h>
28 #include <vtkRenderWindowInteractor.h>
29 #include <vtkDelaunay3D.h>
30 #include <vtkVertexGlyphFilter.h>
31 #include <vtkPolyDataMapper.h>
32 #include <vtkSphereSource.h>
33 #include <vtkUnsignedCharArray.h>
34 #include <vtkPointData.h>
35 #include <vtkGlyph3D.h>
36 #include <vtkCamera.h>
37 
38 #include <vtkWindowToImageFilter.h>
39 #include <vtkPNGWriter.h>
40 
41 #include "OrientationsViewerCLP.h"
42 #include "itkSamplingScheme3D.h"
43 #include "utl.h"
44 
45 #include "utlVTKMacro.h"
46 
47 
52 int
53 main ( int argc, char *argv[] )
54 {
55 
56  PARSE_ARGS;
57 
58  int numberOfShells = _OrientationFile.size();
59  utlGlobalException(numberOfShells==0, "no input orientation file");
60  std::cout << numberOfShells << " shells" << std::endl << std::flush;
61 
62  if (!_OnlyCombined && _OpacityMesh.size()!=numberOfShells)
63  {
64  _OpacityMesh.resize(numberOfShells);
65  for ( int i = 0; i < numberOfShells; i += 1 )
66  _OpacityMesh[i]=1.0- 0.7/(double)(numberOfShells-1)*i;
67  }
68 
69  if (!_OnlyCombined && _OpacitySphere.size()!=numberOfShells)
70  {
71  _OpacitySphere.resize(numberOfShells);
72  for ( int i = 0; i < numberOfShells; i += 1 )
73  _OpacitySphere[i]=_OpacityMesh[i];
74  }
75 
76  std::vector<std::vector<double> > radiusVec;
77  if (_RadiusFileArg.isSet())
78  {
79  double maxRadius = -1;
80  utlGlobalException(_RadiusFile.size()!=numberOfShells, "wrong size of _RadiusFile");
81  for ( int i = 0; i < numberOfShells; i += 1 )
82  {
83  std::vector<double> tmp;
84  utl::ReadVector(_RadiusFile[i], tmp);
85  radiusVec.push_back(tmp);
86  double maxTmp = *std::max_element(tmp.begin(), tmp.end());
87  if (maxTmp>maxRadius)
88  maxRadius = maxTmp;
89  }
90  for ( int i = 0; i < numberOfShells; i += 1 )
91  for ( int j = 0; j < radiusVec[i].size(); j += 1 )
92  radiusVec[i][j] /= maxRadius;
93  }
94  else
95  {
96  if (_Radius.size()!=numberOfShells)
97  {
98  _Radius.resize(numberOfShells);
99  for ( int i = 0; i < numberOfShells; i += 1 )
100  _Radius[i]=i+1;
101  }
102  }
103 
104  if (_Color.size()!=3*numberOfShells)
105  {
106  _Color.resize(3*numberOfShells);
107  if (numberOfShells==1)
108  {
109  _Color[0]=1.0, _Color[1]=1.0, _Color[2]=1.0;
110  }
111  else
112  {
113  if (numberOfShells>1)
114  {
115  _Color[0]=1.0, _Color[1]=0.0, _Color[2]=0.0;
116  if (numberOfShells>=2)
117  _Color[3]=0.0, _Color[4]=1.0, _Color[5]=0.0;
118  if (numberOfShells==3)
119  _Color[6]=0.0, _Color[7]=0.0, _Color[8]=1.0;
120  }
121  if (numberOfShells>3 && numberOfShells<=6)
122  {
123  _Color[9]=0.5, _Color[10]=0.0, _Color[11]=0.0;
124  if (numberOfShells>=7)
125  _Color[12]=0.0, _Color[13]=0.5, _Color[14]=0.0;
126  if (numberOfShells==8)
127  _Color[15]=0.0, _Color[16]=0.0, _Color[17]=0.5;
128  }
129  }
130  utlGlobalException(numberOfShells>7, "too many shells, please manually set _Color");
131  }
132 
133  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
134  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
135  renderer->SetBackground(_BackgroundColor[0], _BackgroundColor[1], _BackgroundColor[2]);
136 
137  typedef itk::SamplingScheme3D<double> SamplingType;
138 
139  std::vector<vtkSmartPointer<vtkPoints> > points(numberOfShells+1);
140  std::vector<vtkSmartPointer<vtkPolyData> > polydata(numberOfShells+1);
141  std::vector<vtkSmartPointer<vtkGlyph3D> > glyph3D(numberOfShells+1);
142 
143  std::vector<vtkSmartPointer<vtkVertexGlyphFilter> > vertexGlyphFilter(numberOfShells+1);
144  std::vector<vtkSmartPointer<vtkPolyDataMapper> > pointsMapper(numberOfShells+1);
145  std::vector<vtkSmartPointer<vtkLODActor> > pointsActor(numberOfShells+1);
146 
147  std::vector<vtkSmartPointer<vtkDelaunay3D> > delaunay3D(numberOfShells+1);
148  std::vector<vtkSmartPointer<vtkDataSetMapper> > delaunayMapper(numberOfShells+1);
149  std::vector<vtkSmartPointer<vtkLODActor> > delaunayActor(numberOfShells+1);
150 
151  std::vector<vtkSmartPointer<vtkSphereSource> > sphereSource(numberOfShells+1);
152  std::vector<vtkSmartPointer<vtkPolyDataMapper> > sphereMapper(numberOfShells+1);
153  std::vector<vtkSmartPointer<vtkLODActor> > sphereActor(numberOfShells+1);
154 
155  unsigned char colorPoint[3] = {0, 0, 0};
156  vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
157  colors->SetName("colors");
158 
159  if (_OnlyCombined)
160  {
161  points.back() = vtkSmartPointer<vtkPoints>::New();
162  polydata.back() = vtkSmartPointer<vtkPolyData>::New();
163  colors->SetNumberOfComponents(3);
164  }
165 
166  vtkSmartPointer<vtkSphereSource> sphereSourceTemp = vtkSmartPointer<vtkSphereSource>::New();
167  sphereSourceTemp->SetCenter(0.0, 0.0, 0.0);
168  sphereSourceTemp->SetRadius(_SpherePointSize);
169  sphereSourceTemp->SetThetaResolution(20);
170  sphereSourceTemp->SetPhiResolution(20);
171 
172 
173  for ( int i = 0; i < _OrientationFile.size(); i += 1 )
174  {
175  SamplingType::Pointer sampling = SamplingType::New();
176 
177  if (_NoSymmetricPoints)
178  sampling->ReadOrientationFile(_OrientationFile[i], DIRECTION_NODUPLICATE);
179  else
180  sampling->ReadOrientationFile(_OrientationFile[i], DIRECTION_DUPLICATE);
181  // utlPrintVar2(true, _NoSymmetricPoints, sampling->GetNumberOfSamples());
182  std::cout << "shell " << i+1 << ", " << sampling->GetNumberOfSamples() << (_NoSymmetricPoints?" points":" antipodal symmetric points") << std::endl << std::flush;
183 
184  points[i] = vtkSmartPointer<vtkPoints>::New();
185 
186  if (_OnlyCombined)
187  {
188  colorPoint[0] = (unsigned char)(255*_Color[3*i]);
189  colorPoint[1] = (unsigned char)(255*_Color[3*i+1]);
190  colorPoint[2] = (unsigned char)(255*_Color[3*i+2]);
191  }
192  // utlPrintVar3(true, colorPoint[0], colorPoint[1], colorPoint[2]);
193  // utlPrintVar3(true, 255*_Color[3*i], 255*_Color[3*i+1], 255*_Color[3*i+2]);
194 
195  for ( unsigned int j = 0; j < sampling->GetNumberOfSamples(); j++ )
196  {
197  double x0 = (*sampling)[j][0];
198  double x1 = (*sampling)[j][1];
199  double x2 = (*sampling)[j][2];
200  // utlPrintVar3(true, x0,x1,x2);
201  if (_RadiusFileArg.isSet())
202  points[i]->InsertNextPoint( radiusVec[i][j]*x0, radiusVec[i][j]*x1, radiusVec[i][j]*x2 );
203  else
204  points[i]->InsertNextPoint( _Radius[i]*x0, _Radius[i]*x1, _Radius[i]*x2 );
205  if (_OnlyCombined)
206  {
207  points.back()->InsertNextPoint( x0, x1, x2 );
208  colors->InsertNextTupleValue(colorPoint);
209  }
210  }
211 
212  polydata[i] = vtkSmartPointer<vtkPolyData>::New();
213  polydata[i]->SetPoints(points[i]);
214  // polydata[i]->Print(std::cout<<"polydata[i]\n");
215 
216  if (!_NoPoints && !_OnlyCombined)
217  {
218  if (_PointType=="SPHERE")
219  {
220  glyph3D[i] = vtkSmartPointer<vtkGlyph3D>::New();
221  glyph3D[i]->SetSourceConnection(sphereSourceTemp->GetOutputPort());
222  vtkSetInputData(glyph3D[i], polydata[i]);
223  glyph3D[i]->Update();
224  pointsMapper[i] = vtkSmartPointer<vtkPolyDataMapper>::New();
225  vtkSetInputData(pointsMapper[i], glyph3D[i]->GetOutput());
226  }
227  else if (_PointType=="VERTEX")
228  {
229  vertexGlyphFilter[i] = vtkSmartPointer<vtkVertexGlyphFilter>::New();
230  vtkSetInputData(vertexGlyphFilter[i], polydata[i]);
231  vertexGlyphFilter[i]->Update();
232  pointsMapper[i] = vtkSmartPointer<vtkPolyDataMapper>::New();
233  vtkSetInputData(pointsMapper[i], vertexGlyphFilter[i]->GetOutput());
234  }
235 
236 
237  // unsigned char red[3] = {255, 0, 0};
238  // vtkSmartPointer<vtkUnsignedCharArray> colors = vtkSmartPointer<vtkUnsignedCharArray>::New();
239  // colors->SetNumberOfComponents(3);
240  // for ( unsigned int i = 0; i < sampling->GetNumberOfSamples(); i++ )
241  // {
242  // colors->InsertNextTupleValue(red);
243  // }
244  // polydata[i]->GetPointData()->SetScalars(colors);
245 
246 
247  pointsActor[i] = vtkSmartPointer<vtkLODActor>::New();
248  pointsActor[i]->SetMapper(pointsMapper[i]);
249  pointsActor[i]->GetProperty()->SetColor(_Color[3*i],_Color[3*i+1],_Color[3*i+2]);
250  if (_PointType=="VERTEX")
251  pointsActor[i]->GetProperty()->SetPointSize(_VertexPointSize);
252 
253  renderer->AddActor(pointsActor[i]);
254  }
255 
256  if (_Mesh && !_OnlyCombined)
257  {
258  delaunay3D[i] = vtkSmartPointer<vtkDelaunay3D>::New();
259  vtkSetInputData(delaunay3D[i], polydata[i]);
260  delaunay3D[i]->SetTolerance(1e-8);
261  delaunay3D[i]->SetOffset(10);
262  delaunay3D[i]->Update();
263  // delaunay3D[i]->Print(std::cout<<"delaunay3D[i]=\n");
264 
265  delaunayMapper[i] = vtkSmartPointer<vtkDataSetMapper>::New();
266  delaunayMapper[i]->SetInputConnection(delaunay3D[i]->GetOutputPort());
267 
268  delaunayActor[i] = vtkSmartPointer<vtkLODActor>::New();
269  delaunayActor[i]->SetMapper(delaunayMapper[i]);
270  delaunayActor[i]->GetProperty()->SetColor(1,1,1);
271  delaunayActor[i]->GetProperty()->SetOpacity(_OpacityMesh[i]);
272 
273  renderer->AddActor(delaunayActor[i]);
274  }
275 
276  if (_Sphere && !_OnlyCombined)
277  {
278  sphereSource[i] = vtkSmartPointer<vtkSphereSource>::New();
279  sphereSource[i]->SetCenter(0.0, 0.0, 0.0);
280  utlGlobalException(_RadiusFileArg.isSet(), "not support to visualize sphere if _RadiusFile is set");
281  sphereSource[i]->SetRadius(_Radius[i]);
282  sphereSource[i]->SetThetaResolution(50);
283  sphereSource[i]->SetPhiResolution(50);
284 
285  sphereMapper[i] = vtkSmartPointer<vtkPolyDataMapper>::New();
286  sphereMapper[i]->SetInputConnection(sphereSource[i]->GetOutputPort());
287 
288  sphereActor[i] = vtkSmartPointer<vtkLODActor>::New();
289  sphereActor[i]->SetMapper(sphereMapper[i]);
290  // sphereActor[i]->GetProperty()->SetColor(_Color[3*i],_Color[3*i+1],_Color[3*i+2]);
291  sphereActor[i]->GetProperty()->SetColor(1,1,1);
292  sphereActor[i]->GetProperty()->SetOpacity(_OpacitySphere[i]);
293 
294  renderer->AddActor(sphereActor[i]);
295  }
296  }
297 
298  if (_OnlyCombined)
299  {
300  polydata.back() = vtkSmartPointer<vtkPolyData>::New();
301  polydata.back()->SetPoints(points.back());
302  polydata.back()->GetPointData()->AddArray(colors);
303 
304  if (!_NoPoints)
305  {
306  vtkSmartPointer<vtkPolyData> polydata2 = vtkSmartPointer<vtkPolyData>::New();
307  if (_PointType=="SPHERE")
308  {
309  glyph3D.back() = vtkSmartPointer<vtkGlyph3D>::New();
310  glyph3D.back()->SetSourceConnection(sphereSourceTemp->GetOutputPort());
311  vtkSetInputData(glyph3D.back(), polydata.back());
312  glyph3D.back()->Update();
313  // polydata2->ShallowCopy(glyph3D.back()->GetOutput());
314 // glyph3D.back()->GetOutput()->GetPointData()->SetScalars(colors);
315  pointsMapper.back() = vtkSmartPointer<vtkPolyDataMapper>::New();
316  vtkSetInputData(pointsMapper.back(), glyph3D.back()->GetOutput());
317  }
318  else if (_PointType=="VERTEX")
319  {
320  vertexGlyphFilter.back() = vtkSmartPointer<vtkVertexGlyphFilter>::New();
321  vtkSetInputData(vertexGlyphFilter.back(), polydata.back());
322  vertexGlyphFilter.back()->Update();
323  // polydata2->ShallowCopy(vertexGlyphFilter.back()->GetOutput());
324 // vertexGlyphFilter.back()->GetOutput()->GetPointData()->SetScalars(colors);
325  pointsMapper.back() = vtkSmartPointer<vtkPolyDataMapper>::New();
326  vtkSetInputData(pointsMapper.back(), vertexGlyphFilter.back()->GetOutput());
327  }
328 
329  pointsMapper.back()->ScalarVisibilityOn();
330  pointsMapper.back()->SetScalarModeToUsePointFieldData();
331  pointsMapper.back()->SetColorModeToDefault();
332  pointsMapper.back()->SelectColorArray("colors");
333 
334  // polydata2->GetPointData()->SetScalars(colors);
335 
336  // pointsMapper.back() = vtkSmartPointer<vtkPolyDataMapper>::New();
337  // vtkSetInputData(pointsMapper.back(), polydata2);
338 
339  pointsActor.back() = vtkSmartPointer<vtkLODActor>::New();
340  pointsActor.back()->SetMapper(pointsMapper.back());
341  if (_PointType=="VERTEX")
342  pointsActor.back()->GetProperty()->SetPointSize(_VertexPointSize);
343 
344  renderer->AddActor(pointsActor.back());
345  }
346 
347  if (_Mesh)
348  {
349  delaunay3D.back() = vtkSmartPointer<vtkDelaunay3D>::New();
350  vtkSetInputData(delaunay3D.back(), polydata.back());
351  delaunay3D.back()->SetTolerance(1e-8);
352  delaunay3D.back()->SetOffset(10);
353  delaunay3D.back()->Update();
354  // delaunay3D.back()->Print(std::cout<<"delaunay3D.back()=\n");
355 
356  delaunayMapper.back() = vtkSmartPointer<vtkDataSetMapper>::New();
357  delaunayMapper.back()->SetInputConnection(delaunay3D.back()->GetOutputPort());
358 
359  delaunayActor.back() = vtkSmartPointer<vtkLODActor>::New();
360  delaunayActor.back()->SetMapper(delaunayMapper.back());
361  delaunayActor.back()->GetProperty()->SetColor(1,1,1);
362  delaunayActor.back()->GetProperty()->SetOpacity(_OpacityMesh[0]);
363 
364  renderer->AddActor(delaunayActor.back());
365  }
366 
367  if (_Sphere)
368  {
369  sphereSource.back() = vtkSmartPointer<vtkSphereSource>::New();
370  sphereSource.back()->SetCenter(0.0, 0.0, 0.0);
371  sphereSource.back()->SetRadius(1.0);
372  sphereSource.back()->SetThetaResolution(50);
373  sphereSource.back()->SetPhiResolution(50);
374 
375  sphereMapper.back() = vtkSmartPointer<vtkPolyDataMapper>::New();
376  sphereMapper.back()->SetInputConnection(sphereSource.back()->GetOutputPort());
377 
378  sphereActor.back() = vtkSmartPointer<vtkLODActor>::New();
379  sphereActor.back()->SetMapper(sphereMapper.back());
380  // sphereActor.back()->GetProperty()->SetColor(_Color[3*i],_Color[3*i+1],_Color[3*i+2]);
381 // sphereActor.back()->GetProperty()->SetColor(1,1,1);
382 // sphereActor.back()->GetProperty()->SetOpacity(_OpacitySphere[0]);
383 
384  renderer->AddActor(sphereActor.back());
385  }
386  }
387 
388 
389  renderWindow->AddRenderer(renderer);
390  renderWindow->SetSize(_Window[0],_Window[1]);
391 
392  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
393  renderWindowInteractor->SetRenderWindow(renderWindow);
394  renderWindowInteractor->SetDesiredUpdateRate(25);
395 
396  renderWindow->Render();
397 
398  renderer->GetActiveCamera()->Roll(_Angle[0]);
399  renderer->GetActiveCamera()->Elevation(_Angle[1]);
400 
401  renderWindow->Render();
402 
403  if (_PNGFile!="")
404  {
405  vtkSmartPointer<vtkWindowToImageFilter> windowToImage = vtkSmartPointer<vtkWindowToImageFilter>::New();
406  windowToImage->SetInput(renderWindow);
407  // windowToImage->SetMagnification(2);
408  vtkSmartPointer<vtkPNGWriter> pngWriter = vtkSmartPointer<vtkPNGWriter>::New();
409  pngWriter->SetInputConnection(windowToImage->GetOutputPort());
410  pngWriter->SetFileName(_PNGFile.c_str());
411  pngWriter->Write();
412  }
413  else
414  renderWindowInteractor->Start();
415 
416 
417 
418  return EXIT_SUCCESS;
419 }
420 
helper functions specifically used in dmritool
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
Created "02-18-2015.
this class describes sampling in a 3D space (Q space or R space).
#define vtkSetInputData(x, y)
Definition: utlVTKMacro.h:18
int main(int argc, char *argv[])
orientations visualization
void ReadVector(const std::string &vectorStr, std::vector< T > &vec, const char *cc=" ")
Definition: utlCore.h:1159
T max_element(const std::vector< T > &v)
Definition: utlCore.h:205