Utilizing just volume reconstruction algorithm in a program

[SerenaIntusurg]

Hi,

I have an internally developed tracking system (provides 6 DOF frames) and also ultrasound images (2D fan) that are already constructed into an existing program that I have running. I wanted to reconstruct these into a 3D volume. Is there a way I can use the Volume Reconstruction algorithm in a C++ program I already have written? Are there examples of how to pack the localization frames, images, and calibration matrix together in a form that the volume reconstruction can work off of?

Thanks so much!
Serena

[Sunderlandkyl]

The volume reconstruction is implemented in IGSIO (vtkIGSIOVolumeReconstructor). Plus uses a subclass called vtkPlusVolumeReconstructor with only minor modifications.

Does your program already use VTK? If so it would be easy to add the algorithm.

The best example of using the class is probably in SlicerIGT: vtkSlicerVolumeReconstructionLogic.

This function in SlicerIGT is responsible for adding a single slice to the volume reconstructor: https://github.com/SlicerIGT/SlicerIGT/blob/d3fd2b26a4c10f6319aa8ae4173abc793c4113ce/VolumeReconstruction/Logic/vtkSlicerVolumeReconstructionLogic.cxx#L321-L407.
The position of the Slice is defined in the above function by specifying the "ImageToROI" transform for each tracked frame.

Let me know if you have any questions about the code.

[SerenaIntusurg]

Ok - I installed VTK 8.2
Then I cloned IGSIO and tried to configure and build (reference the vtk build), but I'm getting some errors when running cmake that I wasn't sure what to do as it seems not to find vtkAddon (see below)? Is there I'm missing in terms of building in the main Vtk build and how do I build it? I also saw IGSIO/IGSIO#14 this post, but I wasn't sure how to get the IGSIO super build to build?

Thanks so much!
Serena

cmake -DVTK_DIR:PATH=/home/swong/dev/VTK/build ../CMakeLists.txt IGSIO_BUILD_VOLUMERECONSTRUCTION=YES
CMake Warning:
Ignoring extra path from command line:

"IGSIO_BUILD_VOLUMERECONSTRUCTION=YES"

CMake Error at Codecs/CMakeLists.txt:6 (FIND_PACKAGE):
By not providing "FindvtkAddon.cmake" in CMAKE_MODULE_PATH this project has
asked CMake to find a package configuration file provided by "vtkAddon",
but CMake did not find one.

Could not find a package configuration file provided by "vtkAddon" with any
of the following names:

vtkAddonConfig.cmake
vtkaddon-config.cmake

Add the installation prefix of "vtkAddon" to CMAKE_PREFIX_PATH or set
"vtkAddon_DIR" to a directory containing one of the above files. If
"vtkAddon" provides a separate development package or SDK, be sure it has
been installed.

-- Configuring incomplete, errors occurred!
See also "/home/swong/dev/IGSIO/build/CMakeFiles/CMakeOutput.log".

[Sunderlandkyl]

You'll need to build vtkAddon as well. You can find it here: https://github.com/Slicer/vtkAddon.

[SerenaIntusurg]

Ok - I got everything to compile and added it to my project.

I was looking at the setup of the reconstructor in vtkSlicerVolumeRecostructionLogic lines 202-283. Sorry if I missed this, but is there a description of the settings - for example the compoundingmode, optimization, interpolation, and fill hole types?

Thanks so much!

[Sunderlandkyl]

Plus still has some documentation on what the various algorithm parameters do:
http://perk-software.cs.queensu.ca/plus/doc/nightly/user/AlgorithmVolumeReconstruction.html

[SerenaIntusurg]

OK - I finally got everything compiling and code not crashing on my machine in the application I was intending to integrate into, but I'm a little confused I've set up everything correctly. So I was trying to do a small test where I have a volume that is 40 X 40 X 40 mm in size with 0.1 mm spacing set up in the reconstructor. I'm trying to add an image that is all a single value (bright) with transform that is the identity matrix, so I think I would see a single bright frame in the xy plane at the origin, but I don't really see anything in tis volume. Do you have any pointers or suggestions of something I'm doing wrong in the setup? And I'm using vtkImageViewer2 to view and scroll through the slices like in one of the stock examples.

Here is my test function setting everything up:

