JAVA 版本 VTK 切割DICOM体数据 任意层厚输出成切片文件

在此之前，需要准备的资料大约有:
1、Java 编译好的VTK版本，请自行搜索Camke 编译VTK有很多教程。有条件的，英文好的请直接梯子去VTK官网搜，有原版的编辑过程。我本人用的是7.1.1版本的VTK。
2、dcm4che-core 5.12.0依赖包。

开发环境 jdk1.8 64位，eclipse。

废话不多说，直接开始正题。代码本人亲测，可以切割指定层厚。

import java.io.File;

import org.dcm4che3.data.Attributes;
import org.dcm4che3.data.Implementation;
import org.dcm4che3.data.Tag;
import org.dcm4che3.data.VR;
import org.dcm4che3.io.DicomEncodingOptions;
import org.dcm4che3.io.DicomOutputStream;
import org.dcm4che3.util.UIDUtils;

import vtk.vtkAlgorithmOutput;
import vtk.vtkDICOMImageReader;
import vtk.vtkImageData;
import vtk.vtkImageReslice;
import vtk.vtkMatrix4x4;
import vtk.vtkNativeLibrary;

/**
 * 切割机
 * 
 * @author Jmx 2020-09-03
 *
 */
public class NewMPRReslicer {

	// 需要静态引用vtk的dll类库等,否则会出错: java.lang.UnsatisfiedLinkLOG.error: VTKInit()
	static {
		if (!vtkNativeLibrary.LoadAllNativeLibraries()) {
			for (vtkNativeLibrary lib : vtkNativeLibrary.values()) {
				if (!lib.IsLoaded()) {
					System.err.println(lib.GetLibraryName() + " not loaded");
				}
			}
		}
		//只是隐藏了在VTK错误或警告出现时可能会弹出的任何调试信息。您还可以提供一个文件路径，以便将任何错误写入磁盘。
		vtkNativeLibrary.DisableOutputWindow(null);
	}

