You are right.
Running the script gives me a core dump.
I will investigate more.
If there is a break in vtkPolyDataEdgeConnectivityFilter please let me know, and if possible provide sample data / script.
See above.
Please download the mesh3D.vtk file with a:
wget -q -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/pyvista/mesh3D.vtk
The script runs with vtkPolyDataConnectivityFilter(), not with vtkPolyDataEdgeConnectivityFilter()
Hello,
So, a core dump means you had an abend. Besides the core dump, was there any higher-level messages (e.g. Python stack trace or OpenGL error messages)? If not, even a core dump can provide some insight to what went wrong depending on where exactly the crash ensued.
Can you identify which executable crashed? Depending on the answer, it’ll be easy or difficult to analyze the dump. For example, if the crash occured inside python or libvtkCommonCore-9.1.dll you can get debug versions of these with relative ease. If the crash took place inside one of the OpenGL backend libraries, you can have a hard time getting the debug version of them because they are supplied by your graphics card vendor.
Once you have the debug version of the binary that crashed, you can analyze the core dump with gdb: linux - How do I analyze a program's core dump file with GDB when it has command-line parameters? - Stack Overflow
regards,
Paulo
Hi,
Back on this issue that has stayed unresolved.
Any tries or hints to share ?
Hello,
Did you try to investigate the core dump?
best,
PC
No. I haven’t been able to investigate the core dump.
Honnestly I was expecting that someone confirms that problem.
The solution to change
connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
by
connectivityFilter = vtk.vtkEdgesPolyDataConnectivityFilter()
is easy to test.
For me it causes a core dumped.
An abend like that does not necessarily happen with others. So, only posting the high-level code that triggers it does not mean that much. Hence, it can be a long time until someone experiences the same problem as yours as these mysterious crashes often have a multitude of causes. Python code depends on piles and piles of lower level APIs. The crash can be starting in any of them. If you wish the community to help you, the least you can do is to get human-readable information on the dump following the steps presented and share them here.
If you don’t feel like getting into the trouble of using a debugger to have an inteligible core dump, I believe your options are:
a) waiting;
b) use a work around;
c) hire someone more experienced to do it for you.
Anyway, have you tried a trivial data set? Maybe the crash is caused by some oddity in the data.
take care,
Paulo
I am investigating trying to understand why changing the filter causes an abend.
Have you tried the code I have proposed ? It is quite simple.
I will try to reduce it more and test indeed this vtkPolyDataEdgeConnectivityFilter that causes problem.
In fact the problem I am facing is how to describe anti-clockwise a polygon. The vtk.vtkPolyDataConnectivityFilter() extract correctly a region. The vtk.vtkStripper() and stripper.SetJoinContiguousSegments(True) joins each segments as one entity but I need to describe the polygon with an anti-clockwise convention.
How to force this ?
I have posted this question related to my problem.
Use the polygon area method, if the returned area is negative, then it is counter-clockwise (for Cartesian system such as VTK) or clockwise (for inverted Y axis such as screen coordinates): https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order
These methods are reliable in case of self intersecting polygons provided there is a great assimetry in the proportion of clockwise and anti-clockwise portions such as your example just above. One can see that it is “mostly” clockwise or counter-clockwise.
If the computed winding order is not that you want, just invert the list of vertexes.
One can see that it is “mostly” clockwise or counter-clockwise.
This is exactly the problem. See orient function should work on folded polygons · Issue #1705 · shapely/shapely · GitHub
The solution could be to separate the polydata into convex polygons that then can be oriented correctly.
What do you think of this ? And how to do it ?
The vtkPolyDataEdgeConnectivityFilter works fine, it is crashing because it needs polygons, not just lines.
Here is a solution in vedo, (uncomment reverse() to have it clockwise):
from vedo import *
m = Mesh("poly1.vtk").clean().triangulate()
ms = m.split(must_share_edge=True) # return a list of meshes
s1 = ms[0].boundaries().join_segments()[0]#.reverse()
print("First polyline", s1.points())
show(ms, s1, s1.labels2d(), axes=1)
First polyline
[[ 1562187.9 -5830366. 2038637.4]
[ 1764714.6 -5772294. 2038637.4]
[ 1965086. -5707191.5 2038637.4]
[ 2163063.2 -5635136. 2038637.4]
[ 2358410.5 -5556212.5 2038637.4]
[ 2331017.8 -5491532.5 2235859.8]
[ 2137944.8 -5569539. 2235859.8]
[ 1942269.2 -5640759. 2235859.8]
[ 1744224.8 -5705107.5 2235859.8]
[ 1544058. -5762504.5 2235859.8]]
pip install -U git+https://github.com/marcomusy/vedo.git
This is what I meant by some oddity in the dataset causing the crash. Hence, I believe @PBrockmann could join the lines with coincident vertexes to form a single polyline and then process it with vtkPolyDataEdgeConnectivityFilter : vtkStripper not joining contiguous lines to form a single vtkPolyLine.
May I ask you how you “join the lines with coincident vertexes to form a single polyline” ?
With vtk python API.
@marcomusy : thank you very much for this real test. I didn’t know vedo. I would try to stay witk vtk/pyvista.
Well, port it. It’s not difficult. The name of the classes and methods are much like the same.
For example:
vtkSmartPointer<vtkPolyData> poly = vtkSmartPointer<vtkPolyData>::New();
poly->SetPoints( points );
poly->SetLines( segments );
poly->GetCellData()->SetScalars( dataIndexes );
Ports to:
poly = vtk.vtkPolyData()
poly.SetPoints( points )
poly.SetLines( segments )
poly.GetCellData().SetScalars( dataIndexes )
I can reproduce what the vedo code from @marcomusy has proposed with the VTK API.
But the split function from the vedo code gives anyway an uncorrect number of regions.
I have prepared 3 vtk files to play with.
https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_11.vtk
https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_07.vtk
https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/poly_03.vtk
In the poly_11.vtk, I expect to find 2 regions.
In the poly_07.vtk, 3 regions.
In the poly_03.vtk, 1 region.
from vedo import *
m = Mesh("poly_07.vtk").clean().triangulate()
ms = m.split(must_share_edge=True) # return a list of meshes
print(len(ms))
Gives 1 instead of 3.
Here is the solution I can proposed.
Described in this answer as well.
Notice that all the different self intersected polygons are separated from their duplicate points into multiple polygons.
import pyvista as pv
import vtk
import random
import numpy as np
! wget -nc https://thredds-su.ipsl.fr/thredds/fileServer/ipsl_thredds/brocksce/tmp/mesh3D_03.vtk
mesh = pv.PolyData('mesh3D_03.vtk')
edges = mesh.extract_feature_edges(boundary_edges=True)
pl = pv.Plotter()
R = 6371E3
pl.add_mesh(pv.Sphere(radius=R*0.999, theta_resolution=360, phi_resolution=180))
pl.add_mesh(mesh, show_edges=True, edge_color="gray")
regions = edges.connectivity()
regCount = len(set(regions.get_array("RegionId")))
connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
stripper = vtk.vtkStripper()
allPolys = []
for r in range(regCount):
connectivityFilter.SetInputData(edges)
connectivityFilter.SetExtractionModeToSpecifiedRegions()
connectivityFilter.InitializeSpecifiedRegionList()
connectivityFilter.AddSpecifiedRegion(r)
connectivityFilter.Update()
stripper.SetInputData(connectivityFilter.GetOutput())
stripper.SetJoinContiguousSegments(True)
stripper.Update()
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(stripper.GetOutput())
cleaner.Update()
reg = cleaner.GetOutput()
print(r, reg.GetNumberOfCells(), reg.GetNumberOfPoints())
#==============================
polys = []
for cellIndex in range(reg.GetNumberOfCells()):
c = reg.GetCell(cellIndex)
print(" Cell #%d Number of points: %d" %(cellIndex, c.GetNumberOfPoints()))
d = c.GetPointIds()
ids = []
for idIndex in range(d.GetNumberOfIds()):
ids.append(d.GetId(idIndex))
unique, count = np.unique(ids, return_counts=True)
dup = unique[count > 1]
for id in dup[::-1]:
select = np.where(ids == id)[0]
polys.append(ids[select[0]:select[1]+1])
idsIndices = list(range(select[0],select[1]))
ids = list(np.delete(ids, idsIndices))
print(" ", len(polys))
for p in polys:
points = []
for id in p:
points.append(reg.GetPoint(id))
m1 = pv.MultipleLines(points)
print(" add poly")
allPolys.append(points)
random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
pl.add_mesh(m1, color=random_color, line_width=4)
viewer = pl.show(jupyter_backend='trame', return_viewer=True)
display(viewer)
Then my intention was to produce maps with filled polygons in various projections.
have a look to
The polygon file polys_03.pickle can be produced from
import shapely
from pyproj import Transformer
#=======================================================
# From a cartesian 3D CS (geocentric) to a global latitude/longitude CS
trans_GPS_to_XYZ = Transformer.from_crs(4978, 4326, always_xy=True)
#=======================================================
# Define a function to make continuous the polygon
# because the projection transformation produces output
# between -180 and 180
def makePolygonContinuous(arrayLonLat):
for i in range(1,len(arrayLonLat)):
lon1 = arrayLonLat[i-1][0]
lon2 = arrayLonLat[i][0]
if (abs(lon2-360-lon1) < abs(lon2-lon1)):
arrayLonLat[i][0] = lon2-360
elif (abs(lon2+360-lon1) < abs(lon2-lon1)):
arrayLonLat[i][0] = lon2+360
return arrayLonLat
#=======================================================
allPolys2 = []
for polyIndex in range(len(allPolys)):
print(polyIndex, len(allPolys[polyIndex]))
points = []
for point in allPolys[polyIndex]:
x,y,z = point
lon, lat, Z = trans_GPS_to_XYZ.transform(x,y,z)
points.append([lon, lat])
points = makePolygonContinuous(points)
polygonShapely = shapely.Polygon(points)
polygonShapelyOriented = shapely.geometry.polygon.orient(polygonShapely)
allPolys2.append(polygonShapelyOriented)
#=======================================================
import pickle
pickle.dump(allPolys2, open("polys_03.pickle","wb"))
Thank you for the different help messages received.
2 Likes