I have several vtkPolyLine objects created from Plot on Intersection Curves filter. How can I use the vtkAppendArcLength filter to add arc_length to its PointData?
For instance, in ParaView, I am trying to get the following script to work (I want to make a nice matplotlib graph):
import numpy as np
import os
from paraview import python_view
from paraview.simple import *
import vtk
import vtk.numpy_interface.dataset_adapter as dsa
from vtk.util import numpy_support
# Disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# Parent directory
root = os.path.dirname(os.path.abspath(__file__))
# Load data
my_data_proxy = LegacyVTKReader(
registrationName="my_data",
FileNames=[
os.path.join(root, "my_data.vtk")
],
)
# Get center and bounds
data = dsa.WrapDataObject(servermanager.Fetch(my_data_proxy))
center = data.GetCenter()
xmin, xmax, ymin, ymax, zmin, zmax = data.GetBounds()
width = xmax - xmin
length = ymax - ymin
height = zmax - zmin
# Create views
data_view = CreateView("RenderView")
python_view = CreateView('PythonView')
# General view properties
view_size = [821, 912]
# Define render view properties
material_library = GetMaterialLibrary()
axes_grid = "GridAxes3DActor"
stereo_type = "Crystal Eyes"
back_end = "OSPRay raycaster"
center_of_rotation = center
view_axes_grid = True
camera_position = [center[0] - width, center[1] - length, center[2]]
camera_focal_point = center
camera_view_up = [0, 0, 1]
camera_focal_disk = 1.0
camera_parallel_scale = 2.5
camera_parallel_projection = True
# Set render view properties
data_view.ViewSize = view_size
data_view.AxesGrid = axes_grid
data_view.CenterOfRotation = center_of_rotation
data_view.StereoType = stereo_type
data_view.BackEnd = back_end
data_view.OSPRayMaterialLibrary = material_library
data_view.AxesGrid.Visibility = 1 if view_axes_grid else 0
data_view.CameraPosition = camera_position
data_view.CameraFocalPoint = camera_focal_point
data_view.CameraViewUp = camera_view_up
data_view.CameraFocalDisk = camera_focal_disk
data_view.CameraParallelScale = camera_parallel_scale
data_view.CameraParallelProjection = 1 if camera_parallel_projection else 0
camera = data_view.GetActiveCamera()
camera.Pitch(-60)
# Set python view properties
python_view.ViewSize = view_size
# Create layout
layout = CreateLayout(name="Analysis")
layout.SplitHorizontal(0, 0.5)
layout.AssignView(1, data_view)
layout.AssignView(2, python_view)
layout.SetSize(1652, 1156)
# Create intersection curves
angle_range = 180
angle_step = 10
for i in range(0, angle_range, angle_step):
curve = PlotOnIntersectionCurves(
registrationName=f"slice_{i:03d}_deg", Input=my_data_proxy
)
curve.SliceType = "Plane"
curve.SliceType.Origin = center
curve.SliceType.Normal = [np.cos(np.radians(90 + i)), np.sin(np.radians(90 + i)), 0.0]
alg = vtk.vtkAppendArcLength()
data = dsa.WrapDataObject(servermanager.Fetch(curve))
for block in range(data.GetNumberOfBlocks()):
obj = data.GetBlock(block)
alg.SetInputData(obj)
alg.Update()
Hide3DWidgets()
SetActiveView(envelope_view)
curve_display = Show(curve, data_view, "GeometryRepresentation")
curve_display.Representation = "Points"
# Show data
SetActiveView(data_view)
Show(my_data_proxy, data_view, "GeometryRepresentation")
data_view.ResetCamera()
python_view.Script = """
import numpy as np
import os
from paraview import python_view
from paraview.simple import *
import vtk
import vtk.numpy_interface.dataset_adapter as dsa
from vtk.util import numpy_support
def setup_data(view):
desired_arrays = ["arc_length", "Data"]
num_visible_objects = view.GetNumberOfVisibleDataObjects()
for ndx in range(num_visible_objects):
obj = view.GetVisibleDataObjectForSetup(ndx)
view.DisableAllAttributeArrays()
for array in desired_arrays:
view.SetAttributeArrayStatus(ndx, vtk.vtkDataObject.POINT, array, 1)
def render(view, width, height):
figure = python_view.matplotlib_figure(width, height)
ax = figure.add_subplot(1,1,1)
ax.minorticks_on()
ax.set_title("My Data")
ax.set_xlabel("Arc Length")
ax.set_ylabel("Data")
num_visible_objects = view.GetNumberOfVisibleDataObjects()
arclengths = np.array([])
data = np.array([])
for ndx in range(num_visible_objects):
num_arrays = view.GetNumberOfAttributeArrays(ndx, vtk.vtkDataObject.POINT)
print(f"Object {ndx} Data Arrays:")
for arr in range(num_arrays):
print(view.GetAttributeArrayName(ndx, vtk.vtkDataObject.POINT, arr))
obj = view.GetVisibleDataObjectForRendering(ndx)
if obj:
if obj.GetDataObjectType() == 13: # MultiBlock
for block in range(obj.GetNumberOfBlocks()):
block_obj = obj.GetBlock(block)
pt_data = block_obj.GetPointData()
s = pt_data.GetArray("arc_length")
d = pt_data.GetArray("Data")
if s:
s_np = vtk_to_numpy(s)
if arclengths.size == 0:
arclengths = s_np
else:
arclengths = np.vstack((arclengths, s_np))
if d:
d_np = vtk_to_numpy(d)
if data.size == 0:
data = d_np
else:
data = np.vstack((data, d_np))
else:
pt_data = obj.GetPointData()
s = pt_data.GetArray("arc_length")
d = pt_data.GetArray("Data")
if s:
s_np = vtk_to_numpy(s)
if arclengths.size == 0:
arclengths = s_np
else:
arclengths = np.vstack((arclengths, s_np))
if d:
d_np = vtk_to_numpy(d)
if data.size == 0:
data = d_np
else:
data = np.vstack((data, d_np))
if arclengths.size > 0 and data.size > 0:
s_avg = np.mean(arclengths, axis=1)
d_avg = np.mean(data, axis=1)
d_std = np.std(data, axis=1)
ax.plot(s_avg, d_avg, label="my_data", color="blue")
ax.fill_between(s_avg, d_avg + 3 * d_std, d_avg - 3 * d_std, facecolor="blue", alpha=0.5)
ax.legend(loc="upper right")
ax.grid()
return python_view.figure_to_image(figure)
"""