ITK 序列图像处理

#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"

#include "itkImageFileReader.h"//读取头文件
#include "itkImageFileWriter.h"//写入头文件
#include "itkGDCMImageIO.h"//ImageIo头文件
#include "itkIntensityWindowingImageFilter.h"//调窗处理头文件
#include "itkBinaryThresholdImageFilter.h"//二值门限处理头文件
#include "itkBinaryBallStructuringElement.h"//基本球形
#include "itkBinaryErodeImageFilter.h"//腐蚀头文件
#include "itkBinaryDilateImageFilter.h"//膨胀头文件
#include "itkGrayscaleFillholeImageFilter.h"//灰度图像孔洞填充
#include "itkBinaryFillholeImageFilter.h"//二值图像孔洞填充
#include "itkRescaleIntensityImageFilter.h"//像素值映射头文件
#include "itkGradientMagnitudeImageFilter.h"//梯度强度头文件
#include "itkMultiplyImageFilter.h"//图像相乘
#include "itkVotingBinaryHoleFillingImageFilter.h"
#include "itkVotingBinaryIterativeHoleFillingImageFilter.h"
#include "itkRecursiveGaussianImageFilter.h"
#include "itkConnectedThresholdImageFilter.h"//区域生长
#include "itkImage.h"
#include "itkCastImageFilter.h"//类型转换
#include "itkCurvatureFlowImageFilter.h"//平滑处理
#include "itkShiftScaleImageFilter.h"//取反