void VolumeWrapper::addFrameTest ()
{
  
  //initialize the vtk structures
  m_list = vtkIGSIOTrackedFrameList::New();
  m_reconstructor = vtkSmartPointer<vtkIGSIOVolumeReconstructor>::New();
  
  m_outputSpacing[0] = 0.1;//0.5;
  m_outputSpacing[1] = 0.1; //0.5;
  m_outputSpacing[2] = 0.1; //0.5;
  
  m_size[0] = 40.0;
  m_size[1] = 30.0;
  m_size[2] = 20.0;
  
  m_bounds[0] = -m_size[0] / 2.0;
  m_bounds[1] = m_size[0] / 2.0;
  m_bounds[2] = -m_size[1] / 2.0;
  m_bounds[3] = m_size[1] / 2.0;
  m_bounds[4] = -m_size[2] / 2.0;
  m_bounds[5] = m_size[2] / 2.0;
  
  m_outputExtent[0] = 0;
  m_outputExtent[1] = static_cast<int>(std::ceil((m_bounds[1] - m_bounds[0]) / m_outputSpacing[0]));
  m_outputExtent[2] = 0;
  m_outputExtent[3] = static_cast<int>(std::ceil((m_bounds[3] - m_bounds[2]) / m_outputSpacing[1]));
  m_outputExtent[4] = 0;
  m_outputExtent[5] = static_cast<int>(std::ceil((m_bounds[5] - m_bounds[4]) / m_outputSpacing[2]));
  
  m_outputOrigin[0] = m_bounds[0];
  m_outputOrigin[1] = m_bounds[2];
  m_outputOrigin[2] = m_bounds[4];
  
  m_reconstructor->SetOutputExtent(m_outputExtent);
  m_reconstructor->SetOutputOrigin(m_outputOrigin);
  m_reconstructor->SetOutputSpacing(m_outputSpacing);
  
  //http://perk-software.cs.queensu.ca/plus/doc/nightly/user/AlgorithmVolumeReconstruction.html
  //https://github.com/IGSIO/IGSIO/blob/master/VolumeReconstruction/vtkIGSIOVolumeReconstructor.cxx
  m_reconstructor->SetCompoundingMode(vtkIGSIOPasteSliceIntoVolume::MEAN_COMPOUNDING_MODE);
  m_reconstructor->SetOptimization(vtkIGSIOPasteSliceIntoVolume::NO_OPTIMIZATION);
  m_reconstructor->SetInterpolation(vtkIGSIOPasteSliceIntoVolume::LINEAR_INTERPOLATION);
  m_reconstructor->SetNumberOfThreads(0);
  
  
  m_reconstructor->SetFillHoles(m_isFillHoles);
  if (m_isFillHoles)
  {
    vtkIGSIOFillHolesInVolume* holeFiller = m_reconstructor->GetHoleFiller();
    holeFiller->SetNumHFElements(1);
    holeFiller->AllocateHFElements();
    FillHolesInVolumeElement hfElement;
    hfElement.setupAsStick(9, 1);
    holeFiller->SetHFElement(0, hfElement);
  }
  m_reconstructor->SetImageCoordinateFrame("Image");
  m_reconstructor->SetReferenceCoordinateFrame("ROI");
  
  
  m_clipRectOrigin [0] = 0;
  m_clipRectOrigin [1] = 0;
  
  m_clipRectSize[0] = 400;
  m_clipRectSize[1] = 400;
  
  m_fanAnglesDeg [0] = -35;
  m_fanAnglesDeg [1] = 35;
  m_fanOriginPixel [0] = 200;
  m_fanOriginPixel [1] = 10;
  
  m_fanRadiusStartPixel = 0.0;
  
  m_reconstructor->SetClipRectangleOrigin(m_clipRectOrigin);
  m_reconstructor->SetClipRectangleSize(m_clipRectSize);
  
  
  m_reconstructor->SetFanAnglesDeg(m_fanAnglesDeg);
  m_reconstructor->SetFanOriginPixel(m_fanOriginPixel);
  m_reconstructor->SetFanRadiusStartPixel(m_fanRadiusStartPixel);

  m_reconstructor->Reset();
  
  //set transform to identity - so the image should be at the origin of the volume in the xy plane?
  vtkTransform* imageToROITransform = vtkTransform::New();  //this is the transform that is associated with an image  
  imageToROITransform->Identity();
  
  vtkMatrix4x4* matrix;
  matrix = imageToROITransform->GetMatrix();  
  LOG_MSG_DBG ("MATRIX of XFORM: %f %f %f %f, %f %f %f %f, %f %f %f %f, %f %f %f %f\n",
    matrix->GetElement (0,0), matrix->GetElement (0,1), matrix->GetElement (0,2), matrix->GetElement (0,3),
    matrix->GetElement (1,0), matrix->GetElement (1,1), matrix->GetElement (1,2), matrix->GetElement (1,3),
    matrix->GetElement (2,0), matrix->GetElement (2,1), matrix->GetElement (2,2), matrix->GetElement (2,3),
    matrix->GetElement (3,0), matrix->GetElement (3,1), matrix->GetElement (3,2), matrix->GetElement (3,3));
  
  m_transformRepository = vtkSmartPointer<vtkIGSIOTransformRepository>::New();  //add to image repository
  m_transformRepository->SetTransform(igsioTransformName("ImageToROI"), imageToROITransform->GetMatrix());
  
  //just create an image of all bright data
  vtkImageData* inputImageData; 
  inputImageData = vtkImageData::New();
  inputImageData->SetDimensions(400,400,1);
  inputImageData->SetSpacing(0.1,0.1,0.1);
  
  vtkInformation* information = vtkInformation::New();
  inputImageData->SetNumberOfScalarComponents(1,information);
  
  int* sizes = inputImageData->GetDimensions();
  int components = inputImageData->GetNumberOfScalarComponents();
  int scalarsize = inputImageData->GetScalarSize();
  
  for (int i = 0; i < sizes[0]; i++)
  {
    for (int j = 0; j < sizes[1]; j++)
    {
      for (int k = 0; k < sizes[2]; k++)
      {
        inputImageData->SetScalarComponentFromDouble(i,j,k,components-1, 800.0);
      }
    }
  }
  
  if (inputImageData)
  {
    m_reconstructor->SetOutputScalarType(inputImageData->GetScalarType());
    LOG_MSG_DBG ("scalar type is %s", inputImageData->GetScalarTypeAsString());
    //https://vtk.org/doc/nightly/html/vtkType_8h_source.html
  }
  
  std::string errorDetail;
  
  m_trackedFrame = new igsioTrackedFrame();
  //igsioVideoFrame* videoframe = new igsioVideoFrame();
  
  FrameSizeType frameSize = { static_cast<unsigned int>(sizes[0]), static_cast<unsigned int>(sizes[1]), static_cast<unsigned int>(sizes[2]) };
  igsioTrackedFrame frame;
  igsioVideoFrame image;
  image.AllocateFrame(frameSize, VTK_DOUBLE, 1);
  image.GetImage()->DeepCopy(inputImageData);
  image.SetImageOrientation(US_IMG_ORIENT_MF);
  image.SetImageType(US_IMG_BRIGHTNESS);
  frame.SetImageData(image);
  m_list->AddTrackedFrame(&frame);


  //insert int volume
  bool isLast = true;
  bool insertedIntoVolume;
  
//  igsioStatus status = m_reconstructor->AddTrackedFrame(m_trackedFrame, m_transformRepository, m_isFirstFrame, isLast, &insertedIntoVolume); // != IGSIO_SUCCESS)
  igsioStatus status = m_reconstructor->AddTrackedFrame(&frame, m_transformRepository, m_isFirstFrame, isLast, &insertedIntoVolume); // != IGSIO_SUCCESS)

  
  if (insertedIntoVolume)
    LOG_MSG_DBG ("SUCCESSFULLY INSERTED INTO VOLUME\n");
  else 
    LOG_MSG_DBG ("COULDN'T INSERT INTO VOLUME\n");
  
  
  if (m_isFirstFrame)
    m_isFirstFrame = false;
    
  //copy the image out and then somewhere else use the vtkImageViewer2 to look through the SetClipRectangleSize
  outputData = vtkImageData::New();
  m_reconstructor->GetReconstructedVolume(outputData);

}

Then I'm setting the input to the viewer and the level and window in another function:

           USVolumeMapper_.addFrameTest ();
            imageViewer_->SetInputData(ionUSVolumeMapper_.changeImageData ()); //gets pointer to the outputData
            imageViewer_->SetColorWindow(900.0);
            imageViewer_->SetColorLevel(450.0);

I basically get dark images like this one through all slices and I've looked in all three cardinal directions too (XY, XZ, and YZ planes)

Thanks so much in advance!