对dicom进行resize处理

1.通过pytdiocm处理

import pydicom
from pydicom.pixel_data_handlers.util import apply_voi_lut
from PIL import Image
import numpy as np
import os


# 重新设置dicom的大小(原始dicom长宽不一致时会有拉伸情况)
def resize_and_save_dicom(input_path, output_path, new_size):
    # Load the DICOM image
    ds = pydicom.dcmread(input_path)
    rows = ds.Rows
    columns = ds.Columns
    print('rows columns:', rows, columns)
    top, bottom, left, right = 0, 0, 0, 0
    if new_size[0] != new_size[1]:
        if new_size[0] > columns or new_size[1] > rows:
            print('can not resize...')
            return
        else:
            delta_lr = columns - new_size[0]
            delta_tb = rows - new_size[1]
            if delta_lr%2 == 0:
                left = (int)(delta_lr/2)
                right = left
            else:
                left = (int)(delta_lr/2)
                right = left + 1

            if delta_tb%2 == 0:
                top = (int)(delta_tb/2)
                bottom = top
            else:
                top = (int)(delta_tb/2)
                bottom = top + 1

    print('info:', top, bottom, left, right)
    # Apply VOI LUT if necessary
    if 'PixelData' in ds and ds.pixel_array is None:
        ds.pixel_array = apply_voi_lut(ds.pixel_array, ds)

    new_pixdata = ds.pixel_array[top:rows-bottom, left:columns-right]
    # Convert pixel_array to Pillow image
    image = Image.fromarray(new_pixdata)
    # Resize the image
    image = image.resize(new_size)
    print('image.size:', image.size)
    # Save the resized image as a DICOM file
    ds.PixelData = image.tobytes()
    ds.Columns, ds.Rows = image.size
    ds.save_as(output_path)


# Example usage
input_path = 'F:/1067_978.dcm'
output_path = 'F:/1067_978_1024.dcm'
new_size1 = (1024, 1024)

resize_and_save_dicom(input_path, output_path, new_size1)

folder_path = 'F:/1111'
for item in os.listdir(folder_path):
    input_path = folder_path + '/' + item
    output_path = folder_path + '/' + os.path.splitext(item)[0] + '_1790_1440.dcm'
    new_size1 = (1790, 1440)
    resize_and_save_dicom(input_path, output_path, new_size1)
    

2.通过itk进行resize,会先将dicom补位长宽一样

python实现

import numpy as np
import matplotlib.pyplot as plt
import itk


#对图像进行翻转
def flipImage(inImage, flipAxes):
    flipFilter = itk.FlipImageFilter[type(inImage)].New()
    flipFilter.SetInput(inImage)
    flipFilter.SetFlipAxes(flipAxes)
    flipFilter.Update()
    return flipFilter.GetOutput()


def readDicomFile(filepath, dataType=0, flipAxes = (False, False, False)):
    ImageType = itk.Image[itk.SS, 3]
    if dataType == 1:
        ImageType = itk.Image[itk.US, 3]
    ReaderType = itk.ImageFileReader[ImageType]
    reader = ReaderType.New()
    reader.SetFileName(filepath)
    reader.Update()
    if dataType == 0:
        castImageFilter = itk.CastImageFilter[type(reader.GetOutput()), itk.Image[itk.US, 3]].New()
        castImageFilter.SetInput(reader.GetOutput())
        castImageFilter.Update()
        return flipImage(castImageFilter.GetOutput(), flipAxes)
    return flipImage(reader.GetOutput(), flipAxes)


def show_image_array(image_array, save_name):
    plt.figure()
    plt.title('show image array:')
    if image_array.dtype == 'uint8':
        plt.imshow(image_array)
    else:
        plt.imshow(image_array, cmap='gray')
    plt.show()


dicom_path = 'F:/1067_978.dcm'
input_image = readDicomFile(dicom_path)
input_image_array = itk.GetArrayFromImage(input_image)[0].copy()
h, w = input_image_array.shape
print("image size: ", h, w, input_image_array.shape)
if h > w:
    pad_w = h-w
    pad_left = int(pad_w/2)
    pad_right = pad_w-pad_left
    print("pad left-right: ", pad_left, pad_right)
    input_image_array = np.pad(input_image_array, ((0, 0), (pad_left, pad_right)), "constant", constant_values=np.min(input_image_array))
elif w > h:
    pad_w = w-h
    pad_up = int(pad_w/2)
    pad_down = pad_w-pad_up
    print("pad left-right: ", pad_up, pad_down)
    input_image_array = np.pad(input_image_array, ((0, 0), (pad_up, pad_down)), "constant", constant_values=np.min(input_image_array))

