Thanks for your response @lassoan.
These are the steps that I want to follow:
Initially, I have a surface model that is closed. I cut a surface from that model and this new open surface model is what I want to study.
I want to select a specific region (a circle with a specific radius) on this open surface model and assign a different label to this region.
I appreciate your help.
Do you want to label the points or cells? What is the preferred behavior at the boundary (subdivide boundary cells to make the intersection boundary smooth; or keep the original cells and have crinkled boundary)? Do you want to measure the distance on the surface or simple Euclidean distance in 3D space?
What are you going to do with the result? Just color it differently or remove cells that are farther away? If you remove cells then do you need to keep the original point and cell IDs?
Thanks for your response @lassoan.
I want to label cells. In the intersection region, I want to keep the original cell and have a crinkled boundary.
I think I did not understand what you mean by measuring the distances. However, I want to fit a closed circular curve with a specific radius to this surface and then assign a different label to cells inside it.
The purpose: I want to use these areas with different labels in my FEM analysis in the next steps. Also, I want to keep both regions outside, and inside of the circular curve. However, I do not need to keep the original cell IDs.
You can create a “region” cell array like this:
centerFid = getNode('F')
modelNode = getNode('MyModel')
distance = 15.0
# Create an array that contains the distances from the chosen center point
import numpy as np
centerPos = centerFid.GetNthControlPointPositionVector(0)
coords = slicer.util.arrayFromModelPoints(modelNode)
vectorsToCenter = coords - centerPos
distances = np.linalg.norm(vectorsToCenter, axis=1)
# Create a new "region" point data array that contains the region ID
regionArrayVtk = vtk.vtkDoubleArray()
regionArrayVtk.SetName('region')
regionArrayVtk.SetNumberOfValues(len(distances))
modelNode.GetPolyData().GetPointData().AddArray(regionArrayVtk)
regionArray = slicer.util.arrayFromModelPointData(modelNode, 'region')
regionArray[distances<distance] = 0.0
regionArray[distances>=distance] = 1.0
slicer.util.arrayFromModelPointDataModified(modelNode, 'region') # indicate that we finished with the updates
# Visualize "region" point data
modelNode.GetDisplayNode().SetActiveScalar("region", vtk.vtkAssignAttribute.POINT_DATA)
modelNode.GetDisplayNode().SetAndObserveColorNodeID('vtkMRMLColorTableNodeFileViridis.txt')
modelNode.GetDisplayNode().SetScalarVisibility(True)
# Create cell data from "region" array
toCellData = vtk.vtkPointDataToCellData()
toCellData.SetInputData(modelNode.GetPolyData())
toCellData.PassPointDataOn()
toCellData.Update()
modelNode.SetAndObservePolyData(toCellData.GetOutput())
1 Like
@lassoan looks like theres a missing piece of the script just before the third from the last line.
Thank you, fixed. It was just an incompletely removed commented out line.
1 Like
Thanks for the code @lassoan. It works great. However, I have two problems:
-
When I run the code, the region is selected and appears in a different color. However, how can I export this region to a new model?
-
The code runs without error when the surface is closed. However, when I try it with an open surface (that is my interest) I get an error when I write the following part of the code:
coords = slicer.util.arrayFromModelPoints(modelNode)
and the error is this:
Traceback (most recent call last):
File "<console>", line 1, in <module>
File "/home/user/Slicer-4.11.20210226-linux-amd64/bin/Python/slicer/util.py", line 1477, in arrayFromModelPoints
pointData = modelNode.GetPolyData().GetPoints().GetData()
AttributeError: 'NoneType' object has no attribute 'GetPoints'
I greatly appreciate it if you can help me solve these two problems.
It is saved in the model when you use a suitable format, such as VTK, VTP/VTU. For example, FEBio studio can read such files directly. Other software may not be able to read VTK file formats, so you need to convert to a suitable one, for example by using meshio Python package.
The code runs without error when the surface is closed. However, when I try it with an open surface (that is my interest) I get an error when I write the following part of the code
Good catch. We should replace GetPolyData by GetMesh in that file. You can make this change locally on your computer (just change the .py file, save it, and restart Slicer).
Thanks for your response @lassoan. I have to use meshio to convert the vtk file to the format that I want. The problem is this: now I can only export my model to vtk (polydata). However, it seems that meshio only reads vtk unstructured grid. How can I convert my polydata vtk to unstructured data vtk?
Also, I replaced the code as you mentioned in my computer, however, I still get the same error:
>>> coords = slicer.util.arrayFromModelPoints(modelNode)
Traceback (most recent call last):
File "<console>", line 1, in <module>
File "/home/user/Slicer-4.11.20210226-linux-amd64/bin/Python/slicer/util.py", line 1477, in arrayFromModelPoints
pointData = modelNode.GetMesh().GetPoints().GetData()
AttributeError: 'NoneType' object has no attribute 'GetPoints'
Your help is greatly appreciated.
You can copy the points and cells of a polydata into an unstructured grid but maybe the meshio has problems with polygon cells. You can try to save as OBJ and PLY, too, and see what point/cell data arrays are preserved.
If GetMesh() returns None then most likely the model node is empty (are you sure you get the correct node from the scene?).
Hello @lassoan. I will work more on the vtk format later. However, now I think I should get the code to work first. I tried several times, however, I still get that error related to the GetPoints for the opne surface model. This is what I do:
I create a point using the Markups module (the point F in the code you sent me) and then I use your code. I do not think the problem is related to this point. Because, in my closed surface model, the code runs even when the point F is outside of the model.
Can you please tell me how can I make sure whether I get the correct node from the scene?
Can you save your scene in mrb file format, upload it somewhere and post the link here? Also copy here the Python script that you are executing.
Thanks for your response @lassoan.I tried several times again and eventually, I found that the problem is solved. I could not understand what was my mistake in the previous runs, however, it is working now.
My only problem right now is the detection of the label. I found that vtk format can preserve this information. I say this because when I open the vtk file in a new slicer window and go to the Models module, I see that number of point scalars is equal to 2 and I assume that this means the region is saved in the file.
Are these regions same as labels now? How can I assign each region to a new model/segmentation?
I have uploaded one of my vtk models here.
Thanks in advance.
You can use vtkThreshold filter to extract a part of the mesh based on point or cell scalar value.
Thanks for your response @lassoan. I am searching for it now, I just wanted to know whether you have any specific example in mind does the same thing or a similar task?
There are tons of VTK Examples on the web. You can get the vtkPolyData object from the model node by calling GetMesh() method, then run the VTK filters, and then call slicer.modules.models.logic().AddModel(...) method to add a new model from the filter output.
Thanks for your response @lassoan. I got this code based on my search:
SMesh=getNode('Model').GetMesh()
tfil=vtk.vtkThreshold()
tfil.SetInputData(SMesh)
tfil.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, vtk.vtkDataSetAttributes.SCALARS );
tfil.ThresholdByLower( 1 )
tfil.Update()
slicer.modules.models.logic().AddModel(tfil.GetOutput())
However, I get an error in the last step. Can you please let me know what is the problme with my code?
This is my sample as well.
Thanks in advance,
Without seeing the error message I have no chance of guessing what could go wrong.
This is the error:
Traceback (most recent call last):
File "<console>", line 1, in <module>
TypeError: arguments do not match any overloaded methods
Thanks
I managed to get two different models using the following code:
meshNode = getNode('Model')
mesh = meshNode.GetMesh()
cellData = mesh.GetCellData()
labelsRange = cellData.GetArray("region").GetRange()
for labelValue in range(int(labelsRange[0]), int(labelsRange[1]+1)):
threshold = vtk.vtkThreshold()
threshold.SetInputData(mesh)
threshold.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, "region")
threshold.ThresholdBetween(labelValue, labelValue)
threshold.Update()
if threshold.GetOutput().GetNumberOfPoints() > 0:
modelNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLModelNode", "{0}_{1}".format(meshNode.GetName(), labelValue))
modelNode.SetAndObserveMesh(threshold.GetOutput())
modelNode.CreateDefaultDisplayNodes()
However, as shown in the picture below, a part of the main model near the region that I had chosen does not exist in the created models. What is the reason for this?
Also, I do not know why the color of my model changes to this rainbow-like colorand I cannot change it in the slicer.
My main model is This.