ITK学习笔记——将处理得到的二维掩码输出为连续序列

本文主要是介绍ITK学习笔记——将处理得到的二维掩码输出为连续序列，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

我们用ITK对图像进行处理的时候往往先会将16位的原始图像转换为0～255的8位无符号整数，但是要想输出为原始连续的序列，还要将这些掩码结果重新转换为16位的原始图像格式。这个过程用ITK实现起来比较复杂，我会一步步进行讲解。

一、读取dicom图像
将读取的signed short图像归一化到0～255的unsigned char图像

ImageType::Pointer readdicom(string filename) {ReaderType::Pointer reader = ReaderType::New();reader->SetFileName(filename);                                  //输入dicom图像ImageIOType::Pointer gdcmImageIO = ImageIOType::New();reader->SetImageIO(gdcmImageIO);try{reader->Update();}catch (itk::ExceptionObject & e){std::cerr << "exception in file reader " << std::endl;std::cerr << e << std::endl;}RescaleFilterType::Pointer rescaler = RescaleFilterType::New();//0～255灰度图rescaler->SetOutputMinimum(0);rescaler->SetOutputMaximum(255);rescaler->SetInput(reader->GetOutput());                 //对读入的dicm图片进行RescaleIntensityImageFilter处理rescaler->Update();//获取图像 ImageType::Pointer image = rescaler->GetOutput();return image;
}

二、处理原图像
假设我们将原图利用OSTU进行二值化处理

int ostu(ImageType::Pointer image)
{int width = image->GetLargestPossibleRegion().GetSize()[0];int heigth = image->GetLargestPossibleRegion().GetSize()[1];int x = 0, y = 0;int pixelCount[256] = { 0 };//每个像素的计数float pixelPro[256] = { 0 };;    //每种像素比例int i, j, pixelSum = width * heigth, threshold = 0;//统计灰度级中每个像素在整幅图像中的个数for (int i = 0; i < width; i++)for (int j = 0; j < heigth; j++){ImageType::IndexType pixelIndex;pixelIndex[0] = i;pixelIndex[1] = j;pixelCount[image->GetPixel(pixelIndex)]++;}//计算每个像素在整幅图像中的比例for (int i = 0; i < 256; i++){pixelPro[i] = (float)pixelCount[i] / (float)pixelSum;}//经典ostu算法,得到前景和背景的分割//遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;for (i = 0; i < 256; i++){w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;for (j = 0; j < 256; j++){if (j <= i) //背景部分{//以i为阈值分类，第一类总的概率w0 += pixelPro[j];u0tmp += j * pixelPro[j];}else       //前景部分{//以i为阈值分类，第二类总的概率w1 += pixelPro[j];u1tmp += j * pixelPro[j];}}u0 = u0tmp / w0;		//第一类的平均灰度u1 = u1tmp / w1;		//第二类的平均灰度u = u0tmp + u1tmp;		//整幅图像的平均灰度//计算类间方差deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);//找出最大类间方差以及对应的阈值if (deltaTmp > deltaMax){deltaMax = deltaTmp;threshold = i;}}ItType it(image, image->GetRequestedRegion());//将迭代器移动到首个元素 it.GoToBegin();//遍历像素，直至结束 while (!it.IsAtEnd()){//获取像素值 ImageType::PixelType value = it.Get();if ((int)value > threshold){it.Set(255);}else{it.Set(0);}//迭代器移动至下一元素 ++it;}//返回最佳阈值;return threshold;

三、输出位原始格式
首先对每张图片进行二值化处理，并与原始图像叠加。
voliter就是关键的三维数据。

void SegmentSingleDCM(int sliceNumber, std::string inputFile, VolumeIteratorType& voliter) {ImageType::Pointer image = readdicom(inputFile);   //将读取的signed short图像归一化到0～255的unsigned char图像int threshold = ostu(image);std::cout << "best_threshold = " << threshold << std::endl;ReaderType::Pointer reader = ReaderType::New();reader->SetFileName(inputFile);                                  //输入dicom图像ImageIOType::Pointer gdcmImageIO = ImageIOType::New();reader->SetImageIO(gdcmImageIO);try{reader->Update();}catch (itk::ExceptionObject & e){std::cerr << "exception in file reader " << std::endl;std::cerr << e << std::endl;}//读取原始图像InputImageType *input = reader->GetOutput();DCMIteratorType input_iter(input, input->GetRequestedRegion());input_iter.GoToBegin();//获取肺实质掩膜ImageType *bin = image;ItType mask_initer(bin, bin->GetRequestedRegion());mask_initer.GoToBegin();//原始图像与肺实质掩膜叠加while ((!input_iter.IsAtEnd()) || (!mask_initer.IsAtEnd())){ImageType::PixelType value_mask_lung = mask_initer.Get();if ((int(value_mask_lung)) == 255){voliter.Set(input_iter.Get());}else{voliter.Set(-1024);}//迭代器移动至下一元素 ++input_iter;++mask_initer;++voliter;}
}

对每张切片进行循环处理，并写入到voliter三维数据中，就能输出连续的和原始图像一样格式的序列了。

void lung_vessel_seg(string directory) {string inputdirectory = "D:/input/";//读入dicom序列图片NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();namesGenerator->SetInputDirectory(inputdirectory);//输入CT_dicm图像系列所在目录const ReaderType_series::FileNamesContainer & filenames =namesGenerator->GetInputFileNames();unsigned int numberOfFilenames = filenames.size();   //文件大小std::cout << numberOfFilenames << std::endl;ImageIOType::Pointer gdcmIO = ImageIOType::New();ReaderType_series::Pointer reader = ReaderType_series::New();reader->SetImageIO(gdcmIO);reader->SetFileNames(filenames);try{reader->Update();}catch (itk::ExceptionObject &excp){std::cerr << "Exception thrown while writing the image" << std::endl;std::cerr << excp << std::endl;}VolumeImageType::Pointer reconstruct = reader->GetOutput();reconstruct->Allocate();VolumeIteratorType voliter(reconstruct, reconstruct->GetRequestedRegion());voliter.GoToBegin();for (int fni = 0; fni < numberOfFilenames; fni++){std::cout << "filename # " << fni << " = ";std::cout << filenames[fni] << std::endl;SegmentSingleDCM(fni, filenames[fni], voliter);//==========================================================}reconstruct->Update();//写进dicom序列string outputDirectory = "D:/output/"itksys::SystemTools::MakeDirectory(outputDirectory);typedef signed short    OutputPixelType;const unsigned int      OutputDimension = 2;typedef itk::Image< OutputPixelType, OutputDimension >    Image2DType;typedef itk::ImageSeriesWriter<VolumeImageType, Image2DType >  SeriesWriterType;SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();seriesWriter->SetInput(reconstruct);seriesWriter->SetImageIO(gdcmIO);namesGenerator->SetOutputDirectory(outputDirectory);seriesWriter->SetFileNames(namesGenerator->GetOutputFileNames());seriesWriter->SetMetaDataDictionaryArray(reader->GetMetaDataDictionaryArray());seriesWriter->Update();
}