h, w = input_image_array.shape
imOrigin = input_image.GetOrigin()
imRes = input_image.GetSpacing()
imRegion = input_image.GetBufferedRegion()
imSize = imRegion.GetSize()
new_image = itk.GetImageFromArray(input_image_array)
new_image.SetSpacing([imRes[0], imRes[1]])
new_image.SetOrigin([imOrigin[0], imOrigin[1]])
print(imRes)
scale = 1024/h
print(scale)
outputSize = [1024, 1024]
outputSpacing = [imRes[0] / scale, imRes[1] / scale]
print(outputSpacing)
outputOrigin = [imOrigin[0], imOrigin[1]]
transformType = itk.IdentityTransform[itk.D, 2].New()
FilterType = itk.ResampleImageFilter[type(new_image), type(new_image)]
filter = FilterType.New()
filter.SetInput(new_image)
filter.SetSize(outputSize)
filter.SetOutputSpacing(outputSpacing)
filter.SetOutputOrigin(outputOrigin)
filter.SetTransform(transformType)
filter.SetInterpolator(itk.NearestNeighborInterpolateImageFunction[type(new_image), itk.D].New())
filter.Update()
resample_image = filter.GetOutput()
show_image_array(itk.GetArrayFromImage(resample_image), '')

c++实现

itk::Image<unsigned short, 3>::Pointer vtkImageToItkImage(vtkImageData* inImage) {
	vtkNew<vtkImageFlip> flipFilter;
	flipFilter->SetFilteredAxis(1); // flip y axis
	flipFilter->SetInputData(inImage);
	flipFilter->Update();

	vtkImageData* flipImage = flipFilter->GetOutput();

	int* dims = flipImage->GetDimensions();
	double* spacings = flipImage->GetSpacing();
	double* origin = flipImage->GetOrigin();

	typedef itk::Image<unsigned short, 3> ResImageType;
	ResImageType::Pointer resImage = ResImageType::New();
	ResImageType::RegionType region;
	ResImageType::IndexType startindex;
	startindex.Fill(0);
	ResImageType::SizeType size;
	size[0] = dims[0];
	size[1] = dims[1];
	size[2] = dims[2];

	region.SetIndex(startindex);
	region.SetSize(size);

	resImage->SetOrigin(origin);
	resImage->SetSpacing(spacings);

	resImage->SetRegions(region);
	resImage->Allocate();
	resImage->FillBuffer(0);

	size_t sliceSize = size[0] * size[1]* size[2];
	memcpy((unsigned short*)resImage->GetBufferPointer(), (unsigned short*)flipImage->GetScalarPointer(), sliceSize * sizeof(unsigned short));
	return resImage;
}

//itkImage转为vtkImage
template<typename TInputImage>
void itkImageToVtkImage(TInputImage* inItkImage, vtkImageData* resData) {
	using FilterType = itk::ImageToVTKImageFilter<TInputImage>;
	auto filter = FilterType::New();
	filter->SetInput(inItkImage);
	try
	{
		filter->Update();
	}
	catch (const itk::ExceptionObject & error)
	{
		qDebug() << "itkImageToVtkImage error happen...";
		return;
	}
	vtkNew<vtkImageFlip> flipXFilter;
	flipXFilter->SetFilteredAxis(1); // flip y axis
	flipXFilter->SetInputData(filter->GetOutput());
	flipXFilter->Update();
	vtkImageData * myvtkImageData = flipXFilter->GetOutput();
	resData->DeepCopy(myvtkImageData);
	int* dims = resData->GetDimensions();
	qDebug() << "dims:" << dims[0] << dims[1] << dims[2];
}