	/**
	 * 切割横断位
	 * seriesDIR 要切割的序列文件夹
	 * outDir 输出的文件夹
	 * thickness 要输出的层厚，不写默认就是dicom序列原始层厚
	 */
	public static void excute(String seriesDIR, String outDir, Double thickness) {
		//加载体数据
		vtkDICOMImageReader imageReader = new vtkDICOMImageReader();
		imageReader.SetDirectoryName(seriesDIR);
		imageReader.Update();
		vtkImageData imageData = imageReader.GetOutput();
		//加载体数据结束

		double bounds[] = imageData.GetBounds();//获取体数据的边界
		double spacing[] = imageData.GetSpacing();//图像的像素间距 
		int extent[] = imageData.GetExtent();//图像的宽高和图像数量
		double origin[] = imageData.GetOrigin();//图像的坐标原点

		// 找到中心点位置(贯穿点)
		double center[] = new double[3];
		center[0] = origin[0] + spacing[0] * 0.5 * (extent[0] + extent[1]);// x
		center[1] = origin[1] + spacing[1] * 0.5 * (extent[2] + extent[3]);// y
		center[2] = origin[2] + spacing[2] * 0.5 * (extent[4] + extent[5]);// z

		// 通过修改切面和贯穿点，可以得到需要的重建图层(冠状、矢状、横断、斜切等)
		int columneIndex = 2;//横断位
		double[] imageOrientationPatient = new double[] { 1, 0, 0, 0, 1, 0 };//
		vtkMatrix4x4 resliceAxes = new vtkMatrix4x4();
		// 此矩阵为横断面(方式1)
		//	double axialElements[] = { 
		//	1, 0, 0, 0, 
		//	0, 1, 0, 0, 
		//	0, 0, 1, 0, 
		//	0, 0, 0, 1
		//	};
		//	resliceAxes.DeepCopy(axialElements);
		//方式2
		resliceAxes.SetElement(0, 0, 1);
		resliceAxes.SetElement(0, 1, 0);
		resliceAxes.SetElement(0, 2, 0);

		resliceAxes.SetElement(1, 0, 0);
		resliceAxes.SetElement(1, 1, 1);
		resliceAxes.SetElement(1, 2, 0);

		resliceAxes.SetElement(2, 0, 0);
		resliceAxes.SetElement(2, 1, 0);
		resliceAxes.SetElement(2, 2, 1);

		if (null != thickness && thickness > 0) {
			spacing[2] = thickness;//以用户传的层厚为准
		}
		int instanceCount = 0;
		//通过层厚计算dicom输出的数量
		double count = (bounds[5] - bounds[4]) / spacing[2] + 1;
		if (count - (int) count != 0) {
			//说明除不尽
			instanceCount = (int) count + 1;
		} else {
			instanceCount = (int) count;
		}

		// 设置中心点
		resliceAxes.SetElement(0, 3, center[0]);
		resliceAxes.SetElement(1, 3, center[1]);
		resliceAxes.SetElement(2, 3, center[2]);
		// 图像重取filter(从体数据中)
		vtkImageReslice imageReslice = new vtkImageReslice();
		vtkAlgorithmOutput algorithmOutput = imageReader.GetOutputPort();
		imageReslice.SetInputConnection(algorithmOutput);
		imageReslice.SetOutputDimensionality(2);// 输出为2D数据
		imageReslice.SetResliceAxes(resliceAxes);
		imageReslice.SetOutputSpacing(spacing);
		//设置插值法
		imageReslice.SetInterpolationModeToLinear();// 双线性插值法
		// 此处从中心点预切割一幅,目的为获取各种参数
		imageReslice.Update();
		//=======================重建结束，准备输出图像======

		vtkImageData imageResliceData = imageReslice.GetOutput();
		extent = imageResliceData.GetExtent();
		spacing = imageResliceData.GetSpacing();
		origin = imageResliceData.GetOrigin();
		bounds = imageResliceData.GetBounds();//获取体数据的边界
		int volume_x = extent[1] + 1;
		int volume_y = extent[3] + 1;

		Attributes fmi = new Attributes(6);
		fmi.setBytes(Tag.FileMetaInformationVersion, VR.OB, new byte[] { 0, 1 });
		fmi.setString(Tag.MediaStorageSOPClassUID, VR.UI, "1.2.840.10008.5.1.4.1.1.128");//这里要动态获取dicom的值
		fmi.setString(Tag.MediaStorageSOPInstanceUID, VR.UI, UIDUtils.createUID("2.25"));
		fmi.setString(Tag.TransferSyntaxUID, VR.UI, "1.2.840.10008.1.2.1");
		fmi.setString(Tag.ImplementationClassUID, VR.UI, Implementation.getClassUID());
		fmi.setString(Tag.ImplementationVersionName, VR.SH, Implementation.getVersionName());

		Attributes data = new Attributes();
		//图像必要的一些tag
		String studuUID = UIDUtils.createUID("1.1");
		String seriesUID = UIDUtils.createUID("1.2");
		//0020 Group
		data.setString(Tag.StudyInstanceUID, VR.UI, studuUID);
		data.setString(Tag.SeriesInstanceUID, VR.UI, seriesUID);

		//0028 Group
		data.setString(Tag.SamplesPerPixel, VR.US, "1");
		data.setString(Tag.PhotometricInterpretation, VR.CS, "MONOCHROME2");
		data.setInt(Tag.Rows, VR.IS, volume_x);
		data.setInt(Tag.Columns, VR.IS, volume_y);
		data.setInt(Tag.BitsAllocated, VR.IS, 16);
		data.setInt(Tag.BitsStored, VR.IS, 16);
		data.setInt(Tag.HighBit, VR.IS, 15);
		data.setInt(Tag.PixelRepresentation, VR.IS, 0);

		data.setInt(Tag.Rows, VR.US, volume_y);//必要的Tag
		data.setInt(Tag.Columns, VR.US, volume_x);//必要的Tag
		data.setInt(Tag.NumberOfSlices, VR.US, instanceCount);
		data.setDouble(Tag.ImageOrientationPatient, VR.DS, imageOrientationPatient);

		data.setDouble(Tag.PixelSpacing, VR.DS, spacing[0], spacing[1]);
		data.setDouble(Tag.SliceThickness, VR.DS, spacing[2]);

		double[] imagePositionPatient = new double[3];
		imagePositionPatient[0] = origin[0];
		imagePositionPatient[1] = origin[1];

		new File(outDir).mkdirs();
		for (int i = 0; i < instanceCount; i++) {
			double sliceLocation = i * spacing[2];
			resliceAxes.SetElement(columneIndex, 3, sliceLocation);//设置第4个切面，来切要的体数据
			imageReslice.SetResliceAxes(resliceAxes);
			imageReslice.Update();
			vtkImageData vtkImageData = imageReslice.GetOutput();
			data.setString(Tag.SeriesDescription, VR.LO, "JMX SC AXIL");
			imagePositionPatient[2] = sliceLocation;
			data.setDouble(Tag.ImagePositionPatient, VR.DS, imagePositionPatient);
			data.setDouble(Tag.SliceLocation, VR.DS, sliceLocation);
			data.setInt(Tag.InstanceNumber, VR.IS, i + 1);
			byte[] dicomarray = new byte[volume_x * volume_y * 2];
			int dicomindex = 0;
			for (int y = 0; y < volume_y; y++) {
				for (int x = 0; x < volume_x; x++) {
					double pixelvalue = vtkImageData.GetScalarComponentAsDouble(x, volume_y - y - 1, 0, 0);
					int intvalue = (int) pixelvalue;
					int highByte = (intvalue & 0xFF00) >> 8;
					int lowByte = intvalue & 0xFF;
					// 低位
					dicomarray[dicomindex++] = (byte) lowByte;
					// 高位
					dicomarray[dicomindex++] = (byte) highByte;
				}
			}
			String outPath = outDir + File.separator + (i + 1) + ".dcm";
			createdicomfile(new File(outPath), dicomarray, fmi, data);
		}
	}

	public static void createdicomfile(File outdcmfile, byte[] buffer, Attributes fmi, Attributes data) {
		try (DicomOutputStream dicomoutstream = new DicomOutputStream(outdcmfile);) {

			dicomoutstream.setEncodingOptions(DicomEncodingOptions.DEFAULT);
			// 使用meta信息,去构建数据集;(meta中含有传输语法)
			dicomoutstream.writeDataset(fmi, data);
			// 写入图像数据部分
			dicomoutstream.writeHeader(Tag.PixelData, VR.OB, buffer.length);
			// offset为输入buffer本身的偏移量
			dicomoutstream.write(buffer, 0, buffer.length);
		} catch (Exception e) {
			e.printStackTrace();
		}
	}
}

vtk切割代码 参考 张晓东 罗火灵 编著的《VTK图形图像开发进阶》