PlusLib  2.9.0
Software library for tracked ultrasound image acquisition, calibration, and processing.
vtkPlusCenterOfRotationCalibAlgo.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 
9 #include "PlusMath.h"
10 #ifdef PLUS_RENDERING_ENABLED
11 #include "PlusPlotter.h"
12 #endif
14 #include "vtkObjectFactory.h"
15 #include "vtkIGSIOTrackedFrameList.h"
16 #include "igsioTrackedFrame.h"
17 #include "vtkPoints.h"
18 #include "vtksys/SystemTools.hxx"
19 #include "vtkPlusHTMLGenerator.h"
20 #include "vtkDoubleArray.h"
21 #include "vtkVariantArray.h"
22 
23 //----------------------------------------------------------------------------
25 
26 //----------------------------------------------------------------------------
28 {
29  this->TrackedFrameList = NULL;
30  this->SetCenterOfRotationPx(0, 0);
31  this->SetSpacing(0, 0);
32  this->ReportTable = NULL;
33  this->ErrorMean = 0.0;
34  this->ErrorStdev = 0.0;
35 }
36 
37 //----------------------------------------------------------------------------
39 {
40  this->SetTrackedFrameList(NULL);
41  this->SetReportTable(NULL);
42 }
43 
44 //----------------------------------------------------------------------------
45 void vtkPlusCenterOfRotationCalibAlgo::PrintSelf(ostream& os, vtkIndent indent)
46 {
47  os << std::endl;
48  this->Superclass::PrintSelf(os, indent);
49  os << indent << "Update time: " << UpdateTime.GetMTime() << std::endl;
50  os << indent << "Spacing: " << this->Spacing[0] << " " << this->Spacing[1] << std::endl;
51  os << indent << "Center of rotation (px): " << this->CenterOfRotationPx[0] << " " << this->CenterOfRotationPx[1] << std::endl;
52  os << indent << "Calibration error: mean=" << this->ErrorMean << " stdev=" << this->ErrorStdev << std::endl;
53 
54  if (this->TrackedFrameList != NULL)
55  {
56  os << indent << "TrackedFrameList:" << std::endl;
57  this->TrackedFrameList->PrintSelf(os, indent);
58  }
59 
60  if (this->ReportTable != NULL)
61  {
62  os << indent << "ReportTable:" << std::endl;
63  this->ReportTable->PrintSelf(os, indent);
64  }
65 }
66 
67 
68 //----------------------------------------------------------------------------
69 void vtkPlusCenterOfRotationCalibAlgo::SetInputs(vtkIGSIOTrackedFrameList* trackedFrameList, std::vector<int>& indices, double spacing[2])
70 {
71  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::SetInputs");
72  this->SetTrackedFrameList(trackedFrameList);
73  this->SetSpacing(spacing);
74  this->SetTrackedFrameListIndices(indices);
75 }
76 
77 //----------------------------------------------------------------------------
79 {
80  this->TrackedFrameListIndices = indices;
81  this->Modified();
82 }
83 
84 
85 //----------------------------------------------------------------------------
87 {
88  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::GetCenterOfRotationPx");
89  // Update calibration result
90  PlusStatus status = this->Update();
91 
92  centerOfRotationPx[0] = this->CenterOfRotationPx[0];
93  centerOfRotationPx[1] = this->CenterOfRotationPx[1];
94 
95  return status;
96 }
97 
98 //----------------------------------------------------------------------------
100 {
101  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::GetError");
102  // Update calibration result
103  PlusStatus status = this->Update();
104 
105  mean = this->ErrorMean;
106  stdev = this->ErrorStdev;
107 
108  return status;
109 }
110 
111 //----------------------------------------------------------------------------
113 {
114  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::Update");
115 
116  if (this->GetMTime() < this->UpdateTime.GetMTime())
117  {
118  LOG_DEBUG("Center of rotation calibration result is up-to-date!");
119  return PLUS_SUCCESS;
120  }
121 
122  // Check if TrackedFrameList is MF oriented BRIGHTNESS image
123  if (vtkIGSIOTrackedFrameList::VerifyProperties(this->TrackedFrameList, US_IMG_ORIENT_MF, US_IMG_BRIGHTNESS) != PLUS_SUCCESS)
124  {
125  LOG_ERROR("Failed to perform calibration - tracked frame list is invalid");
126  return PLUS_FAIL;
127  }
128 
129  // Data containers
130  std::vector<vnl_vector<double> > aMatrix;
131  std::vector<double> bVector;
132 
133  if (this->ConstructLinearEquationForCalibration(aMatrix, bVector) != PLUS_SUCCESS)
134  {
135  LOG_ERROR("Unable to construct linear equation for center of rotation calibration algorithm!");
136  return PLUS_FAIL;
137  }
138 
139  if (aMatrix.size() == 0 || bVector.size() == 0)
140  {
141  LOG_WARNING("Center of rotation calculation failed, no data found!");
142  return PLUS_FAIL;
143  }
144 
145  // The rotation center in original image frame in px
146  vnl_vector<double> centerOfRotationInPx(2, 0);
147  if (PlusMath::LSQRMinimize(aMatrix, bVector, centerOfRotationInPx, &this->ErrorMean, &this->ErrorStdev) != PLUS_SUCCESS)
148  {
149  LOG_WARNING("Failed to run LSQRMinimize!");
150  return PLUS_FAIL;
151  }
152 
153  if (centerOfRotationInPx.empty() || centerOfRotationInPx.size() < 2)
154  {
155  LOG_ERROR("Unable to calibrate center of rotation! Minimizer returned empty result.");
156  return PLUS_FAIL;
157  }
158 
159  // Set center of rotation - don't use set macro, it changes the modification time of the algorithm
160  this->CenterOfRotationPx[0] = centerOfRotationInPx[0] / this->Spacing[0];
161  this->CenterOfRotationPx[1] = centerOfRotationInPx[1] / this->Spacing[1];
162 
163  this->UpdateReportTable();
164 
165  this->UpdateTime.Modified();
166 
167  return PLUS_SUCCESS;
168 
169 }
170 
171 //----------------------------------------------------------------------------
173 {
174  int numberOfNWirePatterns = 0;
175 
176  if (this->TrackedFrameList == NULL)
177  {
178  LOG_WARNING("Unable to compute number of N wire patterns - tracked frame list is NULL!");
179  return numberOfNWirePatterns;
180  }
181 
182  for (unsigned int frame = 0; frame < this->TrackedFrameList->GetNumberOfTrackedFrames(); ++frame)
183  {
184  igsioTrackedFrame* trackedFrame = this->TrackedFrameList->GetTrackedFrame(frame);
185  if (trackedFrame == NULL)
186  {
187  LOG_ERROR("Unable to get tracked frame from the list - tracked frame is NULL (position in the list: " << frame << ")!");
188  continue;
189  }
190 
191  vtkPoints* fiduacialPointsCoordinatePx = trackedFrame->GetFiducialPointsCoordinatePx();
192  if (fiduacialPointsCoordinatePx == NULL)
193  {
194  LOG_ERROR("Unable to get segmented fiducial points from tracked frame - FiducialPointsCoordinatePx is NULL, frame is not yet segmented (position in the list: " << frame << ")!");
195  continue;
196  }
197 
198  if (fiduacialPointsCoordinatePx->GetNumberOfPoints() == 0)
199  {
200  //LOG_DEBUG("Unable to get segmented fiducial points from tracked frame - couldn't segment image (position in the list: " << frame << ")!" );
201  continue;
202  }
203 
204  // Make sure we have 3 points for each N wire pattern
205  if (fiduacialPointsCoordinatePx->GetNumberOfPoints() % 3 != 0)
206  {
207  LOG_WARNING("Unrecognized wire patterns for frame #" << frame);
208  continue;
209  }
210 
211  // Every N fiducial has 3 points on the image: numberOfNFiducials = NumberOfPoints / 3
212  int patterns = fiduacialPointsCoordinatePx->GetNumberOfPoints() / 3;
213 
214  // Make sure all the frames has the same number of patterns
215  if (numberOfNWirePatterns > 0 && patterns != numberOfNWirePatterns)
216  {
217  LOG_ERROR("Number of N-wire patterns different between frames!"); // this shouldn't happen
218  return 0;
219  }
220 
221  numberOfNWirePatterns = patterns;
222  }
223 
224  return numberOfNWirePatterns;
225 }
226 
227 //----------------------------------------------------------------------------
228 PlusStatus vtkPlusCenterOfRotationCalibAlgo::ConstructLinearEquationForCalibration(std::vector<vnl_vector<double> >& aMatrix, std::vector<double>& bVector)
229 {
230  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::ConstructLinearEquationForCalibration");
231  aMatrix.clear();
232  bVector.clear();
233 
234  if (this->TrackedFrameList == NULL)
235  {
236  LOG_ERROR("Failed to construct linear equation for center of rotation calibration - tracked frame list is NULL!");
237  return PLUS_FAIL;
238  }
239 
240  const int numberOfNFiduacials = this->GetNumberOfNWirePatterns();
241 
242  //******************* Count the number of frames used for calibration *****************
243  int numberOfFrames = 0;
244  for (unsigned int index = 0; index < this->TrackedFrameListIndices.size(); ++index)
245  {
246  // Get the next frame index used for calibration
247  int frameNumber = this->TrackedFrameListIndices[index];
248 
249  // Get tracked frame from list
250  igsioTrackedFrame* trackedFrame = this->TrackedFrameList->GetTrackedFrame(frameNumber);
251 
252  // Skip frame if it was not segmented
253  if (trackedFrame->GetFiducialPointsCoordinatePx() == NULL)
254  {
255  continue;
256  }
257 
258  // Skip frame if the segmentation was not successful
259  if (trackedFrame->GetFiducialPointsCoordinatePx()->GetNumberOfPoints() == 0)
260  {
261  continue;
262  }
263 
264  numberOfFrames++;
265  }
266 
267  if (numberOfFrames < 30)
268  {
269  LOG_WARNING("Center of rotation calculation failed - there is not enough data (" << numberOfFrames << " out of at least 30)!");
270  return PLUS_FAIL;
271  }
272 
273  // Reserve memory for wire points
274  std::vector<vtkSmartPointer<vtkPoints> > vectorOfWirePoints;
275  vectorOfWirePoints.reserve(numberOfFrames);
276 
277  for (unsigned int i = 0; i < this->TrackedFrameListIndices.size(); ++i)
278  {
279  // Get the next frame index used for calibration
280  int frameNumber = this->TrackedFrameListIndices[i];
281 
282  // Get tracked frame from list
283  igsioTrackedFrame* trackedFrame = this->TrackedFrameList->GetTrackedFrame(frameNumber);
284 
285  if (trackedFrame->GetFiducialPointsCoordinatePx() == NULL)
286  {
287  LOG_ERROR("Unable to get segmented fiducial points from tracked frame - FiducialPointsCoordinatePx is NULL, frame is not yet segmented (position in the list: " << frameNumber << ")!");
288  continue;
289  }
290 
291  if (trackedFrame->GetFiducialPointsCoordinatePx()->GetNumberOfPoints() == 0)
292  {
293  continue;
294  }
295 
296  vtkSmartPointer<vtkPoints> point = vtkSmartPointer<vtkPoints>::New();
297  point->SetDataTypeToDouble();
298  point->SetNumberOfPoints(numberOfNFiduacials * 2); // Use only the two non-moving points of the N fiducial
299 
300  int vectorID(0); // ID used for point position in vtkPoints for faster execution
301  for (int p = 0; p < trackedFrame->GetFiducialPointsCoordinatePx()->GetNumberOfPoints(); ++p)
302  {
303  if (((p + 1) % 3) != 2) // wire #1,#3,#4,#6... => use only non moving points of the N-wire
304  {
305  double* wireCoordinatePx = trackedFrame->GetFiducialPointsCoordinatePx()->GetPoint(p);
306  point->SetPoint(vectorID++, wireCoordinatePx[0], wireCoordinatePx[1], wireCoordinatePx[2]);
307  }
308  }
309 
310  vectorOfWirePoints.push_back(point);
311  }
312 
313  for (unsigned int i = 0; i <= vectorOfWirePoints.size() - 2; i++)
314  {
315  for (unsigned int j = i + 1; j <= vectorOfWirePoints.size() - 1; j = j + 2)
316  {
317  if (vectorOfWirePoints[i]->GetNumberOfPoints() != vectorOfWirePoints[j]->GetNumberOfPoints())
318  {
319  continue;
320  }
321 
322  for (int point = 0; point < vectorOfWirePoints[i]->GetNumberOfPoints(); point++)
323  {
324  // coordiates of the i-th element
325  double Xi = vectorOfWirePoints[i]->GetPoint(point)[0] * this->Spacing[0];
326  double Yi = vectorOfWirePoints[i]->GetPoint(point)[1] * this->Spacing[1];
327 
328  // coordiates of the j-th element
329  double Xj = vectorOfWirePoints[j]->GetPoint(point)[0] * this->Spacing[0];
330  double Yj = vectorOfWirePoints[j]->GetPoint(point)[1] * this->Spacing[1];
331 
332  // Populate the list of distance
333  vnl_vector<double> rowOfDistance(2, 0);
334  rowOfDistance.put(0, Xi - Xj);
335  rowOfDistance.put(1, Yi - Yj);
336  aMatrix.push_back(rowOfDistance);
337 
338  // Populate the squared distance vector
339  bVector.push_back(0.5 * (Xi * Xi + Yi * Yi - Xj * Xj - Yj * Yj));
340  }
341  }
342  }
343 
344  return PLUS_SUCCESS;
345 }
346 
347 //----------------------------------------------------------------------------
349 {
350  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::UpdateReportTable");
351 
352  // Clear table before update
353  this->SetReportTable(NULL);
354 
355  if (this->ReportTable == NULL)
356  {
357  this->AddNewColumnToReportTable("ProbePosition");
358  this->AddNewColumnToReportTable("ProbeRotation");
359  this->AddNewColumnToReportTable("TemplatePosition");
360 
361  for (int i = 0; i < this->GetNumberOfNWirePatterns(); ++i)
362  {
363  std::ostringstream columnNameRightRadius;
364  columnNameRightRadius << "Wire#" << (i * 3 + 1) << "Radius";
365  this->AddNewColumnToReportTable(columnNameRightRadius.str().c_str());
366 
367  std::ostringstream columnNameRightRadiusDistanceFromMean;
368  columnNameRightRadiusDistanceFromMean << "Wire#" << (i * 3 + 1) << "RadiusDistanceFromMean";
369  this->AddNewColumnToReportTable(columnNameRightRadiusDistanceFromMean.str().c_str());
370 
371  std::ostringstream columnNameRightWireX;
372  columnNameRightWireX << "w" << (i * 3 + 1) << "xPx";
373  this->AddNewColumnToReportTable(columnNameRightWireX.str().c_str());
374 
375  std::ostringstream columnNameRightWireY;
376  columnNameRightWireY << "w" << (i * 3 + 1) << "yPx";
377  this->AddNewColumnToReportTable(columnNameRightWireY.str().c_str());
378 
379  std::ostringstream columnNameLeftRadius;
380  columnNameLeftRadius << "Wire#" << (i * 3 + 3) << "Radius";
381  this->AddNewColumnToReportTable(columnNameLeftRadius.str().c_str());
382 
383  std::ostringstream columnNameLeftRadiusDistanceFromMean;
384  columnNameLeftRadiusDistanceFromMean << "Wire#" << (i * 3 + 3) << "RadiusDistanceFromMean";
385  this->AddNewColumnToReportTable(columnNameLeftRadiusDistanceFromMean.str().c_str());
386 
387  std::ostringstream columnNameLeftWireX;
388  columnNameLeftWireX << "w" << (i * 3 + 3) << "xPx";
389  this->AddNewColumnToReportTable(columnNameLeftWireX.str().c_str());
390 
391  std::ostringstream columnNameLeftWireY;
392  columnNameLeftWireY << "w" << (i * 3 + 3) << "yPx";
393  this->AddNewColumnToReportTable(columnNameLeftWireY.str().c_str());
394  }
395  }
396 
397  const double sX = this->Spacing[0];
398  const double sY = this->Spacing[1];
399  std::vector<std::vector<double> > wireRadiusVector(this->GetNumberOfNWirePatterns() * 2); // each wire has two points
400  std::vector<std::vector<double> > wirePositions(this->GetNumberOfNWirePatterns() * 4); // each wire has 2 point and each point has 2 coordinates
401 
402  std::vector<double> probePosVector;
403  std::vector<double> probeRotVector;
404  std::vector<double> templatePosVector;
405 
406  for (unsigned int i = 0; i < this->TrackedFrameListIndices.size(); i++)
407  {
408  const int frameNumber = this->TrackedFrameListIndices[i];
409  igsioTrackedFrame* frame = this->TrackedFrameList->GetTrackedFrame(frameNumber);
410 
411  if (frame->GetFiducialPointsCoordinatePx() == NULL
412  || frame->GetFiducialPointsCoordinatePx()->GetNumberOfPoints() == 0)
413  {
414  // This frame was not segmented
415  continue;
416  }
417 
418  double probePos(0), probeRot(0), templatePos(0);
419  if (!igsioTrackedFrameEncoderPositionFinder::GetStepperEncoderValues(frame, probePos, probeRot, templatePos))
420  {
421  LOG_WARNING("Unable to get probe position from tracked frame info for frame #" << frameNumber);
422  continue;
423  }
424  probePosVector.push_back(probePos);
425  probeRotVector.push_back(probeRot);
426  templatePosVector.push_back(templatePos);
427 
428  // Compute radius
429  for (int w = 0; w < this->GetNumberOfNWirePatterns(); ++w)
430  {
431  double w1x = frame->GetFiducialPointsCoordinatePx()->GetPoint(w * 3)[0];
432  double w1y = frame->GetFiducialPointsCoordinatePx()->GetPoint(w * 3)[1];
433  double w1Radius = sqrt(pow((w1x - this->CenterOfRotationPx[0]) * sX, 2) + pow((w1y - this->CenterOfRotationPx[1]) * sY, 2));
434 
435  wirePositions[4 * w].push_back(w1x);
436  wirePositions[4 * w + 1].push_back(w1y);
437  wireRadiusVector[2 * w].push_back(w1Radius);
438 
439 
440  double w3x = frame->GetFiducialPointsCoordinatePx()->GetPoint(w * 3 + 2)[0];
441  double w3y = frame->GetFiducialPointsCoordinatePx()->GetPoint(w * 3 + 2)[1];
442  double w3Radius = sqrt(pow((w3x - this->CenterOfRotationPx[0]) * sX, 2) + pow((w3y - this->CenterOfRotationPx[1]) * sY, 2));
443 
444  wirePositions[4 * w + 2].push_back(w3x);
445  wirePositions[4 * w + 3].push_back(w3y);
446  wireRadiusVector[2 * w + 1].push_back(w3Radius);
447  }
448  }
449 
450  const int numberOfElements = wireRadiusVector[0].size();
451  std::vector<double> wireRadiusMean(this->GetNumberOfNWirePatterns() * 2, 0); // each wire has two points
452 
453  for (int i = 0; i < numberOfElements; ++i)
454  {
455  for (int w = 0; w < this->GetNumberOfNWirePatterns() * 2; ++w)
456  {
457  wireRadiusMean[w] += wireRadiusVector[w][i] / (1.0 * numberOfElements);
458  }
459  }
460 
461  std::vector<std::vector<double> > wireDistancesFromMeanRadius(this->GetNumberOfNWirePatterns() * 2); // each wire has two points
462  for (int i = 0; i < numberOfElements; ++i)
463  {
464  for (int w = 0; w < this->GetNumberOfNWirePatterns() * 2; ++w)
465  {
466  wireDistancesFromMeanRadius[w].push_back(wireRadiusMean[w] - wireRadiusVector[w][i]);
467  }
468  }
469 
470  for (int row = 0; row < numberOfElements; ++row)
471  {
472  vtkSmartPointer<vtkVariantArray> tableRow = vtkSmartPointer<vtkVariantArray>::New();
473 
474  tableRow->InsertNextValue(probePosVector[row]); //ProbePosition
475  tableRow->InsertNextValue(probeRotVector[row]); //ProbeRotation
476  tableRow->InsertNextValue(templatePosVector[row]); //TemplatePosition
477 
478  for (int w = 0; w < this->GetNumberOfNWirePatterns() * 2; ++w)
479  {
480  tableRow->InsertNextValue(wireRadiusVector[w][row]); //Wire Radius
481  tableRow->InsertNextValue(wireDistancesFromMeanRadius[w][row]); //Wire RadiusDistanceFromMean
482  tableRow->InsertNextValue(wirePositions[2 * w][row]); //wire pixel coordinate x
483  tableRow->InsertNextValue(wirePositions[2 * w + 1][row]); //wire pixel coordinate y
484  }
485 
486  if (tableRow->GetNumberOfTuples() == this->ReportTable->GetNumberOfColumns())
487  {
488  this->ReportTable->InsertNextRow(tableRow);
489  }
490  else
491  {
492  LOG_WARNING("Unable to insert new row to center of rotation error table, number of columns are different ("
493  << tableRow->GetNumberOfTuples() << " vs. " << this->ReportTable->GetNumberOfColumns() << ").");
494  }
495 
496  }
497 
498 #ifdef PLUS_RENDERING_ENABLED
499  PlusPlotter::WriteTableToFile(*this->ReportTable, "RotationAxisCalibrationErrorReport.txt");
500 #endif
501 
502  return PLUS_SUCCESS;
503 }
504 
505 //----------------------------------------------------------------------------
507 {
508  if (this->ReportTable == NULL)
509  {
510  vtkSmartPointer<vtkTable> table = vtkSmartPointer<vtkTable>::New();
511  this->SetReportTable(table);
512  }
513 
514  if (columnName == NULL)
515  {
516  LOG_ERROR("Failed to add new column to table - column name is NULL!");
517  return PLUS_FAIL;
518  }
519 
520  if (this->ReportTable->GetColumnByName(columnName) != NULL)
521  {
522  LOG_WARNING("Column name " << columnName << " already exist!");
523  return PLUS_FAIL;
524  }
525 
526  // Create table column
527  vtkSmartPointer<vtkDoubleArray> col = vtkSmartPointer<vtkDoubleArray>::New();
528  col->SetName(columnName);
529  this->ReportTable->AddColumn(col);
530 
531  return PLUS_SUCCESS;
532 }
533 
534 //----------------------------------------------------------------------------
536 {
537  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::GenerateReport");
538 
539  // Update result before report generation
540  if (this->Update() != PLUS_SUCCESS)
541  {
542  LOG_ERROR("Unable to generate report - center of rotation axis calibration failed!");
543  return PLUS_FAIL;
544  }
545 
547 }
548 
549 //----------------------------------------------------------------------------
550 //static
552  vtkPlusHTMLGenerator* htmlReport, vtkTable* reportTable, double centerOfRotationPx[2])
553 {
554  LOG_TRACE("vtkPlusCenterOfRotationCalibAlgo::GenerateCenterOfRotationReport");
555 
556  if (htmlReport == NULL)
557  {
558  LOG_ERROR("Caller should define HTML report generator before report generation!");
559  return PLUS_FAIL;
560  }
561 
562  if (reportTable == NULL)
563  {
564  LOG_ERROR("Unable to generate report - input report table is NULL!");
565  return PLUS_FAIL;
566  }
567 
568 
569  std::string reportFileName = vtkPlusConfig::GetInstance()->GetApplicationStartTimestamp() + ".CenterOfRotationCalculationError.txt";
570  std::string reportFile = vtkPlusConfig::GetInstance()->GetOutputPath(reportFileName);
571 #ifdef PLUS_RENDERING_ENABLED
572  if (PlusPlotter::WriteTableToFile(*reportTable, reportFile.c_str()) != PLUS_SUCCESS)
573  {
574  LOG_ERROR("Failed to dump translation axis calibration report table to " << reportFile);
575  return PLUS_FAIL;
576  }
577 #endif
578  htmlReport->AddText("Center of Rotation Calculation Analysis", vtkPlusHTMLGenerator::H1);
579 
580  std::ostringstream report;
581  report << "Center of rotation (px): " << centerOfRotationPx[0] << " " << centerOfRotationPx[1] << "</br>" ;
582  htmlReport->AddParagraph(report.str().c_str());
583 
584  // Every N wire has 2 plots, one for each side wire
585  const int numOfPlots = numberOfNWirePatterns * 2;
586 
587  // Create wire labels
588  std::vector<int> wires(numOfPlots); // = {1, 3, 4, 6, ... };
589  for (int i = 0; i < numberOfNWirePatterns; ++i)
590  {
591  wires[2 * i] = i * 3 + 1;
592  wires[2 * i + 1] = i * 3 + 3;
593  }
594 
595  for (int i = 0; i < numOfPlots; i++)
596  {
597  std::ostringstream wireName;
598  wireName << wires[i];
599  std::string wireFullName = std::string("Wire #") + wireName.str();
600  htmlReport->AddText(wireFullName.c_str(), vtkPlusHTMLGenerator::H3);
601 
602  double valueRangeMin = -2.0;
603  double valueRangeMax = 2.0;
604  int numberOfBins = 41;
605  int imageSize[2] = {800, 400};
606 
607  std::string outputImageFilename = std::string("w") + wireName.str() + "-PositionError.png";
608  std::string outputImagePath = htmlReport->AddImageAutoFilename(outputImageFilename.c_str(), "Wire distance from transducer");
609 #ifdef PLUS_RENDERING_ENABLED
610  PlusPlotter::WriteScatterChartToFile("Probe position (mm)", "Wire distance from transducer (mm)", *reportTable, 0 /* ProbePosition */, 3 + i * 4 /* Wire#xRadius */, imageSize, outputImagePath.c_str());
611 
612  outputImageFilename = std::string("w") + wireName.str() + "-AngleError.png";
613  outputImagePath = htmlReport->AddImageAutoFilename(outputImageFilename.c_str(), "Wire distance from transducer");
614  PlusPlotter::WriteScatterChartToFile("Probe angle (deg)", "Wire distance from transducer (mm)", *reportTable, 1 /* ProbeAngle */, 3 + i * 4 /* Wire#xRadius */, imageSize, outputImagePath.c_str());
615 
616  outputImageFilename = std::string("w") + wireName.str() + "-ErrorHistogram.png";
617  outputImagePath = htmlReport->AddImageAutoFilename(outputImageFilename.c_str(), "Wire distance from transducer - error histogram");
618  PlusPlotter::WriteHistogramChartToFile("Wire distance from transducer - error histogram", *reportTable, 4 + i * 4 /* Wire#xRadius */, valueRangeMin, valueRangeMax, numberOfBins, imageSize, outputImagePath.c_str());
619 #endif
620  }
621 
622  htmlReport->AddHorizontalLine();
623 
624  return PLUS_SUCCESS;
625 }
static PlusStatus WriteTableToFile(vtkTable &table, const std::string &filename)
virtual void AddText(const char *text, HEADINGS h, const char *style=NULL)
std::string GetOutputPath(const std::string &subPath)
static PlusStatus LSQRMinimize(const std::vector< std::vector< double > > &aMatrix, const std::vector< double > &bVector, vnl_vector< double > &resultVector, double *mean=NULL, double *stdev=NULL, vnl_vector< unsigned int > *notOutliersIndices=NULL)
Definition: PlusMath.cxx:34
virtual PlusStatus GenerateReport(vtkPlusHTMLGenerator *htmlReport)
igsioStatus PlusStatus
Definition: PlusCommon.h:40
for i
Algorithm for computing the center of in-plane rotation of ultrasound images of a calibration phantom...
virtual void SetTrackedFrameList(vtkIGSIOTrackedFrameList *)
#define PLUS_FAIL
Definition: PlusCommon.h:43
virtual void PrintSelf(ostream &os, vtkIndent indent) VTK_OVERRIDE
static vtkPlusConfig * GetInstance()
virtual void SetCenterOfRotationPx(int, int)
Image slice number p
Definition: algo4.m:14
static PlusStatus GenerateCenterOfRotationReport(int numberOfNWirePatterns, vtkPlusHTMLGenerator *htmlReport, vtkTable *reportTable, double centerOfRotationPx[2])
#define PLUS_SUCCESS
Definition: PlusCommon.h:44
virtual void AddHorizontalLine()
virtual void SetTrackedFrameListIndices(std::vector< int > &indices)
virtual void SetReportTable(vtkTable *)
vtkStandardNewMacro(vtkPlusCenterOfRotationCalibAlgo)
virtual PlusStatus GetCenterOfRotationPx(double centerOfRotationPx[2])
virtual void AddParagraph(const char *paragraph)
static PlusStatus WriteHistogramChartToFile(const std::string &chartTitle, vtkTable &inputTable, int inputColumnIndex, double valueRangeMin, double valueRangeMax, int numberOfBins, int imageSize[2], const std::string &outputImageFilename)
Definition: PlusPlotter.cxx:95
virtual std::string AddImageAutoFilename(const char *filenamePostfix, const char *description, const int widthPx=0, const int heightPx=0)
virtual PlusStatus ConstructLinearEquationForCalibration(std::vector< vnl_vector< double > > &aMatrix, std::vector< double > &bVector)
virtual void SetSpacing(double, double)
class for generating basic html tags
PlusStatus AddNewColumnToReportTable(const char *columnName)
static PlusStatus WriteScatterChartToFile(const std::string &chartTitle, const std::string &yAxisText, vtkTable &inputTable, int xColumnIndex, int yColumnIndex, int imageSize[2], const std::string &outputImageFilename)
Definition: PlusPlotter.cxx:31
virtual void SetInputs(vtkIGSIOTrackedFrameList *trackedFrameList, std::vector< int > &indices, double spacing[2])
virtual PlusStatus GetError(double &mean, double &stdev)