- У меня есть сегментированный стек изображений, хранящийся в vtkImageData
- Я создаю сетку изоконтур с помощью vtkSurfaceNets3D (сглаженную, с заполненными отверстиями и т. д.)
- Я делаю это трафарет, чтобы создать новый вокселизированный объем, отслеживая при этом, какие вокселы находятся внутри и снаружи сетки SurfaceNets
- Теперь начинается сложная часть: я пытаюсь определить, какие вокселы на самом деле пересекаются сеткой, и пытаюсь вычислить их относительную долю «снаружи/внутри». Для этого я рассматриваю для каждого вокселя метод Монте-Карло, рассматривая N точек внутри вокселя и вычисляя их расстояние до сетки с помощью vtkImplicitPolyDataDistance.
Я не могу определить, вызвано ли это смещение несоответствием начала координат/расстояния между сеткой и изображением или это связано с особым поведением vtkImplicitPolyDataDistance метод.
Вот мой код ниже для справки. Вычисленное мной значение дроби — frac, как определено функцией voxel_fraction_mc.
import tifffile
import vtk
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from vtk.util.numpy_support import vtk_to_numpy, numpy_to_vtk
plt.close('all')
#%% Parametres
filename = 'sphere_stack.tiff'
resol = 0.02 # taille physique d'un pixel en mm
add_mat = False # rajouter de la matière autour du volume du pore
matId_matiere = 2
# épaisseurs à rajouter en mm dans chaque direction (symmétrique)
ep_z = 2
ep_y = 2
ep_x = 2
convert2VTK = False
#%% Fetch data
# lecture data
img_stack = (tifffile.imread(filename)-1)
if add_mat:
img_stack = np.pad(img_stack,((int(ep_z/resol),int(ep_z/resol)),(int(ep_y/resol),int(ep_y/resol)),(int(ep_x/resol),int(ep_x/resol))),'constant',constant_values=matId_matiere)
img_stack = img_stack.transpose((2,1,0))
# visualitation
id_img = img_stack.shape[2]//2
single_img = img_stack[:,:,id_img].astype(np.int8)
x1 = np.linspace(0.5*resol, img_stack.shape[0]*resol-0.5*resol, img_stack.shape[0])
y1 = np.linspace(0.5*resol, img_stack.shape[1]*resol-0.5*resol, img_stack.shape[1])
z1 = np.linspace(0.5*resol, img_stack.shape[2]*resol-0.5*resol, img_stack.shape[2])
xx1, yy1, zz1 = np.meshgrid(x1,y1,z1, indexing='ij')
centers1 = np.stack([xx1, yy1, zz1], axis=-1)
#%% SurfaceNets3D
from vtk.numpy_interface import dataset_adapter as dsa
def getMeshes(stack,resol,hole_fill_size=10*resol):
"""
Uses VTK SurfaceNets3D to generate pixel coordinate meshes.
stack: 3 dimensional labeled array. Expecting z, y, x dimensions, but
shouldn't matter.
"""
#VTK seems to use x as the first index.
nx,ny,nz = stack.shape
img = vtk.vtkImageData();
img.SetDimensions(nx,ny,nz)
img.SetSpacing(resol,resol,resol)
img.SetOrigin(0,0,0)
flat = stack.ravel(order='F')
vtk_array = numpy_to_vtk(flat, deep=True)
img.GetPointData().SetScalars(vtk_array)
snets = vtk.vtkSurfaceNets3D()
snets.SetInputData(img)
snets.SetOutputMeshTypeToTriangles()
# snets.SmoothingOff()
snets.Update()
# Hole fill
fill = vtk.vtkFillHolesFilter()
fill.SetInputData(snets.GetOutput())
fill.SetHoleSize(hole_fill_size)
fill.Update()
# Smoothing
smooth = vtk.vtkSmoothPolyDataFilter()
smooth.SetInputData(fill.GetOutput())
smooth.SetNumberOfIterations(20)
smooth.Update()
#The normals are not all the same direction.
nrm = vtk.vtkPolyDataNormals()
nrm.ConsistencyOn()
nrm.AutoOrientNormalsOn()
nrm.SetSplitting(False)
nrm.SetInputDataObject(smooth.GetOutputDataObject(0) )
nrm.Update()
polys = nrm.GetOutput()
#polys = snets.GetOutput()
pda = dsa.WrapDataObject(polys)
points = np.array(pda.GetPoints())
return polys, pda, points, nrm
def vtk_to_triangles(pda):
cells = pda.GetPolys()
cell_array = cells.GetData()
polys = vtk_to_numpy(cell_array)
faces = []
i = 0
while i < len(polys):
n = polys
face = polys[i+1:i+1+n]
faces.append(face)
i += n + 1
return np.array(faces)
poly, pda, points, normals = getMeshes(img_stack,resol)
#%% Voxelize
# Définition de la nouvelle voxelisation
voxel_size = 0.8*resol
nx = int(img_stack.shape[0]*resol / voxel_size)
ny = int(img_stack.shape[1]*resol / voxel_size)
nz = int(img_stack.shape[2]*resol / voxel_size)
x = np.linspace(0.5*voxel_size, nx*voxel_size-0.5*voxel_size, nx)
y = np.linspace(0.5*voxel_size, ny*voxel_size-0.5*voxel_size, ny)
z = np.linspace(0.5*voxel_size, nz*voxel_size-0.5*voxel_size, nz)
xx, yy, zz = np.meshgrid(x,y,z, indexing='ij')
centers2 = np.stack([xx, yy, zz], axis=-1)
# Identification des pixels à l'intérieur, à l'extérieur et à l'interface
image = vtk.vtkImageData()
image.SetDimensions(nx, ny, nz)
image.SetOrigin(0,0,0)
image.SetSpacing(voxel_size, voxel_size, voxel_size)
image.AllocateScalars(vtk.VTK_UNSIGNED_CHAR, 1)
image.GetPointData().GetScalars().Fill(1)
pol2stenc = vtk.vtkPolyDataToImageStencil()
pol2stenc.SetInputData(poly)
pol2stenc.SetOutputWholeExtent(image.GetExtent())
pol2stenc.SetOutputSpacing(voxel_size, voxel_size, voxel_size)
pol2stenc.Update()
inside = vtk.vtkImageStencil()
inside.SetInputData(image)
inside.SetStencilConnection(pol2stenc.GetOutputPort())
inside.ReverseStencilOff()
inside.SetBackgroundValue(0) # extérieur = 0
inside.Update()
insideoutput = inside.GetOutput()
inside = insideoutput.GetPointData().GetScalars()
inside = vtk_to_numpy(inside)
inside = inside.reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)
outside = vtk.vtkImageStencil()
outside.SetInputData(image)
outside.SetStencilConnection(pol2stenc.GetOutputPort())
outside.ReverseStencilOn()
outside.SetBackgroundValue(0) # extérieur = 0
outside.Update()
outsideoutput = outside.GetOutput()
outside = outsideoutput.GetPointData().GetScalars()
outside = vtk_to_numpy(outside)
outside = outside.reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)
implicit = vtk.vtkImplicitPolyDataDistance()
implicit.SetInput(poly)
dilate = vtk.vtkImageDilateErode3D()
dilate.SetInputData(insideoutput)
dilate.SetDilateValue(1)
dilate.SetErodeValue(0)
dilate.SetKernelSize(3,3,3)
dilate.Update()
dilate2 = vtk.vtkImageDilateErode3D()
dilate2.SetInputData(outsideoutput)
dilate2.SetDilateValue(1)
dilate2.SetErodeValue(0)
dilate2.SetKernelSize(3,3,3)
dilate2.Update()
erode = vtk.vtkImageDilateErode3D()
erode.SetInputData(insideoutput)
erode.SetDilateValue(0)
erode.SetErodeValue(1)
erode.SetKernelSize(3,3,3)
erode.Update()
boundary = vtk.vtkImageMathematics()
boundary.SetOperationToMultiply()
boundary.SetInput1Data(dilate2.GetOutput())
boundary.SetInput2Data(dilate.GetOutput())
boundary.Update()
intersection = vtk_to_numpy(boundary.GetOutput().GetPointData().GetScalars()).reshape((nz, ny, nx)).transpose(2,1,0).astype(bool)
#%% Labelisation des voxels
# Fraction volumiques des voxels composite
def voxel_fraction_mc(center, voxel_size, implicit, N=50): # fraction volumique monte carlo
pts = np.random.uniform(-0.5,0.5,(N,3))*voxel_size + center
d = [implicit.EvaluateFunction(p) for p in pts]
return sum(di > 0 for di in d)/N # fraction DANS du mesh
# Identification de la nature des voxels composite
frac = np.array([voxel_fraction_mc(pt, voxel_size, implicit) for pt in centers2[intersection].reshape(-1,3)])
centers2_label = np.ones(centers2.shape[:-1])
centers2_label[inside] = 0
centers2_label[outside] = 2
centers2_label[intersection] = 2*frac
#%% Visualisation
# 2D
plane = vtk.vtkPlane()
plane.SetOrigin(0,0,centers2_label.shape[-1]//2*voxel_size)
plane.SetNormal(0, 0, 1)
# --- couper le mesh ---
cutter = vtk.vtkCutter()
cutter.SetInputData(poly) # ton mesh
cutter.SetCutFunction(plane)
cutter.Update()
# --- récupérer les points ---
cut_poly = cutter.GetOutput()
points_vtk = cut_poly.GetPoints()
points = vtk_to_numpy(points_vtk.GetData())
fig,axs = plt.subplots(1,2,figsize=(12,8))
axs[0].pcolormesh(single_img, edgecolors='black', linewidth=0.5, shading='nearest')
axs[0].set_aspect('equal')
axs[1].pcolormesh(centers2_label[:,:,centers2_label.shape[2]//2], edgecolors='black', linewidth=0.5, shading='nearest')
axs[1].set_aspect('equal')
axs[1].scatter(points[:,1]/voxel_size,points[:,0]/voxel_size,color='r',s=1)
axs[0].set_xticks(np.arange(-0.5, img_stack.shape[1], 1), minor=True)
axs[0].set_yticks(np.arange(-0.5, img_stack.shape[0], 1), minor=True)
axs[0].grid(which='minor', color='black', linestyle='-', linewidth=1,alpha=0.5)
axs[0].tick_params(which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
axs[1].scatter(points[:,1]/voxel_size,points[:,0]/voxel_size,color='r',s=1)
axs[1].set_xticks(np.arange(-0.5, centers2_label.shape[1], 1), minor=True)
axs[1].set_yticks(np.arange(-0.5, centers2_label.shape[0], 1), minor=True)
axs[1].grid(which='minor', color='black', linestyle='-', linewidth=1,alpha=0.5)
axs[1].tick_params(which='both', bottom=False, left=False, labelbottom=False, labelleft=False)
