Is there a complete example to use itk::CartesianToPolarTransform to transform an image?

ma1282029525 · August 28, 2025, 1:56am

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;
}

dzenanz · August 28, 2025, 3:14pm

Is this helpful?

github.com/KitwareMedical/SlicerITKUltrasound

ScanConvertCurvilinearArray/ScanConvertCurvilinearArray.cxx

master

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

This file has been truncated. show original

It heavily relies on this:

github.com/KitwareMedical/ITKUltrasound

include/itkCurvilinearArraySpecialCoordinatesImage.h

master

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkCurvilinearArraySpecialCoordinatesImage_h
#define itkCurvilinearArraySpecialCoordinatesImage_h

This file has been truncated. show original

ma1282029525 · August 29, 2025, 5:40am

I looked into this class and found that it converts from polar coordinates back to Cartesian coordinates. What I want is the opposite conversion. Can this class accomplish that? Also, are there no examples available? Is the itk::CartesianToPolarTransformclass a complete class? Why have all the methods I tried resulted in empty generated images?

ma1282029525 · August 29, 2025, 8:18am

It is working correctly; I had the transformation reversed.

dzenanz · August 29, 2025, 1:30pm

If you share your updated source code, it will be useful to people finding this discussion in the future.

ma1282029525 · September 1, 2025, 7:09am

ringartifact.cpp (7.7 KB)
Although the code is functioning correctly, the actual effectiveness in removing ring artifacts is unsatisfactory.