C# vtk 3D网格模型嵌入Dicom序列中

将一个3D 网格模型vtkPolyData  嵌入到 Dicom序列 vtkImageData 。

如图将一个球模型与Dicom数据层面进行融合的体绘制。

1.读取Dicom 序列，读取3D模型

vtk 读取Dicom序列不太好用，我用的ITK 读的然后转为 vtkImageData （不知道怎么转的可以看我的另一篇讲解）。

读取vtkPolyData 。

2. 模型 vtkPolyData 转换为 vtkImageData

注意：vtkImaData 的spacing越小 效果越好，但是受到  Dicom的spacing的限制。

3.对模型vtkImageData进行 Z 轴切片 然后在 对应 Dicom 切片 的位置修改数据。

  vtkMatrix4x4 axialElements = new vtkMatrix4x4();
  for (int n = 0; n < 4; n++)
  {
      for (int m = 0; m < 4; m++)
      {
          if (n == m)
          {
              axialElements.SetElement(n, m, 1);
          }
          else
          {
              axialElements.SetElement(n, m, 0);
          }
      }
  }
  vtkPoints points = new vtkPoints();

  vtkMatrix4x4 resliceAxes = new vtkMatrix4x4();
         

  for (int i = 0; i < imgSize[2]; i++)
  { 
      double[] point = new double[3];
      point[0] = itkCenter[0];
      point[1] = itkCenter[1];
      point[2] = itkOrgin[2] + i * itkSpacing[2];

      resliceAxes.DeepCopy(axialElements);
      resliceAxes.SetElement(0, 3, point[0]);
      resliceAxes.SetElement(1, 3, point[1]);
      resliceAxes.SetElement(2, 3, point[2]);
      

      vtkImageReslice itkImageReslice = vtkImageReslice.New();
      itkImageReslice.SetInputData(itkImageData);
      itkImageReslice.SetOutputDimensionality(2);
      itkImageReslice.SetInterpolationModeToLinear();
      itkImageReslice.SetResliceAxes(resliceAxes);
      itkImageReslice.Update();

      vtkImageData itkResliceData = new vtkImageData();
      itkResliceData.DeepCopy(itkImageReslice.GetOutput());
    
      
      vtkImageReslice stlImageReslice = vtkImageReslice.New();
      stlImageReslice.SetInputData(stlImageData);
      stlImageReslice.SetOutputDimensionality(2);
      stlImageReslice.SetInterpolationModeToLinear();
      stlImageReslice.SetResliceAxes(resliceAxes);
      stlImageReslice.Update();

      vtkImageData stlResliceData = new vtkImageData();
      stlResliceData.DeepCopy(stlImageReslice.GetOutput());

      ImageDataBlend(ref itkImageData, itkResliceData,  stlResliceData, i);

      Console.WriteLine(" i = " + i.ToString());
  }

 // 切片融合
 public void ImageDataBlend(ref vtkImageData itkImage, vtkImageData originImage, vtkImageData blendImage, int nz)
 {
     try
     {
         // 确定范围
         double[] originBound = originImage.GetBounds();
         double oriMinWidth = originBound[0];
         double oriMaxWidth = originBound[1];
         double oriMinHeight = originBound[2];
         double oriMaxHeight = originBound[3];
         int[] oriDim = originImage.GetDimensions();
         double[] oriSpacing = originImage.GetSpacing();
         double[] oriOrigin = originImage.GetOrigin();
          
         vtkDataArray oriDataArray = originImage.GetPointData().GetScalars();
     
         double[] blendBound = blendImage.GetBounds();
         double bleMinWidth = blendBound[0];
         double bleMaxWidth = blendBound[1];
         double bleMinHeight = blendBound[2];
         double bleMaxHeight = blendBound[3];
         int[] bleDim = blendImage.GetDimensions();
         double[] bleSpacing = blendImage.GetSpacing();
         double[] bleorigin = blendImage.GetOrigin();

         vtkDataArray bleDataArray = blendImage.GetPointData().GetScalars();
         int ny = 0;
         for (double y = oriOrigin[1]; y < oriMaxHeight; y = y + oriSpacing[1])
         {
             ny++;
             if (y < bleMinHeight || y > bleMaxHeight)
             {
                 continue;
             }
             int nx = 0;
             for (double x = oriOrigin[0]; x < oriMaxWidth; x = x + oriSpacing[0])
             {
                 nx++;
                 if (x < bleMinWidth || x > bleMaxWidth)
                 {
                     continue;
                 }
                 // 通过坐标获取 blendImage 像素 index
                 int dy = (int)Math.Ceiling((y - bleMinHeight) / bleSpacing[1]);
                 int dx = (int)Math.Ceiling((x - bleMinWidth) / bleSpacing[0]);

                 int bleIndex = (int)Math.Ceiling((y - bleMinHeight) / bleSpacing[1]) * bleDim[0]
                    + (int)Math.Ceiling((x - bleMinWidth) / bleSpacing[0]);


                 int oriIndex =  ny * oriDim[1] + nx;
                 if (bleIndex > bleDataArray.GetNumberOfValues() - 1 || bleIndex <0)
                 {
                     continue;
                 }
                 
                 int bleValue = (int)blendImage.GetScalarComponentAsDouble(dx, dy, 0, 0);

                 if (bleValue == 0)
                 {
                     continue;
                 }
                 if (oriIndex > oriDataArray.GetNumberOfValues() - 1)
                 {
                     continue;
                 }
                 itkImage.SetScalarComponentFromDouble(nx, ny, nz, 0, 51.0);
             }
         }
     }
     catch (Exception ex)
     {
         Console.WriteLine(  "err");
         
     }
    
 }