[text] Smoothing operatiobn

Viewer

copydownloadembedprintName: Smoothing operatiobn
  1. MeshBuilder& MeshBuilder::AddGaussianSmoothingPreprocessOperation(
  2.    const SmoothingImageConfiguration& configuration)
  3. {
  4.    mImagesOperations.emplace_back(
  5.       std::make_shared<GaussianSmoothingImageOperation>(configuration));
  6.    return *this;
  7. }
  8.  
  9. vtkSmartPointer<vtkImageData> GaussianSmoothingImageOperation::Exec(
  10.    const vtkSmartPointer<vtkImageData>& imageData)
  11. {
  12.    ORTHO_LOG_INFO("GaussianSmoothingImageOperation - Start");
  13.    // The padding is now adjusted automatically 
  14.    double spacingImage[3];
  15.    imageData->GetSpacing(spacingImage);
  16.  
  17.    // Compute the image padding
  18.    double paddingImage[3] = { 
  19.       mConfiguration.mGaussianSmoothingActualRadiusX / spacingImage[0],
  20.       mConfiguration.mGaussianSmoothingActualRadiusY / spacingImage[1],
  21.       mConfiguration.mGaussianSmoothingActualRadiusZ / spacingImage[2] };
  22.  
  23.    // pad the input volume because if the segmentation touches the edge of the dataset,
  24.    // the generated mesh will contain holes
  25.    vtkSmartPointer<vtkImageConstantPad> padImageFilter = vtkImageConstantPad::New();
  26.    padImageFilter->SetInputData(imageData);
  27.    padImageFilter->SetConstant(0);
  28.    int extent[6];
  29.    imageData->GetExtent(extent);
  30.    double space[3];
  31.    imageData->GetSpacing(space);
  32.    std::cout << "image data" << space[0] << std::endl;
  33.  
  34.    // Add one voxel for the volume boundary in order to ensure an empty space
  35.    // pad the x dimension
  36.    extent[0] -= (paddingImage[0] + 1);
  37.    std::cout << extent[0] << std::endl;
  38.    extent[1] += (paddingImage[0] + 1);
  39.    std::cout << extent[1] << std::endl;
  40.    // pad the y dimension
  41.    extent[2] -= (paddingImage[1] + 1);
  42.    std::cout << extent[2] << std::endl;
  43.    extent[3] += (paddingImage[1] + 1);
  44.    std::cout << extent[3] << std::endl;
  45.    // pad the z dimension
  46.    extent[4] -= (paddingImage[2] + 1);
  47.    std::cout << extent[4] << std::endl;
  48.    extent[5] += (paddingImage[2] + 1);
  49.    std::cout << extent[5] << std::endl;
  50.    padImageFilter->SetOutputWholeExtent(extent);
  51.    padImageFilter->Update();
  52.    std::cout << "pad image" << padImageFilter->GetOutput()->GetSpacing() << std::endl;
  53.  
  54.    // Convert to a floating point volume. 
  55.    // Note from proulxc: This is unfortunate, memory wise, but the VTK gaussian filter
  56.    // on integer types performs weirdly. It appears to apply the mask in three successive
  57.    // passes, rounding each time and causing non-symmetries and numerical errors in the
  58.    // process. 
  59.    const double MAX_VALUE = 1.0;
  60.    vtkSmartPointer<vtkImageCast> castImage =
  61.       vtkSmartPointer<vtkImageCast>::New();
  62.    castImage->SetInputConnection(padImageFilter->GetOutputPort());
  63.    castImage->SetOutputScalarTypeToFloat();
  64.    castImage->Update();
  65.    double space2[3];
  66.    castImage->GetOutput()->GetSpacing(space2);
  67.    
  68.    std::cout << "cast image" << space2[0] << std::endl;
  69.    vtkSmartPointer<vtkImageGaussianSmooth> gaussianSmoothingFilter =
  70.       vtkSmartPointer<vtkImageGaussianSmooth>::New();
  71.    gaussianSmoothingFilter->SetInputConnection(castImage->GetOutputPort());
  72.  
  73.    gaussianSmoothingFilter->SetDimensionality(3);
  74.  
  75.    double sigmaX = mConfiguration.mGaussianSmoothingFilterStdX;
  76.    double sigmaY = mConfiguration.mGaussianSmoothingFilterStdY;
  77.    double sigmaZ = mConfiguration.mGaussianSmoothingFilterStdZ;
  78.  
  79.    // Note that radius factor is multiplied by standard deviation to obtain
  80.    // the actual radius. Eg: sigma = 2 and radius factor = 2 -> actual radius = 4
  81.    // i.e. 9x9x9 kernel. We compute the radius factor so that the actual radius
  82.    // will always be the one set in the Algo Parameters config file, regardless of sigma.
  83.    double radiusFactorX = 1.0;
  84.    double radiusFactorY = 1.0;
  85.    double radiusFactorZ = 1.0;
  86.    double actualRadiusX = mConfiguration.mGaussianSmoothingActualRadiusX;
  87.    double actualRadiusY = mConfiguration.mGaussianSmoothingActualRadiusY;
  88.    double actualRadiusZ = mConfiguration.mGaussianSmoothingActualRadiusZ;
  89.    if (sigmaX >= std::numeric_limits<double>::epsilon())
  90.    {
  91.       radiusFactorX = actualRadiusX / sigmaX;
  92.    }
  93.    if (sigmaY >= std::numeric_limits<double>::epsilon())
  94.    {
  95.       radiusFactorY = actualRadiusY / sigmaY;
  96.    }
  97.    if (sigmaZ >= std::numeric_limits<double>::epsilon())
  98.    {
  99.       radiusFactorZ = actualRadiusZ / sigmaZ;
  100.    }
  101.  
  102.    gaussianSmoothingFilter->SetRadiusFactors(radiusFactorX, radiusFactorY, radiusFactorZ);
  103.    gaussianSmoothingFilter->SetStandardDeviations(sigmaX, sigmaY, sigmaZ);
  104.    gaussianSmoothingFilter->Update();
  105.    ORTHO_LOG_INFO("GaussianSmoothingImageOperation - End");
  106.    const auto mhdWriter = vtkSmartPointer<vtkMetaImageWriter>::New();
  107.    auto im = gaussianSmoothingFilter->GetOutput();
  108.    double spac[3];
  109.    im->GetSpacing(spac);
  110.    std::cout << "smoothed image" << spac[0] << std::endl;
  111.    mhdWriter->SetInputData(im);
  112.    mhdWriter->SetCompression(false);
  113.    mhdWriter->SetFileName("C:/temp/vtk9_smoothingwithdirMat.mhd");
  114.    mhdWriter->Write();
  115.    return im;
  116. }

Editor

You can edit this paste and save as new:


File Description
  • Smoothing operatiobn
  • Paste Code
  • 17 Feb-2023
  • 4.88 Kb
You can Share it: