I’d like to replicate the Maximum Intensity Projection (‘enhance high values’) feature of Leapfrog, which brings sparse points with high scalar values to the foreground. According to Cowan (2014) this is ‘easily’ achieved by using the z-index ordering capabilities of OpenGL, reordering de z index according to grade value instead of distance from the point of view, even when interacting with the scene. ¿Can z-ordering be accessed through the vtkOpenGLActor class, or is there a better way to achieve the same results?. If this can be done through Pyvista, much better.
Below are some figures that describe what I want to achieve.
Thanks!
From Cowan (2014):
Dense point cloud: not enhanced, low values block higher values at depth:
Dense point cloud: enhanced (MIP), higher values are brought to the foreground:
I’m also very interested in this application! As maximum intensity projection is already available for volume rendering (e.g., 3D Slicer) in vtk, it seems as though the basic machinery is in place and that it wouldn’t be too challenging to port over the same approach for 3D scatters/point clouds or objects rendered from the vtkPolyData class. However, it appears that, even two years after the post above, this has not been done. Might anyone have suggestions on where to start?
I finally made it work! A simple but obscure solution is to access and edit the OpenGL vertex Shader through VTK’s ShaderProperty class. Replacing the z screen coordinate with -scalar value for each vertex did the trick.
If opacity is not 1, depth peeling is required.
Thanks for your interest!
import vtk
import numpy as np
import pyvista as pv
from pyvista import examples
####################################################################
######## MAXIMUM INTENSITY PROJECTION (MIP) FOR POINT CLOUDS #######
######## Ramón Aguirre M. 2024 #######
####################################################################
#Donwload Knee dataset
vol = examples.download_knee_full()
# Get a Sample of vertices from the dataset
sample_n=100000
randix=np.random.randint(0,vol.n_points,sample_n)
pts=vol.extract_points(randix,adjacent_cells=False,include_cells=False).cast_to_pointset()
#Start pyvista plotter
p=pv.Plotter(shape=(1,2))
#Define actor parameters
cmap='jet'
clim=[0,100]
opacity=1
#Plot normally for comparison
p.subplot(0,0)
#add point cloud actor to render window
point_actor=p.add_mesh(pts,cmap=cmap,opacity=opacity,clim=clim)
p.add_title('Normal Projection\n(point cloud)',font_size=12)
#Implement MIP
p.subplot(0,1)
#add point cloud actor to render window
point_actor=p.add_mesh(pts,cmap=cmap,opacity=opacity,clim=clim)
# Maximimum Intensity Projection MIP does not work with opacity <1
# unless Depth Peeling is enabled
if point_actor.prop.GetOpacityMinValue()<1:
p.enable_depth_peeling(number_of_peels=None,occlusion_ratio=0)
# Get the Vertex Shader Program
# it runs once for each vertex every time the render is updated
# MAIN TRICK: we reset the z screen coordinates so vertices
# with higher values are closer to the screen (Cowan, 2014)
property=point_actor.GetShaderProperty()
property.AddShaderReplacement(
vtk.vtkShader.Vertex,
"//VTK::LineWidthGLES30::Impl", # We replace this text in the shader program
True, # before the standard replacements
'''
gl_Position.z = -colorTCoord[0]; //We add this line to set vertices with highest scalar values closer to the screen
//VTK::LineWidthGLES30::Impl", //Keeping the original text
''',
False # only do it once
)
p.add_title('Maximum Intensity Projection\n(point cloud)',font_size=12)
# Final settings and render
p.enable_parallel_projection()
p.link_views()
p.show()
#Reference:
# 2014: Cowan, E.J., ‘X-ray Plunge Projection’— Understanding Structural Geology from Grade Data, AusIMM Monograph 30:
# Mineral Resource and Ore Reserve Estimation — The AusIMM Guide to Good Practice, second edition, 207-220
