Hello,
I have two sources that are overlapping that have common global ids (they are coming from a common source but divided), to better exemplify, here is a small example to achieve this:
from paraview.simple import *
wavelet1 = Wavelet(registrationName='Wavelet1')
globalPointAndCellIds1 = GlobalPointAndCellIds(registrationName='GlobalPointAndCellIds1', Input=wavelet1)
coordinates1 = Coordinates(registrationName='Coordinates1', Input=globalPointAndCellIds1)
UpdatePipeline()
threshold1 = Threshold(registrationName='Threshold1', Input=coordinates1)
threshold1.Scalars = ['POINTS', 'PointLocations']
threshold1.LowerThreshold = 0.0
threshold1.UpperThreshold = 10.0
threshold2 = Threshold(registrationName='Threshold2', Input=coordinates1)
threshold2.Scalars = ['POINTS', 'PointLocations']
threshold2.LowerThreshold = 0.0
threshold2.UpperThreshold = 10.0
threshold2.Invert = 1
this creates two filters (threshold1 and threshold2) that have inherently the same global ids.
I would like to create a programmable filter that:
- takes the two thresholds as inputs
- gets the shared points, using global Points IDs
- gets the cells that have this points on threshold1 and stores all other points of this cells
- gets another ‘layer’ of cells that have this specific points
I got a prototype working of it:
input1 = self.GetInputDataObject(0, 0) # target mesh (removedLayer)
input2 = self.GetInputDataObject(0, 1) # reference mesh (finalChannel)
output = self.GetOutputDataObject(0)
output.ShallowCopy(input2)
arr1 = input1.GetPointData().GetArray("GlobalPointIds")
gpids1 = {arr1.GetValue(i) for i in range(input1.GetNumberOfPoints())}
nCells = input2.GetNumberOfCells()
import vtk
shared = vtk.vtkIntArray()
shared.SetName("flag")
shared.SetNumberOfTuples(nCells)
IDsToRecheck=set()
cellIDsToNotRecheck=[]
for cid in range(nCells):
cell = input2.GetCell(cid)
ids = cell.GetPointIds()
count = 0
for i in range(ids.GetNumberOfIds()):
pid = ids.GetId(i)
gpid = input2.GetPointData().GetArray("GlobalPointIds").GetValue(pid)
if gpid in gpids1:
count += 1
if count!=0:
shared.SetValue(cid, count)
for i in range(ids.GetNumberOfIds()):
IDsToRecheck.add(ids.GetId(i))
else:
shared.SetValue(cid, 0)
#for i in dir(shared):
# print(i)
#second layer
for cid in range(nCells):
cell = input2.GetCell(cid)
val=shared.GetValue(cid)
if val==0:
ids = cell.GetPointIds()
pid = [ids.GetId(i) for i in range(ids.GetNumberOfIds())]
if any([i in IDsToRecheck for i in pid]):
shared.SetValue(cid, -11)
output.GetCellData().AddArray(shared)
my question is if there is a faster way to do this, without the need to go thought the hole nCells each time. as I am looking to identify the interface and i have a bunch of cells on the bulk that i know that will not be there.
I am thinking of two different approaches (this is where i am looking for help) that not sure how to achieve it:
- instead of going thought the hole range of the nCells in a for loop which can be quite slow in python, i would like to select specific ids inside a bounding box that would be a slight larger (~3 voxels in each direction, of the input1). this is far from optimal as in case that input1 is not a box (as it is in the example) the number of cells in the loop can be unnecessary checked.
- i was thinking on:
a. find the shared GlobalPointIds between the two sources.
b. select using get the internal point ids corresponding to these values of GlobalPointIds (found in a.), ideally not needing this step.
c. select the cells using the information in a.(ideally) or b., something like ‘hey, return me the cell that has this specific value, or that it has this point id inside of itself’
would like some input on how to improve this script. lastly, if someone comes to mind a good learning material for this kind of vtk manipulation in python would appreciate.