I want to complete the task of removing ring artifacts. I also came across the FourierStripeArtifactImageFilter class, but I failed at the previous step—I am unable to convert the Cartesian coordinate system image to a polar coordinate image. Is there any available example for reference? This is my code.
#include “itkFourierStripeArtifactImageFilter.h”
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkCartesianToPolarTransform.h"
#include "itkPolarToCartesianTransform.h"
#include "itkRescaleIntensityImageFilter.h"
#include <iostream>
#include <cmath>
#include <random>
void CreateTestPatternImage()
{
try
{
// 定义图像类型
constexpr unsigned int Dimension = 2;
using PixelType = unsigned char;
using ImageType = itk::Image<PixelType, Dimension>;
// 创建图像
ImageType::Pointer image = ImageType::New();
// 设置图像大小和原点
ImageType::SizeType size;
size[0] = 512;
size[1] = 512;
ImageType::IndexType start;
start.Fill(0);
ImageType::RegionType region;
region.SetSize(size);
region.SetIndex(start);
image->SetRegions(region);
image->Allocate();
image->FillBuffer(0);
// 设置图像间距和原点
ImageType::SpacingType spacing;
spacing[0] = 1.0;
spacing[1] = 1.0;
image->SetSpacing(spacing);
ImageType::PointType origin;
origin[0] = 0.0;
origin[1] = 0.0;
image->SetOrigin(origin);
// 计算图像中心
ImageType::IndexType centerIndex;
centerIndex[0] = size[0] / 2;
centerIndex[1] = size[1] / 2;
// 创建多个不同粗细的同心圆环
const int numRings = 30;
const int minRingWidth = 1; // 最小圆环宽度改为1像素
const int maxRingWidth = 1; // 最大圆环宽度改为3像素
// 随机数生成器
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_int_distribution<> widthDist(minRingWidth, maxRingWidth);
// 存储圆环半径和宽度
std::vector<std::pair<double, double>> rings;
double currentRadius = 30.0; // 起始半径
for (int i = 0; i < numRings; ++i)
{
int ringWidth = widthDist(gen);
rings.push_back({ currentRadius, currentRadius + ringWidth });
currentRadius += ringWidth + 2; // 减小圆环之间的间隔
}
// 遍历所有像素
itk::ImageRegionIterator<ImageType> imageIterator(image, region);
while (!imageIterator.IsAtEnd())
{
// 获取当前像素索引
ImageType::IndexType index = imageIterator.GetIndex();
// 计算到中心的距离
double dx = index[0] - centerIndex[0];
double dy = index[1] - centerIndex[1];
double distance = std::sqrt(dx * dx + dy * dy);
// 绘制圆环
for (const auto& ring : rings)
{
if (distance >= ring.first && distance < ring.second)
{
imageIterator.Set(255); // 白色圆环
break;
}
}
++imageIterator;
}
// 添加矩形
ImageType::IndexType rectStart;
rectStart[0] = 100;
rectStart[1] = 100;
ImageType::SizeType rectSize;
rectSize[0] = 80;
rectSize[1] = 50;
ImageType::RegionType rectRegion;
rectRegion.SetIndex(rectStart);
rectRegion.SetSize(rectSize);
itk::ImageRegionIterator<ImageType> rectIterator(image, rectRegion);
while (!rectIterator.IsAtEnd())
{
rectIterator.Set(200); // 灰色矩形
++rectIterator;
}
// 添加小圆形
ImageType::IndexType circleCenter;
circleCenter[0] = 400;
circleCenter[1] = 400;
double circleRadius = 25.0;
ImageType::RegionType circleRegion;
ImageType::IndexType circleStart;
circleStart[0] = circleCenter[0] - circleRadius - 1;
circleStart[1] = circleCenter[1] - circleRadius - 1;
ImageType::SizeType circleSize;
circleSize[0] = 2 * circleRadius + 2;
circleSize[1] = 2 * circleRadius + 2;
circleRegion.SetIndex(circleStart);
circleRegion.SetSize(circleSize);
itk::ImageRegionIterator<ImageType> circleIterator(image, circleRegion);
while (!circleIterator.IsAtEnd())
{
ImageType::IndexType idx = circleIterator.GetIndex();
double dx = idx[0] - circleCenter[0];
double dy = idx[1] - circleCenter[1];
double distance = std::sqrt(dx * dx + dy * dy);
if (distance <= circleRadius)
{
circleIterator.Set(150); // 浅灰色圆形
}
++circleIterator;
}
// 保存图像
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("ring.mha");
writer->SetInput(image);
writer->Update();
std::cout << "测试图案图像已保存为: ring.mha" << std::endl;
}
catch (const itk::ExceptionObject& error)
{
std::cerr << "创建或保存测试图案图像时出错: " << error << std::endl;
throw;
}
}
int main(int argc, char* argv[])
{
try
{
if (argc < 3)
{
std::cerr << "Usage: " << argv[0] << " inputImage outputImage" << std::endl;
return EXIT_FAILURE;
}
CreateTestPatternImage();
constexpr unsigned int Dimension = 2;
using PixelType = float;
using ImageType = itk::Image< PixelType, Dimension >;
// Read input image
using ReaderType = itk::ImageFileReader< ImageType >;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(argv[1]);
reader->Update();
ImageType::Pointer inputImage = reader->GetOutput();
// 保存原始图像（转换为unsigned char以便查看）
using RescaleFilterType = itk::RescaleIntensityImageFilter<ImageType, ImageType>;
RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
rescaleFilter->SetInput(inputImage);
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("original_image.mha");
writer->SetInput(rescaleFilter->GetOutput());
writer->Update();
// Get image properties
ImageType::RegionType region = inputImage->GetLargestPossibleRegion();
ImageType::SizeType size = region.GetSize();
ImageType::SpacingType spacing = inputImage->GetSpacing();
ImageType::PointType origin = inputImage->GetOrigin();
// Set transform center to image center
ImageType::IndexType centerIndex;
centerIndex[0] = size[0] / 2;
centerIndex[1] = size[1] / 2;
ImageType::PointType centerPoint;
inputImage->TransformIndexToPhysicalPoint(centerIndex, centerPoint);
// Convert to polar coordinates
using PolarTransformType = itk::CartesianToPolarTransform< double, Dimension >;
PolarTransformType::Pointer polarTransform = PolarTransformType::New();
polarTransform->SetCenter(centerPoint);
using ResampleFilterType = itk::ResampleImageFilter< ImageType, ImageType >;
ResampleFilterType::Pointer polarResampler = ResampleFilterType::New();
polarResampler->SetTransform(polarTransform);
polarResampler->SetInput(inputImage);
// Set output parameters for polar image
ImageType::SizeType polarSize;
polarSize[0] = 360; // Angular resolution
polarSize[1] = std::sqrt(centerIndex[0] * centerIndex[0] + centerIndex[1] * centerIndex[1]); // Radial resolution
// 设置极坐标图像的原点和间距
ImageType::PointType polarOrigin;
polarOrigin[0] = 0.0; // 角度从0开始
polarOrigin[1] = 0.0; // 半径从0开始
ImageType::SpacingType polarSpacing;
polarSpacing[0] = 2 * itk::Math::pi / polarSize[0]; // 角度间距
polarSpacing[1] = 1.0; // 半径间距
using LinearInterpolatorType = itk::LinearInterpolateImageFunction<ImageType, ImageType>;
auto interpolator = LinearInterpolatorType::New();
polarResampler->SetSize(polarSize);
polarResampler->SetOutputOrigin(polarOrigin);
polarResampler->SetOutputSpacing(polarSpacing);
polarResampler->SetOutputDirection(inputImage->GetDirection());
polarResampler->SetInterpolator(itk::LinearInterpolateImageFunction< ImageType, double >::New());
polarResampler->Update();
ImageType::Pointer polarImage = polarResampler->GetOutput();
ImageType::IndexType testIndex;
testIndex[0] = 0;
testIndex[1] = 30;
std::cout << "Polar image at (0,30): " << polarImage->GetPixel(testIndex) << std::endl;
testIndex[0] = 90;
testIndex[1] = 30;
std::cout << "Polar image at (90,30): " << polarImage->GetPixel(testIndex) << std::endl;
testIndex[0] = 180;
testIndex[1] = 30;
std::cout << "Polar image at (180,30): " << polarImage->GetPixel(testIndex) << std::endl;
testIndex[0] = 0;
testIndex[1] = 100;
std::cout << "Polar image at (0,100): " << polarImage->GetPixel(testIndex) << std::endl;
// 检查输入图像是否有效
std::cout << "Input image size: " << size[0] << " x " << size[1] << std::endl;
std::cout << "Input image spacing: " << spacing[0] << " x " << spacing[1] << std::endl;
std::cout << "Center point: " << centerPoint[0] << ", " << centerPoint[1] << std::endl;
std::cout << "Polar image size: " << polarSize[0] << " x " << polarSize[1] << std::endl;
std::cout << "Polar image spacing: " << polarSpacing[0] << " x " << polarSpacing[1] << std::endl;
WriterType::Pointer polarWriter = WriterType::New();
polarWriter->SetFileName("polar_image.mha");
polarWriter->SetInput(polarResampler->GetOutput());
polarWriter->Update();
// Apply Fourier stripe artifact filter
using FourierFilterType = itk::FourierStripeArtifactImageFilter< ImageType >;
FourierFilterType::Pointer fourierFilter = FourierFilterType::New();
fourierFilter->SetInput(polarResampler->GetOutput());
fourierFilter->SetDirection(0); // Filter along angular direction
fourierFilter->Update();
ImageType::Pointer filteredPolarImage = fourierFilter->GetOutput();
// 保存滤波后的极坐标图像
RescaleFilterType::Pointer filteredPolarRescale = RescaleFilterType::New();
filteredPolarRescale->SetInput(filteredPolarImage);
filteredPolarRescale->SetOutputMinimum(0);
filteredPolarRescale->SetOutputMaximum(255);
WriterType::Pointer filteredPolarWriter = WriterType::New();
filteredPolarWriter->SetFileName("filtered_polar_image.mha");
filteredPolarWriter->SetInput(filteredPolarRescale->GetOutput());
filteredPolarWriter->Update();
// Convert back to Cartesian coordinates
using CartesianTransformType = itk::PolarToCartesianTransform< double, Dimension >;
CartesianTransformType::Pointer cartesianTransform = CartesianTransformType::New();
cartesianTransform->SetCenter(centerPoint);
ResampleFilterType::Pointer cartesianResampler = ResampleFilterType::New();
cartesianResampler->SetTransform(cartesianTransform);
cartesianResampler->SetInput(filteredPolarImage);
// Set output parameters to match original image
cartesianResampler->SetSize(inputImage->GetLargestPossibleRegion().GetSize());
cartesianResampler->SetOutputOrigin(inputImage->GetOrigin());
cartesianResampler->SetOutputSpacing(inputImage->GetSpacing());
cartesianResampler->SetOutputDirection(inputImage->GetDirection());
cartesianResampler->SetInterpolator(itk::LinearInterpolateImageFunction< ImageType, double >::New());
cartesianResampler->Update();
// Write output image
WriterType::Pointer outputWriter = WriterType::New();
outputWriter->SetFileName(argv[2]);
outputWriter->SetInput(cartesianResampler->GetOutput());
outputWriter->Update();
std::cout << "处理完成，输出已保存到: " << argv[2] << std::endl;
std::cout << "中间图像已保存: original_image.mha, polar_image.mha, filtered_polar_image.mha" << std::endl;
}
catch (const itk::ExceptionObject& error)
{
std::cerr << "ITK异常: " << error << std::endl;
return EXIT_FAILURE;
}
catch (const std::exception& error)
{
std::cerr << "标准异常: " << error.what() << std::endl;
return EXIT_FAILURE;
}
catch (...)
{
std::cerr << "未知异常发生" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}