PlusLib  2.9.0
Software library for tracked ultrasound image acquisition, calibration, and processing.
vtkUsSimulatorAlgoTest.cxx
Go to the documentation of this file.
1 /*=Plus=header=begin======================================================
2 Program: Plus
3 Copyright (c) Laboratory for Percutaneous Surgery. All rights reserved.
4 See License.txt for details.
5 =========================================================Plus=header=end*/
6 
7 #include "PlusConfigure.h"
8 #include "igsioMath.h"
9 #include "igsioTrackedFrame.h"
10 #include "vtkAppendPolyData.h"
11 #include "vtkCubeSource.h"
12 #include "vtkImageData.h"
13 #include "vtkMatrix4x4.h"
14 #include "vtkPointData.h"
15 #include "vtkSTLWriter.h"
16 #include "vtkPlusSequenceIO.h"
17 #include "vtkSmartPointer.h"
18 #include "vtkTimerLog.h"
19 #include "vtkIGSIOTrackedFrameList.h"
20 #include "vtkTransform.h"
21 #include "vtkTransformPolyDataFilter.h"
22 #include "vtkIGSIOTransformRepository.h"
23 #include "vtkPlusUsSimulatorAlgo.h"
24 #include "vtkXMLImageDataWriter.h"
25 #include "vtkXMLUtilities.h"
26 #include "vtksys/CommandLineArguments.hxx"
27 #include <iomanip>
28 #include <iostream>
29 
30 //display
31 #include "vtkImageActor.h"
32 #include "vtkActor.h"
33 #include "vtkRenderWindow.h"
34 #include "vtkRenderer.h"
35 #include "vtkRenderWindowInteractor.h"
36 #include "vtkPolyDataMapper.h"
37 #include "vtkProperty.h"
38 #include "vtkInteractorStyleImage.h"
39 
40 //-----------------------------------------------------------------------------
41 void CreateSliceModels(vtkIGSIOTrackedFrameList* trackedFrameList, vtkIGSIOTransformRepository* transformRepository, igsioTransformName& imageToReferenceTransformName, vtkPolyData* outputPolyData)
42 {
43  // Prepare the output polydata.
44  vtkSmartPointer< vtkAppendPolyData > appender = vtkSmartPointer<vtkAppendPolyData>::New();
45 
46  // Loop over each tracked image slice.
47  for (unsigned int frameIndex = 0; frameIndex < trackedFrameList->GetNumberOfTrackedFrames(); ++ frameIndex)
48  {
49  igsioTrackedFrame* frame = trackedFrameList->GetTrackedFrame(frameIndex);
50 
51  // Update transform repository
52  if (transformRepository->SetTransforms(*frame) != PLUS_SUCCESS)
53  {
54  LOG_ERROR("Failed to set repository transforms from tracked frame!");
55  continue;
56  }
57 
58  vtkSmartPointer<vtkMatrix4x4> tUserDefinedMatrix = vtkSmartPointer<vtkMatrix4x4>::New();
59  if (transformRepository->GetTransform(imageToReferenceTransformName, tUserDefinedMatrix) != PLUS_SUCCESS)
60  {
61  std::string strTransformName;
62  imageToReferenceTransformName.GetTransformName(strTransformName);
63  LOG_ERROR("Failed to get transform from repository: " << strTransformName);
64  continue;
65  }
66 
67  vtkSmartPointer< vtkTransform> tUserDefinedTransform = vtkSmartPointer< vtkTransform >::New();
68  tUserDefinedTransform->SetMatrix(tUserDefinedMatrix);
69 
70  FrameSizeType frameSize = frame->GetFrameSize();
71 
72  vtkSmartPointer<vtkTransform> tCubeToImage = vtkSmartPointer<vtkTransform>::New();
73  tCubeToImage->Scale(frameSize[ 0 ], frameSize[ 1 ], 1);
74  tCubeToImage->Translate(0.5, 0.5, 0.5); // Moving the corner to the origin.
75 
76  vtkSmartPointer<vtkTransform> tCubeToTracker = vtkSmartPointer< vtkTransform >::New();
77  tCubeToTracker->Identity();
78  tCubeToTracker->Concatenate(tUserDefinedTransform);
79  tCubeToTracker->Concatenate(tCubeToImage);
80 
81  vtkSmartPointer<vtkTransformPolyDataFilter > cubeToTracker = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
82  cubeToTracker->SetTransform(tCubeToTracker);
83  vtkSmartPointer<vtkCubeSource> source = vtkSmartPointer<vtkCubeSource>::New();
84  cubeToTracker->SetInputConnection(source->GetOutputPort());
85  cubeToTracker->Update();
86 
87  appender->AddInputConnection(cubeToTracker->GetOutputPort());
88  }
89 
90  appender->Update();
91  outputPolyData->DeepCopy(appender->GetOutput());
92 }
93 
94 //-----------------------------------------------------------------------------
95 void ShowResults(vtkIGSIOTrackedFrameList* trackedFrameList, vtkIGSIOTransformRepository* transformRepository, igsioTransformName imageToReferenceTransformName, std::string intersectionFile)
96 {
97  // Setup Renderer to visualize surface model and ultrasound planes
98  vtkSmartPointer<vtkRenderer> rendererPoly = vtkSmartPointer<vtkRenderer>::New();
99  vtkSmartPointer<vtkRenderWindow> renderWindowPoly = vtkSmartPointer<vtkRenderWindow>::New();
100  renderWindowPoly->AddRenderer(rendererPoly);
101  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractorPoly = vtkSmartPointer<vtkRenderWindowInteractor>::New();
102  renderWindowInteractorPoly->SetRenderWindow(renderWindowPoly);
103  vtkSmartPointer<vtkInteractorStyleTrackballCamera> style = vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
104  renderWindowInteractorPoly->SetInteractorStyle(style);
105 
106  // Visualization of the surface model
107 
108  /*
109  TODO: add surface model of each SpatialModel in the simulator
110  vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
111  mapper->SetInputData((vtkPolyData*)usSimulator->GetInput());
112  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
113  actor->SetMapper(mapper);
114  rendererPoly->AddActor(actor);
115  */
116 
117  // Visualization of the image planes
118  vtkSmartPointer< vtkPolyData > slicesPolyData = vtkSmartPointer< vtkPolyData >::New();
119  CreateSliceModels(trackedFrameList, transformRepository, imageToReferenceTransformName, slicesPolyData);
120 
121  if (!intersectionFile.empty())
122  {
123  vtkSmartPointer<vtkSTLWriter> surfaceModelWriter = vtkSmartPointer<vtkSTLWriter>::New();
124  surfaceModelWriter->SetFileName(intersectionFile.c_str());
125  surfaceModelWriter->SetInputData(slicesPolyData);
126  surfaceModelWriter->Write();
127  }
128 
129  vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
130  mapper->SetInputData(slicesPolyData);
131  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
132  actor->SetMapper(mapper);
133  rendererPoly->AddActor(actor);
134 
135  renderWindowPoly->Render();
136  renderWindowInteractorPoly->Start();
137 }
138 
139 //-----------------------------------------------------------------------------
140 void ShowImage(vtkImageData* simOutput)
141 {
142  vtkSmartPointer<vtkImageData> usImage = vtkSmartPointer<vtkImageData>::New();
143  usImage->DeepCopy(simOutput);
144 
145  // Display output of filter
146  vtkSmartPointer<vtkImageActor> redImageActor = vtkSmartPointer<vtkImageActor>::New();
147  redImageActor->SetInputData(simOutput);
148 
149  // Visualize
150  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
151 
152  // Red image is displayed
153  renderer->AddActor(redImageActor);
154  renderer->ResetCamera();
155 
156  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
157  renderWindow->AddRenderer(renderer);
158  renderer->SetBackground(0, 72, 0);
159 
160  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
161  vtkSmartPointer<vtkInteractorStyleImage> style = vtkSmartPointer<vtkInteractorStyleImage>::New();
162 
163  renderWindowInteractor->SetInteractorStyle(style);
164  renderWindowInteractor->SetRenderWindow(renderWindow);
165 }
166 
167 //-----------------------------------------------------------------------------
168 int main(int argc, char** argv)
169 {
170  bool printHelp = false;
171  std::string inputTransformsFile;
172  std::string inputConfigFileName;
173  std::string outputUsImageFile;
174  std::string intersectionFile;
175  bool showResults = false;
176  bool useCompression(true);
177 
178  int verboseLevel = vtkPlusLogger::LOG_LEVEL_UNDEFINED;
179 
180  vtksys::CommandLineArguments args;
181  args.Initialize(argc, argv);
182 
183  args.AddArgument("--help", vtksys::CommandLineArguments::NO_ARGUMENT, &printHelp, "Print this help.");
184  args.AddArgument("--verbose", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &verboseLevel, "Verbose level (1=error only, 2=warning, 3=info, 4=debug, 5=trace)");
185  args.AddArgument("--config-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &inputConfigFileName, "Config file containing the image to probe and phantom to reference transformations");
186  args.AddArgument("--transforms-seq-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &inputTransformsFile, "Input file containing coordinate frames and the associated model to image transformations");
187  args.AddArgument("--use-compression", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &useCompression, "Use compression when outputting data");
188  args.AddArgument("--output-us-img-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &outputUsImageFile, "File name of the generated output ultrasound image");
189  args.AddArgument("--output-slice-model-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &intersectionFile, "Name of STL output file containing the model of all the frames (optional)");
190  args.AddArgument("--show-results", vtksys::CommandLineArguments::NO_ARGUMENT, &showResults, "Show the simulated image on the screen");
191 
192  // Input arguments error checking
193  if (!args.Parse())
194  {
195  std::cerr << "Problem parsing arguments" << std::endl;
196  std::cout << "Help: " << args.GetHelp() << std::endl;
197  exit(EXIT_FAILURE);
198  }
199 
200  if (printHelp)
201  {
202  std::cout << args.GetHelp() << std::endl;
203  exit(EXIT_SUCCESS);
204  }
205 
206  vtkPlusLogger::Instance()->SetLogLevel(verboseLevel);
207 
208  if (inputConfigFileName.empty())
209  {
210  std::cerr << "--config-file required " << std::endl;
211  exit(EXIT_FAILURE);
212  }
213  if (inputTransformsFile.empty())
214  {
215  std::cerr << "--transforms-seq-file required" << std::endl;
216  exit(EXIT_FAILURE);
217  }
218  if (outputUsImageFile.empty())
219  {
220  std::cerr << "--output-us-img-file required" << std::endl;
221  exit(EXIT_FAILURE);
222  }
223 
224  // Read transformations data
225  LOG_DEBUG("Reading input meta file...");
226  vtkSmartPointer< vtkIGSIOTrackedFrameList > trackedFrameList = vtkSmartPointer< vtkIGSIOTrackedFrameList >::New();
227  if (vtkPlusSequenceIO::Read(inputTransformsFile, trackedFrameList) != PLUS_SUCCESS)
228  {
229  LOG_ERROR("Unable to load input sequences file.");
230  exit(EXIT_FAILURE);
231  }
232  LOG_DEBUG("Reading input meta file completed");
233 
234  // Create repository for ultrasound images correlated to the iput tracked frames
235  vtkSmartPointer<vtkIGSIOTrackedFrameList> simulatedUltrasoundFrameList = vtkSmartPointer<vtkIGSIOTrackedFrameList>::New();
236 
237  // Read config file
238  LOG_DEBUG("Reading config file...")
239  vtkSmartPointer<vtkXMLDataElement> configRootElement = vtkSmartPointer<vtkXMLDataElement>::New();
240  if (PlusXmlUtils::ReadDeviceSetConfigurationFromFile(configRootElement, inputConfigFileName.c_str()) == PLUS_FAIL)
241  {
242  LOG_ERROR("Unable to read configuration from file " << inputConfigFileName.c_str());
243  return EXIT_FAILURE;
244  }
245  LOG_DEBUG("Reading config file finished.");
246 
247  // Create transform repository
248  vtkSmartPointer<vtkIGSIOTransformRepository> transformRepository = vtkSmartPointer<vtkIGSIOTransformRepository>::New();
249  if (transformRepository->ReadConfiguration(configRootElement) != PLUS_SUCCESS)
250  {
251  LOG_ERROR("Failed to read transforms for transform repository!");
252  exit(EXIT_FAILURE);
253  }
254 
255  // Create simulator
256  vtkSmartPointer<vtkPlusUsSimulatorAlgo> usSimulator = vtkSmartPointer<vtkPlusUsSimulatorAlgo>::New();
257  if (usSimulator->ReadConfiguration(configRootElement) != PLUS_SUCCESS)
258  {
259  LOG_ERROR("Failed to read US simulator configuration!");
260  exit(EXIT_FAILURE);
261  }
262  usSimulator->SetTransformRepository(transformRepository);
263  igsioTransformName imageToReferenceTransformName(usSimulator->GetImageCoordinateFrame(), usSimulator->GetReferenceCoordinateFrame());
264 
265  // Write slice model file
266  if (!intersectionFile.empty())
267  {
268  vtkSmartPointer< vtkPolyData > slicesPolyData = vtkSmartPointer< vtkPolyData >::New();
269  CreateSliceModels(trackedFrameList, transformRepository, imageToReferenceTransformName, slicesPolyData);
270  vtkSmartPointer<vtkSTLWriter> surfaceModelWriter = vtkSmartPointer<vtkSTLWriter>::New();
271  surfaceModelWriter->SetFileName(intersectionFile.c_str());
272  surfaceModelWriter->SetInputData(slicesPolyData);
273  surfaceModelWriter->Write();
274  }
275 
276  if (showResults)
277  {
278  ShowResults(trackedFrameList, transformRepository, imageToReferenceTransformName, intersectionFile);
279  }
280 
281  std::vector<double> timeElapsedPerFrameSec;
282  double startTimeSec = 0;
283  double endTimeSec = 0;
284 
285  // TODO: remove this, it's just for testing
286  trackedFrameList->RemoveTrackedFrameRange(15, trackedFrameList->GetNumberOfTrackedFrames() - 1);
287 
288  for (unsigned int i = 0; i < trackedFrameList->GetNumberOfTrackedFrames(); i++)
289  {
290  startTimeSec = vtkTimerLog::GetUniversalTime();
291 
292  LOG_DEBUG("Processing frame " << i);
293  igsioTrackedFrame* frame = trackedFrameList->GetTrackedFrame(i);
294 
295  // Update transform repository
296  if (transformRepository->SetTransforms(*frame) != PLUS_SUCCESS)
297  {
298  LOG_ERROR("Failed to set repository transforms from tracked frame!");
299  return EXIT_FAILURE;
300  }
301 
302  // TODO: would be nice to implement an automatic mechanism to find out if the transforms have been changed
303  usSimulator->Modified(); // Signal that the transforms have changed so we need to recompute
304  usSimulator->Update();
305  vtkImageData* simOutput = usSimulator->GetOutput();
306 
307  igsioVideoFrame* simulatorOutputigsioVideoFrame = new igsioVideoFrame();
308  simulatorOutputigsioVideoFrame->DeepCopyFrom(simOutput);
309 
310  frame->SetImageData(*simulatorOutputigsioVideoFrame);
311 
312  /*
313  double origin[3]=
314  {
315  imageToReferenceMatrix->Element[0][3],
316  imageToReferenceMatrix->Element[1][3],
317  imageToReferenceMatrix->Element[2][3]
318  };
319  simulatedUsImage->SetOrigin(origin);
320 
321  double spacing[3]=
322  {
323  sqrt(imageToReferenceMatrix->Element[0][0]*imageToReferenceMatrix->Element[0][0]+
324  imageToReferenceMatrix->Element[1][0]*imageToReferenceMatrix->Element[1][0]+
325  imageToReferenceMatrix->Element[2][0]*imageToReferenceMatrix->Element[2][0]),
326  sqrt(imageToReferenceMatrix->Element[0][1]*imageToReferenceMatrix->Element[0][1]+
327  imageToReferenceMatrix->Element[1][1]*imageToReferenceMatrix->Element[1][1]+
328  imageToReferenceMatrix->Element[2][1]*imageToReferenceMatrix->Element[2][1]),
329  1.0
330  };
331  simulatedUsImage->SetSpacing(spacing);
332  */
333 
334  if (showResults)
335  {
336  ShowImage(simOutput);
337  }
338 
339  endTimeSec = vtkTimerLog::GetUniversalTime();
340  timeElapsedPerFrameSec.push_back(endTimeSec - startTimeSec);
341  }
342 
343  if (vtkPlusSequenceIO::Write(outputUsImageFile, trackedFrameList, trackedFrameList->GetImageOrientation(), useCompression) != PLUS_SUCCESS)
344  {
345  // Error has already been logged
346  return EXIT_FAILURE;
347  }
348 
349  LOG_INFO("Computation time for the fits frame (not included in the statistics because of BSP Tree Building, in sec): " << timeElapsedPerFrameSec.at(0));
350 
351  // Remove the first frame from the statistics computation because extra processing is done for the first frame
352  timeElapsedPerFrameSec.erase(timeElapsedPerFrameSec.begin());
353 
354  double meanTimeElapsedPerFrameSec = 0;
355  double stdevTimeElapsedPerFrameSec = 0;
356  igsioMath::ComputeMeanAndStdev(timeElapsedPerFrameSec, meanTimeElapsedPerFrameSec, stdevTimeElapsedPerFrameSec);
357  LOG_INFO(" Average computation time per frame (sec): " << meanTimeElapsedPerFrameSec) ;
358  LOG_INFO(" Standard dev computation time per frame (sec): " << stdevTimeElapsedPerFrameSec) ;
359  LOG_INFO(" Average fps: " << 1 / meanTimeElapsedPerFrameSec) ;
360 
361  return EXIT_SUCCESS;
362 }
const char * source
Definition: phidget22.h:2461
void ShowImage(vtkImageData *simOutput)
void CreateSliceModels(vtkIGSIOTrackedFrameList *trackedFrameList, vtkIGSIOTransformRepository *transformRepository, igsioTransformName &imageToReferenceTransformName, vtkPolyData *outputPolyData)
for i
#define PLUS_FAIL
Definition: PlusCommon.h:43
static igsioStatus Write(const std::string &filename, igsioTrackedFrame *frame, US_IMAGE_ORIENTATION orientationInFile=US_IMG_ORIENT_MF, bool useCompression=true, bool EnableImageDataWrite=true)
#define PLUS_SUCCESS
Definition: PlusCommon.h:44
static igsioStatus Read(const std::string &filename, vtkIGSIOTrackedFrameList *frameList)
int main(int argc, char **argv)
void ShowResults(vtkIGSIOTrackedFrameList *trackedFrameList, vtkIGSIOTransformRepository *transformRepository, igsioTransformName imageToReferenceTransformName, std::string intersectionFile)
static vtkIGSIOLogger * Instance()
static PlusStatus ReadDeviceSetConfigurationFromFile(vtkXMLDataElement *config, const char *filename)
Definition: PlusXmlUtils.h:23