void resizeVtkImageData(vtkImageData* inImage, vtkImageData* outImage, int newImageSize[3]) {
	if (!inImage || !outImage)
	{
		qDebug() << "inImage or outImage is nullptr...";
		return;
	}
	int* inDims = inImage->GetDimensions();
	if ( (inDims[0] == newImageSize[0]) && (inDims[1] == newImageSize[1]) && (inDims[2] == newImageSize[2]))
	{
		qDebug() << "no need resize...";
		outImage->DeepCopy(inImage);
		return;
	}
	itk::Image<unsigned short, 3>::Pointer inItkImage = vtkImageToItkImage(inImage);
	{
		using ImageType = itk::Image<unsigned short, 3>;

		float inImageOrigin[3] = { 0.0,0.0,0.0 };
		float inImageSpacings[3] = { 1.0,1.0,1.0 };
		inItkImage->SetSpacing(inImageSpacings);
		inItkImage->SetOrigin(inImageOrigin);
		const typename ImageType::SizeType resImageSize = inItkImage->GetBufferedRegion().GetSize();

		ImageType::Pointer newResImage = ImageType::New();
		bool addNewPixel = false;
		int newSize[2] = { 0 };
		if (resImageSize[0] != resImageSize[1])
		{
			qDebug() << "inequality:" << resImageSize[0] << resImageSize[1];
			addNewPixel = true;
			double minValue, maxValue;
			getImageMaxMinValue<ImageType>(inItkImage, minValue, maxValue);
			int topAdd, bottomAdd, leftAdd, rightAdd;
			topAdd = bottomAdd = leftAdd = rightAdd = 0;
			//处理长宽不同的情况
			bool wGreaterh = true;
			if (resImageSize[0] > resImageSize[1]) //长比宽大
			{
				wGreaterh = true;
				newSize[0] = newSize[1] = resImageSize[0];
				int delta = resImageSize[0] - resImageSize[1];
				if (delta % 2 == 0)
				{
					topAdd = delta / 2;
					bottomAdd = topAdd;
				}
				else {
					topAdd = delta / 2;
					bottomAdd = topAdd + 1;
				}
			}
			else {  //宽比长大
				wGreaterh = false;
				newSize[0] = newSize[1] = resImageSize[1];
				int delta = resImageSize[1] - resImageSize[0];
				if (delta % 2 == 0)
				{
					leftAdd = delta / 2;
					rightAdd = leftAdd;
				}
				else {
					leftAdd = delta / 2;
					rightAdd = leftAdd + 1;
				}
			}

			qDebug() << "new info:" << newSize[0] << newSize[1] << topAdd << bottomAdd << leftAdd << rightAdd << wGreaterh;
			ImageType::RegionType region;
			ImageType::IndexType startindex;
			startindex.Fill(0);
			int dims[3];
			dims[0] = newSize[0];
			dims[1] = newSize[1];
			dims[2] = 1;
			region.SetIndex(startindex);
			region.SetSize(0, dims[0]);
			region.SetSize(1, dims[1]);
			region.SetSize(2, dims[2]);

			newResImage->SetOrigin(inImageOrigin);
			newResImage->SetSpacing(inImageSpacings);

			newResImage->SetRegions(region);
			newResImage->Allocate();
			newResImage->FillBuffer(signed short(minValue));

			ImageType::IndexType tmpIndex;
			if (wGreaterh)
			{
				for (int j = topAdd; j < resImageSize[1] + topAdd; j++)
				{
					for (int i = 0; i < newSize[0]; i++) {
						tmpIndex[2] = 0;
						tmpIndex[1] = j - topAdd;
						tmpIndex[0] = i;
						signed short oldValue = inItkImage->GetPixel(tmpIndex);
						tmpIndex[2] = 0;
						tmpIndex[1] = j;
						tmpIndex[0] = i;
						newResImage->SetPixel(tmpIndex, oldValue);
					}
				}
			}
			else {
				for (int j = 0; j < newSize[1]; j++)
				{
					for (int i = leftAdd; i < resImageSize[0] + leftAdd; i++) {
						tmpIndex[2] = 0;
						tmpIndex[1] = j;
						tmpIndex[0] = i - leftAdd;
						signed short oldValue = inItkImage->GetPixel(tmpIndex);
						tmpIndex[2] = 0;
						tmpIndex[1] = j;
						tmpIndex[0] = i;
						newResImage->SetPixel(tmpIndex, oldValue);
					}
				}
			}
			newResImage->Update();
		}

		// Set the desired size of the output image
		ImageType::SizeType outputSize;
		outputSize[0] = newImageSize[0];  // Desired width
		outputSize[1] = newImageSize[1];  // Desired height
		outputSize[2] = 1;

		// Set the desired spacing of the output image
		ImageType::SpacingType outputSpacing;
		if (addNewPixel)
		{
			outputSpacing[0] = float(newSize[0]) / newImageSize[0];
			outputSpacing[1] = float(newSize[1]) / newImageSize[1];
		}
		else {
			outputSpacing[0] = float(resImageSize[0]) / newImageSize[0];
			outputSpacing[1] = float(resImageSize[1]) / newImageSize[1];
		}
		outputSpacing[2] = 1.0;  // Desired pixel spacing in the z direction

		// Set the desired origin of the output image
		ImageType::PointType outputOrigin;
		outputOrigin[0] = 0.0;  // Desired x coordinate of the origin
		outputOrigin[1] = 0.0;  // Desired y coordinate of the origin
		outputOrigin[2] = 0.0;  // Desired z coordinate of the origin

		// Set up the resampling filter
		typedef itk::ResampleImageFilter<ImageType, ImageType> ResampleFilterType;
		ResampleFilterType::Pointer resampler = ResampleFilterType::New();
		if (addNewPixel)
		{
			resampler->SetInput(newResImage);
		}
		else {
			resampler->SetInput(inItkImage);
		}
		itk::NearestNeighborInterpolateImageFunction<ImageType>::Pointer interpolator =
			itk::NearestNeighborInterpolateImageFunction<ImageType>::New();
		resampler->SetInterpolator(interpolator);
		resampler->SetSize(outputSize);
		resampler->SetOutputSpacing(outputSpacing);
		resampler->SetOutputOrigin(outputOrigin);
		// Perform the resampling
		resampler->Update();

		// Get the resized output image
		ImageType::Pointer outputImage = resampler->GetOutput();
		outputImage->SetSpacing(inImageSpacings);
		itkImageToVtkImage<itk::Image<unsigned short, 3>>(outputImage, outImage);
		return;
	}
}