pyvista / pyvista

3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK)
https://docs.pyvista.org
MIT License
2.73k stars 501 forks source link

Points Marker Style #165

Closed banesullivan closed 2 years ago

banesullivan commented 5 years ago

matplotlib figure of some hard data points and an interpolated grid

download

vtki figure of the same data (note how truth points are barely visible)

download

banesullivan commented 2 years ago

AFAIK, this isn't possible

ramaguirre commented 1 month ago

Hello @banesullivan

yes to all three, this can be done by updating the Fragment Shader. I was referred to this post (https://discourse.vtk.org/t/different-shaped-vertex-points/10136).

Example below, combining Maximum Intensity Projection with circular points, with and without boders. It works by discarding pixels located outside a radius of 0.5, and optionally assigning ambientColor = diffuseColor = (0,0,0) for a radius of 0.3 for a black outline.

There is a lot of potential in Shader replacement.

image

import vtk
import numpy as np
import pyvista as pv
from pyvista import examples
import matplotlib.pyplot as plt

####################################################################
######## 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,3))
p.enable_parallel_projection()

#Define actors parameters
cmap='jet'
clim=[0,100]
opacity=1
rpas=False
lighting=False

#Plot normally for comparison
p.subplot(0,0)
#add point cloud actor to render window
point_actor=p.add_mesh(pts.flip_z(),cmap=cmap,opacity=opacity,clim=clim,
                       render_points_as_spheres=rpas,lighting=lighting,point_size=10)

p.add_title('Knee dataset\nNormal Projection\n(point cloud)',font_size=8)

#Implement MIP
p.subplot(0,1)
#add point cloud actor to render window
point_actor=p.add_mesh(pts.flip_z(),cmap=cmap,opacity=opacity,clim=clim,
                       render_points_as_spheres=rpas,lighting=lighting,point_size=10)

# 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()

max_scalar=np.nanmax(point_actor.mapper.dataset.active_scalars)
min_scalar=np.nanmin(point_actor.mapper.dataset.active_scalars)

property.AddShaderReplacement(
    vtk.vtkShader.Vertex,
    "//VTK::LineWidthGLES30::Impl", # We replace this text in the shader program
    True, # before the standard replacements
    f'''
    gl_Position.z = -(abs(colorTCoord[0])-{min_scalar})/({max_scalar}-{min_scalar});  //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
)

property.AddShaderReplacement(
    vtk.vtkShader.Fragment,  # in the fragment shader
     "//VTK::Color::Impl", # replace the color block
     True, # before the standard replacements
     '''
     //VTK::Color::Impl // we still want the default calc
     //error
     float dist = length(vec2(gl_PointCoord.x - 0.5, gl_PointCoord.y - 0.5)); //Distance to centroid
     if (dist > 0.3) { ambientColor = vec3(0,0,0); //Set border colors to black
     diffuseColor = vec3(0,0,0);}   //Set border colors to black
     if (dist > 0.5) {discard;}''',  #but we add this: remove points outside of radius
     False # only do it once
)

p.add_title('Knee dataset\nMaximum Intensity Projection\n(point cloud)',font_size=8)

p.subplot(0,2)
#add point cloud actor to render window
point_actor=p.add_mesh(pts.flip_z(),cmap=cmap,opacity=opacity,clim=clim,
                       render_points_as_spheres=rpas,lighting=lighting,point_size=10)

# 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()

max_scalar=np.nanmax(point_actor.mapper.dataset.active_scalars)
min_scalar=np.nanmin(point_actor.mapper.dataset.active_scalars)

#Implement MIP
property.AddShaderReplacement(
    vtk.vtkShader.Vertex,
    "//VTK::LineWidthGLES30::Impl", # We replace this text in the shader program
    True, # before the standard replacements
    f'''
    gl_Position.z = -(abs(colorTCoord[0])-{min_scalar})/({max_scalar}-{min_scalar});  //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
)

#Make markers circular
property.AddShaderReplacement(
    vtk.vtkShader.Fragment,  # in the fragment shader
     "//VTK::Color::Impl", # replace the color block
     True, # before the standard replacements
     '''
     //VTK::Color::Impl // we still want the default calc
     //error
     float dist = length(vec2(gl_PointCoord.x - 0.5, gl_PointCoord.y - 0.5)); //Distance to centroid
     if (dist > 0.5) {discard;}''', #but we add this: remove points outside of radius
     False # only do it once
)

p.add_title('Knee dataset\nMaximum Intensity Projection\n(point cloud)',font_size=8)

# Final settings and render

p.link_views()

#p.show(jupyter_backend) # MIP doas not work with 'html' backend
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