int main(int argc, char* argv[])
{


    //选择像素类型和图像维度,生成读写图片类型
    typedef   float           InternalPixelType;
    const     unsigned int    Dimension = 3;
    typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;
    typedef itk::Image< signed short, Dimension >  InternalImageType1;
    
    typedef  float          OutputPixelType;
    const unsigned int      OutputDimension = 2;
    typedef itk::Image< OutputPixelType, OutputDimension > OutputImageType;


    //用图片类型作为参数实例化一个滤波器类型,通过调用实例化对象的 New()方式来创建的，并将结果指向itk::SmartPointers
    typedef  itk::ImageSeriesReader< InternalImageType > ReaderType;
    typedef  itk::ImageSeriesWriter<  InternalImageType1, OutputImageType  > WriterType;//自动在输入和输出之间转换

    ReaderType::Pointer reader = ReaderType::New();
    WriterType::Pointer writer = WriterType::New();
    //智能指针:


    //GDCMImageIO对象被作为读写滤波器使用的ImageIO的对象连接,找到合适的ImageIO类去执行IO操作的itk::ImageIOFactory机构
    typedef itk::GDCMImageIO             ImageIOType;
    ImageIOType::Pointer gdcmIO = ImageIOType::New();

    //文件生成器 GDCMSeriesFileNames将会搜索这个目录并生成DICOM文件名的序列
    typedef itk::GDCMSeriesFileNames     NamesGeneratorType;
    NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
    namesGenerator->SetInputDirectory("D:/leg_dcm/");//输入目录

    //文件目录(地址)通过GetFileNames()返回
    const ReaderType::FileNamesContainer& filenames = namesGenerator->GetInputFileNames();
    reader->SetImageIO(gdcmIO);
    reader->SetFileNames(filenames);


    //异常处理是C++语言中的一个标准部分 出现异常课跳出程序
    //异常处理中加入输出,可现实异常类型*
    try
    {
        reader->Update();
    }
    catch (itk::ExceptionObject& excp)
    {
        std::cerr << "Exception thrown while writing the image" << std::endl;
        std::cerr << excp << std::endl;
        std::cout << "succed1 " << std::endl;
        system("pause");
        return EXIT_FAILURE;
    }


    //调窗处理
    //1.窗宽窗位是CT图像特有的概念，MRI图像中可没有此概念
    //2.CT 图像必须先转换成 HU值再做窗宽窗位调整
    //3.nii格式的图片自动将灰度值转化为HU值;dcm格式的图不会将原始数据自动转化为Hu
    //4.CT值和灰度值属于不同的概念.CT值属于医学领域的概念，通常称亨氏单位(hounsfield unit ，HU)，反映了组织对X射线吸收程度
    //5.窗宽是CT图像上显示的CT值范围,窗位是窗的中心位置
    typedef itk::IntensityWindowingImageFilter<InternalImageType, InternalImageType >  IntensityFilterType;
    IntensityFilterType::Pointer intensityWindowing = IntensityFilterType::New();
    intensityWindowing->SetWindowMinimum(124);//最小窗值
    intensityWindowing->SetWindowMaximum(126);//最大窗值
    intensityWindowing->SetOutputMinimum(0);//
    intensityWindowing->SetOutputMaximum(255); // 
    intensityWindowing->SetInput(reader->GetOutput());
    intensityWindowing->Update();

    //二值门限
    typedef itk::BinaryThresholdImageFilter<InternalImageType, InternalImageType > ThresholdingFilterType;
    ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
    thresholder->SetInput(intensityWindowing->GetOutput());
    thresholder->SetUpperThreshold(255);//上下阈值参数设定
    thresholder->SetLowerThreshold(250);
    thresholder->SetOutsideValue(0);//输出值设定，背景为0，前景为255
    thresholder->SetInsideValue(255);
    thresholder->Update();


    //获取图像尺寸
    InternalImageType::Pointer image = thresholder->GetOutput();
    InternalImageType::RegionType region = image->GetLargestPossibleRegion();
    InternalImageType::SizeType size = region.GetSize();
    std::cout << size << std::endl;
    std::cout << size[0] << std::endl;//三个维度尺寸
    std::cout << size[1] << std::endl;
    std::cout << size[2] << std::endl;

    //创建掩膜
    using ImageType5 = itk::Image<float, 3>;
    ImageType5::Pointer image1 = ImageType5::New();

    ImageType5::IndexType start; //创建itk::Index对象,用来指定图像起点位置
    start[0] = 0;
    start[1] = 0;
    start[2] = 0;
    ImageType5::SizeType size1;   //创建itk::Size对象,指定图像各方向大小
    size1[0] = size[0];
    size1[1] = size[1];
    size1[2] = size[2];
    ImageType5::RegionType region1; //创建图像区域，并设置起点和大小
    region1.SetSize(size1);
    region1.SetIndex(start);
    image1->SetRegions(region1);
    ImageType5::SpacingType spacing;  //定义像素间距
    spacing[0] = 0.589844;
    spacing[1] = 0.589844;
    spacing[2] = 5.0;
    image1->SetSpacing(spacing);
    ImageType5::PointType origin;     //定义像素原点
    origin[0] = -151;
    origin[1] = -151;
    origin[2] = -325;
    image1->SetOrigin(origin);

    image1->Allocate();   //分配内存
    image1->FillBuffer(itk::NumericTraits<signed short>::Zero); //初始化图像缓冲区

    //对像素赋值
    for (int i = 0; i < size1[0]; i++) //像素所在行
    {
        start[0] = i;
        for (int k = 0; k < size1[2]; k++)
        {

            start[2] = k;
            for (int j = 0; j < size1[1] / 100 * 80; j++) //像素所在列
            {
                start[1] = j;
                image1->SetPixel(start, 1);  //给指定位置设置像素值

            }
            for (int j = size1[1] / 100 * 80; j < size1[1]; j++) //像素所在列
            {
                start[1] = j;
                image1->SetPixel(start, 0);  //给指定位置设置像素值
            }



        }
    }
    for (int i = 0; i < size1[0]; i++) //像素所在行
    {
        start[0] = i;
        for (int k = 0; k < size1[2] - 78; k++)
        {

            start[2] = k;
            for (int j = 0; j < size1[1] / 100 * 84; j++) //像素所在列
            {
                start[1] = j;
                image1->SetPixel(start, 1);  //给指定位置设置像素值

            }
            for (int j = size1[1] / 100 * 84; j < size1[1]; j++) //像素所在列
            {
                start[1] = j;
                image1->SetPixel(start, 0);  //给指定位置设置像素值
            }



        }
    }

    //掩膜处理,提取感兴趣区域
    typedef itk::MultiplyImageFilter<InternalImageType, ImageType5, InternalImageType >FilterType2;
    FilterType2::Pointer multfilter = FilterType2::New();
    multfilter->SetInput1(image);
    multfilter->SetInput2(image1);
    try
    {
        multfilter->Update();
    }
    catch (itk::ExceptionObject& excep)
    {
        std::cerr << "Exception caught !" << std::endl;
        std::cerr << excep << std::endl;
        std::cout << "succed " << std::endl;
        system("pause");
        return EXIT_FAILURE;
    }

    //01取反
    typedef itk::ShiftScaleImageFilter< InternalImageType, InternalImageType >ShiftScaleType0;
    ShiftScaleType0::Pointer shiftScale_mask = ShiftScaleType0::New();
    shiftScale_mask->SetInput(multfilter->GetOutput());
    shiftScale_mask->SetShift(-1);
    shiftScale_mask->SetScale(-1);
    shiftScale_mask->Update();

        
    //平滑处理
    typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType > CurvatureFlowImageFilterType;
    CurvatureFlowImageFilterType::Pointer smoothing = CurvatureFlowImageFilterType::New();
    smoothing->SetInput(shiftScale_mask->GetOutput());
    smoothing->SetNumberOfIterations(6);
    smoothing->SetTimeStep(0.625);
    smoothing->Update();


    //区域生长法
    typedef itk::ConnectedThresholdImageFilter< InternalImageType, InternalImageType > ConnectedFilterType;
    ConnectedFilterType::Pointer connectedThreshold = ConnectedFilterType::New();
    const InternalPixelType lowerThreshold = -1;//提取区域上下门限阈值
    const InternalPixelType upperThreshold = 1;
    connectedThreshold->SetLower(lowerThreshold);
    connectedThreshold->SetUpper(upperThreshold);
    connectedThreshold->SetReplaceValue(255);//设置提取区域像素值,其余区域默认为0
    InternalImageType::IndexType  index;//选取种子点
    index[0] = 0;
    index[1] = 0;
    index[2] = 0;
    connectedThreshold->SetSeed(index);
    connectedThreshold->SetInput(smoothing->GetOutput());
    connectedThreshold->Update();


    //01取反
    typedef itk::ShiftScaleImageFilter< InternalImageType, InternalImageType >ShiftScaleType0;
    ShiftScaleType0::Pointer shiftScale_mask1 = ShiftScaleType0::New();
    shiftScale_mask1->SetInput(connectedThreshold->GetOutput());
    shiftScale_mask1->SetShift(-1);
    shiftScale_mask1->SetScale(-1);
    shiftScale_mask1->Update();


    //形态学处理
    //用于二值图像的构造成员
    typedef itk::BinaryBallStructuringElement< InternalPixelType, 3  > StructuringElementType;
    StructuringElementType  structuringElement;
    structuringElement.SetRadius(1);   //领域大小为3*3*3
    structuringElement.CreateStructuringElement();

    typedef itk::BinaryErodeImageFilter <InternalImageType, InternalImageType, StructuringElementType >  ErodeFilterType;// 腐蚀，需要用输入、输出图像类型和构造成员实例化
    ErodeFilterType::Pointer  binaryErode = ErodeFilterType::New();
    binaryErode->SetInput(shiftScale_mask1->GetOutput());
    binaryErode->SetKernel(structuringElement);
    binaryErode->SetErodeValue(255);
    binaryErode->Update();


    //类型转换
    typedef itk::CastImageFilter< InternalImageType, InternalImageType1 > CastingFilterType;
    CastingFilterType::Pointer caster = CastingFilterType::New();
    caster->SetInput(binaryErode->GetOutput());
    caster->Update();


    //输出文件夹
    const char* outputDirectory = "D:/sus7/";
    itksys::SystemTools::MakeDirectory(outputDirectory);
    //设置输出文件夹
    namesGenerator->SetOutputDirectory(outputDirectory);
    writer->SetFileNames(namesGenerator->GetOutputFileNames());
    //用原来的文件头
    writer->SetMetaDataDictionaryArray(reader->GetMetaDataDictionaryArray());
    writer->SetInput(caster->GetOutput());
    writer->SetImageIO(gdcmIO);
    try
    {
        writer->Update();
        //writer->UpdateLargestPossibleRegion();
    }
    catch (itk::ExceptionObject& excep)
    {
        std::cerr << "Exception caught !" << std::endl;
        std::cerr << excep << std::endl;
        std::cout << "succed " << std::endl;
        system("pause");
        return EXIT_FAILURE;
    }

    system("pause");
    return EXIT_SUCCESS;
}