PlusLib  2.9.0
Software library for tracked ultrasound image acquisition, calibration, and processing.
vtkFreehandCalibrationStatisticalEvaluation.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 
13 #include "PlusConfigure.h"
14 
16 #include "PlusMath.h"
17 #include "igsioTrackedFrame.h"
18 #include "vtkCallbackCommand.h"
19 #include "vtkCommand.h"
20 #include "vtkMath.h"
21 #include "vtkMatrix4x4.h"
23 #include "vtkIGSIOSequenceIO.h"
24 #include "vtkSmartPointer.h"
25 #include "vtkIGSIOTrackedFrameList.h"
26 #include "vtkTransform.h"
27 #include "vtkIGSIOTransformRepository.h"
28 #include "vtkIGSIOTransformRepository.h"
29 #include "vtkXMLDataElement.h"
30 #include "vtkXMLUtilities.h"
31 #include "vtksys/CommandLineArguments.hxx"
32 #include "vtksys/SystemTools.hxx"
33 #include <iostream>
34 #include <stdlib.h>
35 
36 PlusStatus SubSequenceMetafile(vtkIGSIOTrackedFrameList* aTrackedFrameList, std::vector<unsigned int> selectedFrames);
37 PlusStatus SetOptimizationMethod(vtkPlusProbeCalibrationAlgo* freehandCalibration, std::string method);
38 
40 {
46 };
47 
48 int main(int argc, char* argv[])
49 {
50  std::string inputCalibrationSeqMetafile;
51  std::string inputValidationSeqMetafile;
52 
53  std::string inputConfigFileName;
54  std::string resultConfigFileName = "";
55  std::string saveResultsFilename;
56  std::string strOperation;
57 
58 #ifndef _WIN32
59  //double inputTranslationErrorThreshold(LINUXTOLERANCE);
60  //double inputRotationErrorThreshold(LINUXTOLERANCE);
61 #else
62  //double inputTranslationErrorThreshold(0);
63  //double inputRotationErrorThreshold(0);
64 #endif
65 
66  int verboseLevel = vtkPlusLogger::LOG_LEVEL_UNDEFINED;
67 
68  vtksys::CommandLineArguments cmdargs;
69  cmdargs.Initialize(argc, argv);
70 
71  cmdargs.AddArgument("--calibration-seq-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &inputCalibrationSeqMetafile, "Sequence metafile name of input calibration dataset.");
72  cmdargs.AddArgument("--validation-seq-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &inputValidationSeqMetafile, "Sequence metafile name of input validation dataset.");
73 
74  cmdargs.AddArgument("--config-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &inputConfigFileName, "Configuration file name prefix");
75  cmdargs.AddArgument("--save-results-file", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &saveResultsFilename, "Save results file name");
76  cmdargs.AddArgument("--operation", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &strOperation, "Type of experiment to perform");
77 
78  cmdargs.AddArgument("--verbose", vtksys::CommandLineArguments::EQUAL_ARGUMENT, &verboseLevel, "Verbose level (1=error only, 2=warning, 3=info, 4=debug, 5=trace)");
79 
80  if (!cmdargs.Parse())
81  {
82  std::cerr << "Problem parsing arguments" << std::endl;
83  std::cout << "Help: " << cmdargs.GetHelp() << std::endl;
84  exit(EXIT_FAILURE);
85  }
86 
87  vtkPlusLogger::Instance()->SetLogLevel(verboseLevel);
88 
89  LOG_INFO("Initialize");
90 
91  // Read configuration
92  vtkSmartPointer<vtkXMLDataElement> configRootElement = vtkSmartPointer<vtkXMLDataElement>::New();
93  if (PlusXmlUtils::ReadDeviceSetConfigurationFromFile(configRootElement, inputConfigFileName.c_str()) == PLUS_FAIL)
94  {
95  LOG_ERROR("Unable to read configuration from file " << inputConfigFileName.c_str());
96  return EXIT_FAILURE;
97  }
99 
100  // Read coordinate definitions
101  vtkSmartPointer<vtkIGSIOTransformRepository> transformRepository = vtkSmartPointer<vtkIGSIOTransformRepository>::New();
102  if (transformRepository->ReadConfiguration(configRootElement) != PLUS_SUCCESS)
103  {
104  LOG_ERROR("Failed to read CoordinateDefinitions!");
105  return EXIT_FAILURE;
106  }
107 
108  vtkSmartPointer<vtkPlusProbeCalibrationAlgo> freehandCalibration = vtkSmartPointer<vtkPlusProbeCalibrationAlgo>::New();
109  freehandCalibration->ReadConfiguration(configRootElement);
110 
111  PlusFidPatternRecognition patternRecognition;
113  patternRecognition.ReadConfiguration(configRootElement);
114 
115  bool debugOutput = vtkPlusLogger::Instance()->GetLogLevel() >= vtkPlusLogger::LOG_LEVEL_TRACE;
116  patternRecognition.GetFidSegmentation()->SetDebugOutput(debugOutput);
117 
118  // Load and segment calibration image
119  vtkSmartPointer<vtkIGSIOTrackedFrameList> calibrationTrackedFrameList = vtkSmartPointer<vtkIGSIOTrackedFrameList>::New();
120  if (vtkIGSIOSequenceIO::Read(inputCalibrationSeqMetafile, calibrationTrackedFrameList) != PLUS_SUCCESS)
121  {
122  LOG_ERROR("Reading calibration images from '" << inputCalibrationSeqMetafile << "' failed!");
123  return EXIT_FAILURE;
124  }
125 
126  int numberOfSuccessfullySegmentedCalibrationImages = 0;
127  std::vector<unsigned int> segmentedCalibrationFramesIndices;
128  if (patternRecognition.RecognizePattern(calibrationTrackedFrameList, error, &numberOfSuccessfullySegmentedCalibrationImages, &segmentedCalibrationFramesIndices) != PLUS_SUCCESS)
129  {
130  LOG_ERROR("Error occured during segmentation of calibration images!");
131  return EXIT_FAILURE;
132  }
133 
134  LOG_INFO("Segmentation success rate of calibration images: " << numberOfSuccessfullySegmentedCalibrationImages << " out of " << calibrationTrackedFrameList->GetNumberOfTrackedFrames());
135 
136  // Load and segment validation image
137  vtkSmartPointer<vtkIGSIOTrackedFrameList> validationTrackedFrameList = vtkSmartPointer<vtkIGSIOTrackedFrameList>::New();
138  if (vtkIGSIOSequenceIO::Read(inputValidationSeqMetafile, validationTrackedFrameList) != PLUS_SUCCESS)
139  {
140  LOG_ERROR("Reading validation images from '" << inputValidationSeqMetafile << "' failed!");
141  return EXIT_FAILURE;
142  }
143 
144  int numberOfSuccessfullySegmentedValidationImages = 0;
145  if (patternRecognition.RecognizePattern(validationTrackedFrameList, error, &numberOfSuccessfullySegmentedValidationImages) != PLUS_SUCCESS)
146  {
147  LOG_ERROR("Error occured during segmentation of validation images!");
148  return EXIT_FAILURE;
149  }
150 
151  LOG_INFO("Segmentation success rate of validation images: " << numberOfSuccessfullySegmentedValidationImages << " out of " << validationTrackedFrameList->GetNumberOfTrackedFrames());
152 
153  // Keep only the images properly segmented
154  SubSequenceMetafile(calibrationTrackedFrameList, segmentedCalibrationFramesIndices);
155 
156 
157  std::ofstream outputFile;
158  outputFile.open(saveResultsFilename.c_str(), std::ios_base::app);
159 
160 
161  OperationType operation(NO_OPERATION);
162  // Set operation
163  if (strOperation.empty())
164  {
165  LOG_INFO("No modification operation has been specified (specify --operation parameter to change the input sequence).");
166  }
167  else if (STRCASECMP(strOperation.c_str(), "INCREMENTAL_FRAME_DISTANCE") == 0)
168  {
169  operation = INCREMENTAL_FRAME_DISTANCE;
170  }
171  else if (STRCASECMP(strOperation.c_str(), "INCREMENTAL_NUMBER_OF_FRAMES") == 0)
172  {
173  operation = INCREMENTAL_NUMBER_OF_FRAMES;
174  }
175  else if (STRCASECMP(strOperation.c_str(), "SAME_NUMBER_OF_FRAMES") == 0)
176  {
177  operation = SAME_NUMBER_OF_FRAMES;
178  }
179  else if (STRCASECMP(strOperation.c_str(), "MOBILE_WINDOW") == 0)
180  {
181  operation = MOBILE_WINDOW;
182  }
183 
184  int minFrame = 0;
185  int maxFrame = segmentedCalibrationFramesIndices.size();
186  int numberOfConfigurations = 15;
187  outputFile << "Number of configurations = " << numberOfConfigurations << "\n";
188  std::vector<int> numberOfCalibrationFrames;
189  std::vector< std::vector<int> > selectedFrames;
190 
191  if (operation == INCREMENTAL_NUMBER_OF_FRAMES)
192  {
193  for (int i = 0; i < numberOfConfigurations; i++)
194  {
195  numberOfCalibrationFrames.push_back(5 + i * 5);
196  std::vector<int> frames;
197  for (int j = 0; j < numberOfCalibrationFrames.at(i) ; j++)
198  {
199  int randomIndex = rand() % (maxFrame - minFrame) + minFrame;
200  while (std::find(frames.begin(), frames.end(), randomIndex) != frames.end())
201  {
202  /* frames contains randomIndex */
203  randomIndex = rand() % (maxFrame - minFrame) + minFrame;
204  }
205  frames.push_back(randomIndex);
206  }
207  selectedFrames.push_back(frames);
208  }
209  }
210 
211 
212  if (operation == SAME_NUMBER_OF_FRAMES)
213  {
214  int numberOfFrames = 80;
215  for (int i = 0; i < numberOfConfigurations; i++)
216  {
217  numberOfCalibrationFrames.push_back(numberOfFrames);
218  std::vector<int> frames;
219  for (int j = 0; j < numberOfCalibrationFrames.at(i) ; j++)
220  {
221  int randomIndex = rand() % (maxFrame - minFrame) + minFrame;
222  while (std::find(frames.begin(), frames.end(), randomIndex) != frames.end())
223  {
224  /* frames contains randomIndex */
225  randomIndex = rand() % (maxFrame - minFrame) + minFrame;
226  }
227  frames.push_back(randomIndex);
228  }
229  selectedFrames.push_back(frames);
230  }
231  }
232 
233  if (operation == INCREMENTAL_FRAME_DISTANCE)
234  {
235  int numberOfFrames = 10;
236  for (int i = 0; i < numberOfConfigurations; i++)
237  {
238  maxFrame = (2 + i) * numberOfFrames;
239  std::vector<int> frames;
240  for (int j = 0; j < numberOfFrames ; j++)
241  {
242  int randomIndex = rand() % (maxFrame - minFrame) + minFrame;
243  while (std::find(frames.begin(), frames.end(), randomIndex) != frames.end())
244  {
245  /* frames contains randomIndex */
246  randomIndex = rand() % (maxFrame - minFrame) + minFrame;
247  }
248  frames.push_back(randomIndex);
249  }
250  selectedFrames.push_back(frames);
251  }
252  }
253 
254  if (operation == MOBILE_WINDOW)
255  {
256  int width = 40;
257  int numberOfFrames = 2 * width + 1;
258  for (int i = 0; i < numberOfConfigurations; i++)
259  {
260  int center = rand() % (maxFrame - minFrame - numberOfFrames) + width ;
261  std::vector<int> frames;
262  for (int j = 0; j < numberOfFrames ; j++)
263  {
264  frames.push_back(center - width + j);
265  }
266  selectedFrames.push_back(frames);
267  }
268  }
269 
270  vtkSmartPointer<vtkIGSIOTrackedFrameList> sequenceTrackedFrameList = vtkSmartPointer<vtkIGSIOTrackedFrameList>::New();
271 
272  std::string methods[] = {"NO_OPTIMIZATION", "7param2D", "7param3D", "8param2D", "8param3D"};
273  int numberOfMethods = sizeof(methods) / sizeof(methods[0]);
274  outputFile << "Number of methods = " << numberOfMethods << "\n";
275  std::vector<double> calibError;
276  std::vector<double> validError;
277  vnl_matrix_fixed<double, 4, 4> imageToProbeTransformMatrix;
278  FrameSizeType frameSize = calibrationTrackedFrameList->GetTrackedFrame(0)->GetFrameSize();
279  outputFile << "Frame size = " << frameSize[0] << " " << frameSize[1] << "\n";
280  for (int i = 0; i < numberOfConfigurations; i++)
281  {
282  bool calibrationFail = false;
283  int numberOfFramesOfTheSequence = selectedFrames.at(i).size();
284  outputFile << "Frames of the sequence = ";
285  for (int j = 0; j < numberOfFramesOfTheSequence; j++)
286  {
287  int indexInCalibrationSequence = selectedFrames.at(i).at(j);
288  sequenceTrackedFrameList->AddTrackedFrame(calibrationTrackedFrameList->GetTrackedFrame(indexInCalibrationSequence));
289  outputFile << indexInCalibrationSequence << " ";
290  }
291  outputFile << " \n ";
292 
293  for (int k = 0; k < numberOfMethods; k++)
294  {
295  outputFile << "Method = " << methods[k] << " \n ";
296  if (!calibrationFail)
297  {
298  SetOptimizationMethod(freehandCalibration, methods[k]);
299 
300  // Calibrate
301  if (freehandCalibration->Calibrate(validationTrackedFrameList, sequenceTrackedFrameList, transformRepository, patternRecognition.GetFidLineFinder()->GetNWires()) != PLUS_SUCCESS)
302  {
303  LOG_ERROR("Calibration failed!");
304  calibError.push_back(-1);
305  calibError.push_back(-1);
306  validError.push_back(-1);
307  validError.push_back(-1);
308  calibrationFail = true;
309  //return EXIT_FAILURE;
310  }
311 
312  freehandCalibration->GetCalibrationReport(&calibError, &validError, &imageToProbeTransformMatrix);
313  // TODO: double-check if the reported values matches the expected values
314 
315  outputFile << "Calibration error = ";
316  outputFile << calibError.at(0) << " " << calibError.at(1) << " ";
317  outputFile << validError.at(0) << " " << validError.at(1) << " \n ";
318 
319  outputFile << "\n ";
320  outputFile << "Image to Probe transform matrix = \n ";
321  for (int m = 0; m < 4; m++)
322  {
323  for (int n = 0; n < 4; n++)
324  {
325  outputFile << imageToProbeTransformMatrix(m, n) << " ";
326  }
327  outputFile << "\n ";
328  }
329 
330  }
331 
332  calibError.clear();
333  validError.clear();
334 
335  }
336 
337  sequenceTrackedFrameList->Clear();
338  }
339 
340  outputFile.close();
341  std::cout << "Exit success!!!" << std::endl;
342  return EXIT_SUCCESS;
343 }
344 
345 //-------------------------------------------------------------------------------------------------
346 
347 PlusStatus SubSequenceMetafile(vtkIGSIOTrackedFrameList* aTrackedFrameList, std::vector<unsigned int> selectedFrames)
348 {
349  LOG_INFO("Create a sub sequence using" << selectedFrames.size() << " frames");
350  std::sort(selectedFrames.begin(), selectedFrames.end());
351  unsigned int FirstFrameIndex = selectedFrames.at(0);
352  unsigned int LastFrameIndex = selectedFrames.at(selectedFrames.size() - 1);
353  if (LastFrameIndex >= aTrackedFrameList->GetNumberOfTrackedFrames() || FirstFrameIndex > LastFrameIndex)
354  {
355  LOG_ERROR("Invalid input range: (" << FirstFrameIndex << ", " << LastFrameIndex << ")" << " Permitted range within (0, " << aTrackedFrameList->GetNumberOfTrackedFrames() - 1 << ")");
356  return PLUS_FAIL;
357  }
358 
359  if (LastFrameIndex != aTrackedFrameList->GetNumberOfTrackedFrames() - 1)
360  {
361  aTrackedFrameList->RemoveTrackedFrameRange(LastFrameIndex + 1, aTrackedFrameList->GetNumberOfTrackedFrames() - 1);
362  }
363 
364  for (int i = selectedFrames.size() - 2; i > 0; i--)
365  {
366  aTrackedFrameList->RemoveTrackedFrameRange(selectedFrames.at(i) + 1, selectedFrames.at(i + 1) - 1);
367  }
368 
369  if (FirstFrameIndex != 0)
370  {
371  aTrackedFrameList->RemoveTrackedFrameRange(0, FirstFrameIndex - 1);
372  }
373  return PLUS_SUCCESS;
374 }
375 
376 //-------------------------------------------------------------------------------------------------
377 PlusStatus SetOptimizationMethod(vtkPlusProbeCalibrationAlgo* freehandCalibration, std::string method)
378 {
379  vtkPlusProbeCalibrationOptimizerAlgo* optimizer = freehandCalibration->GetOptimizer();
380 
381  if (STRCASECMP(method.c_str(), "NO_OPTIMIZATION") == 0)
382  {
384  }
385  else if (STRCASECMP(method.c_str(), "7param2D") == 0)
386  {
388  optimizer->SetIsotropicPixelSpacing(true);
389  }
390  else if (STRCASECMP(method.c_str(), "7param3D") == 0)
391  {
393  optimizer->SetIsotropicPixelSpacing(true);
394  }
395  else if (STRCASECMP(method.c_str(), "8param2D") == 0)
396  {
398  optimizer->SetIsotropicPixelSpacing(false);
399  }
400  else if (STRCASECMP(method.c_str(), "8param3D") == 0)
401  {
403  optimizer->SetIsotropicPixelSpacing(false);
404  }
405  return PLUS_SUCCESS;
406 }
Refines the image to probe transform using non-linear optimization.
PlusFidSegmentation * GetFidSegmentation()
PlusStatus SetOptimizationMethod(vtkPlusProbeCalibrationAlgo *freehandCalibration, std::string method)
void SetIsotropicPixelSpacing(bool isotropicPixelSpacing)
igsioStatus PlusStatus
Definition: PlusCommon.h:40
for i
#define PLUS_FAIL
Definition: PlusCommon.h:43
static vtkPlusConfig * GetInstance()
std::vector< PlusNWire > GetNWires()
#define PLUS_SUCCESS
Definition: PlusCommon.h:44
vtkPlusProbeCalibrationOptimizerAlgo * GetOptimizer()
int main(int argc, char *argv[])
PhidgetLCD_Font int * width
Definition: phidget22.h:4275
void SetDebugOutput(bool value)
static vtkIGSIOLogger * Instance()
PlusFidLineFinder * GetFidLineFinder()
PlusStatus SubSequenceMetafile(vtkIGSIOTrackedFrameList *aTrackedFrameList, std::vector< unsigned int > selectedFrames)
void SetOptimizationMethod(OptimizationMethodType optimizationMethod)
PlusStatus RecognizePattern(vtkIGSIOTrackedFrameList *trackedFrameList, PatternRecognitionError &patternRecognitionError, int *numberOfSuccessfullySegmentedImages=NULL, std::vector< unsigned int > *segmentedFramesIndices=NULL)
void SetDeviceSetConfigurationData(vtkXMLDataElement *deviceSetConfigurationData)
PlusStatus ReadConfiguration(vtkXMLDataElement *rootConfigElement)
static PlusStatus ReadDeviceSetConfigurationFromFile(vtkXMLDataElement *config, const char *filename)
Definition: PlusXmlUtils.h:23
Probe calibration algorithm class.