"""Color palettes, opacity curves and HU lookup tables for rendering medical volumes."""
import numpy as np
import matplotlib.pyplot as plt
import pyvcad as pv
__all__ = [
"raw_color_palettes",
"get_color_palette",
"lerp",
"interp_color",
"linear_opacity_base",
"sigmoid_opacity_base",
"re_range_opacity",
"create_linear_gradient_hu_map",
"plot_opacity_function",
"plot_hu_color_map",
]
#: Named color palettes. Each value is a list of ``(r, g, b)`` tuples with channels in 0-255,
#: ordered from the low end of the scale to the high end.
raw_color_palettes = {
"grayscale":[(80,80,80),(255,255,255)],
"jet": [(37,94,199),(5,255,255),(255,255,0)],
"rainbow":[(55,39,130),(72,107,222),(78,184,231),(74,221,139),(166,244,67),(241,226,53),(246,144,18),(218,61,0),(126,0,0)],
"viridis":[(127,14,156),(77,108,186),(0,255,242),(94,201,98)],
"cividis":[(0,34,77),(67,85,130),(128,129,145),(190,172,129)],
"magma":[(0,0,4),(84,17,68),(187,55,84),(249,142,8)],
}
[docs]
def get_color_palette(name):
"""
Returns a named palette from :data:`raw_color_palettes`.
Arguments:
name (str): Palette name: "grayscale", "jet", "rainbow", "viridis", "cividis" or "magma".
Returns:
list[tuple[int, int, int]]: The palette's ``(r, g, b)`` stops, with channels in 0-255.
Raises:
ValueError: If no palette has this name.
Example:
>>> import pyvcad_medical as med
>>> palette = med.get_color_palette("viridis")
"""
if name in raw_color_palettes:
return raw_color_palettes[name]
else:
raise ValueError(f"Color palette '{name}' not found.")
[docs]
def lerp(a,b,t):
"""
Linearly interpolates between two numbers.
Arguments:
a (float): Value at ``t = 0``.
b (float): Value at ``t = 1``.
t (float): Interpolation parameter. Values outside [0, 1] extrapolate.
Returns:
float: ``a + (b - a) * t``.
"""
return a+(b-a)*t
[docs]
def interp_color(palette,t):
"""
Samples a palette at a normalized position.
The palette stops are spaced evenly over [0, 1] and adjacent stops are linearly interpolated.
Arguments:
palette (list[tuple[int, int, int]]): ``(r, g, b)`` stops with channels in 0-255, such as a
value from :data:`raw_color_palettes`.
t (float): Position in [0, 1]. Values of 1 or more return the last stop.
Returns:
tuple[float, float, float]: The ``(r, g, b)`` color with channels in [0, 1].
Example:
>>> import pyvcad_medical as med
>>> med.interp_color(med.get_color_palette("grayscale"), 0.5)
"""
n=len(palette)-1; x=t*n; i=int(x); f=x-i
if i>=n: return tuple(v/255 for v in palette[-1])
c1,c2=palette[i],palette[i+1]
return tuple((lerp(c1[j],c2[j],f))/255 for j in range(3))
[docs]
def linear_opacity_base(t: float, ymin: float = 0.0, ymax: float = 1.0) -> float:
"""
Linear opacity curve.
Arguments:
t (float): Normalized position. It is clamped to [0, 1].
ymin (float): Opacity at ``t = 0``.
ymax (float): Opacity at ``t = 1``.
Returns:
float: ``ymin + (ymax - ymin) * t``.
"""
t_clamped = max(0.0, min(1.0, t))
return ymin + (ymax - ymin) * t_clamped
[docs]
def sigmoid_opacity_base(t: float, k: float = 10.0, x0: float = 0.5, ymin: float = 0.0, ymax: float = 1.0) -> float:
"""
Sigmoid (S-curve) opacity curve.
Arguments:
t (float): Normalized position in [0, 1].
k (float): Steepness of the transition. Larger values give a sharper step.
x0 (float): Midpoint of the transition in [0, 1].
ymin (float): Lower bound of the output.
ymax (float): Upper bound of the output.
Returns:
float: ``ymin + (ymax - ymin) / (1 + exp(-k * (t - x0)))``.
"""
import math
s = 1.0 / (1.0 + math.exp(-k * (t - x0)))
return ymin + (ymax - ymin) * s
[docs]
def re_range_opacity(t: float, function, lower_opacity, upper_opacity) -> float:
"""
Rescales an opacity curve to a new output range.
Arguments:
t (float): Normalized position passed to ``function``.
function (callable): Base opacity curve ``f(t)``. Its result is clamped to [0, 1].
lower_opacity (float): Output when ``function(t)`` is 0.
upper_opacity (float): Output when ``function(t)`` is 1.
Returns:
float: ``lower_opacity + (upper_opacity - lower_opacity) * clamp(function(t))``.
Example:
>>> import pyvcad_medical as med
>>> opacity = lambda t: med.re_range_opacity(t, med.sigmoid_opacity_base, 0.1, 0.9)
"""
t_clamped = max(0.0, min(1.0, function(t)))
return lower_opacity + (upper_opacity - lower_opacity) * t_clamped
[docs]
def create_linear_gradient_hu_map(lower_hu, upper_hu, palette=raw_color_palettes["jet"], steps=20, opacity_function=lambda t: 1.0):
"""
Builds a Hounsfield-unit lookup table that sweeps through a color palette.
The range [lower_hu, upper_hu] is split into ``steps`` contiguous bins. Bin ``i`` takes the palette
color at ``t = i / (steps - 1)`` and the opacity ``opacity_function(t)``.
Arguments:
lower_hu (float): Lower HU bound of the first bin.
upper_hu (float): Upper HU bound of the last bin.
palette (list[tuple[int, int, int]]): ``(r, g, b)`` stops with channels in 0-255. Defaults to
the "jet" palette.
steps (int): Number of bins.
opacity_function (callable): Opacity curve ``f(t)`` for ``t`` in [0, 1]. Defaults to a
constant 1.0.
Returns:
list[LookupTableEntry]: One entry per bin, mapping an HU range to an RGBA ``Vec4``.
Example:
>>> import pyvcad as pv
>>> import pyvcad_medical as med
>>> entries = med.create_linear_gradient_hu_map(100, 1000, palette=med.get_color_palette("rainbow"), steps=30)
>>> converter = pv.LookupTableConverter(["hu"], ["color_rgba"], entries, pv.InterpolationMode.STEP)
"""
entries = []
hu_range = upper_hu - lower_hu
step_size = hu_range / steps
for i in range(steps):
# Calculate geometry: contiguous blocks
hu_range_lower = lower_hu + i * step_size
hu_range_upper = lower_hu + (i + 1) * step_size
# Calculate color/opacity parameter: distributed 0..1
t = i / (steps - 1) if steps > 1 else 0.0
color = interp_color(palette, t)
# Map t -> [0,1] via the provided function
opacity = opacity_function(t)
entries.append(pv.LookupTableEntry(hu_range_lower, hu_range_upper, pv.Vec4(color[0], color[1], color[2], opacity)))
return entries
[docs]
def plot_opacity_function(opacity_function, steps=100):
"""
Plots an opacity curve over [0, 1] with matplotlib and shows the figure.
Arguments:
opacity_function (callable): Opacity curve ``f(t)``.
steps (int): Number of evenly spaced samples. Must be at least 2.
"""
ts = [i / (steps - 1) for i in range(steps)]
opacities = [opacity_function(t) for t in ts]
plt.plot(ts, opacities)
plt.xlabel('Normalized t (0 to 1)')
plt.ylabel('Opacity')
plt.title('Opacity Function')
plt.grid(True)
plt.show()
[docs]
def plot_hu_color_map(hu_color_map, lower_hu=None, upper_hu=None, use_transparency=True):
"""
Plots an HU lookup table as a horizontal color bar and shows the figure.
Each entry's HU range is drawn at its true position and width; gaps between ranges are drawn
transparent.
Arguments:
hu_color_map (list): Entries ordered by HU, either ``LookupTableEntry`` objects or legacy
``([lower_hu, upper_hu], color)`` pairs.
lower_hu (float, optional): Left end of the bar. Defaults to the first entry's lower bound.
upper_hu (float, optional): Right end of the bar. Defaults to the last entry's upper bound.
use_transparency (bool): Draw each color with its alpha over a checkerboard. If False, colors
are drawn opaque.
"""
import matplotlib.colors as mcolors
# Determine data bounds
first = hu_color_map[0]
last = hu_color_map[-1]
if hasattr(first, 'input_low'):
map_min = first.input_low
map_max = last.input_high
else:
map_min = first[0][0]
map_max = last[0][1]
# Use provided bounds or fallback to data extent
vmin = lower_hu if lower_hu is not None else map_min
vmax = upper_hu if upper_hu is not None else map_max
# Safety check for empty range
if vmax <= vmin:
vmax = vmin + 1.0
# Construct the colormap segments (list of (position 0..1, color))
# This approach allows for arbitrary gaps and widths.
segments = []
trans_white = (1.0, 1.0, 1.0, 0.0) # Fully transparent white for gaps
def get_norm(val):
# Normalize value to 0..1 relative to the plot range
return max(0.0, min(1.0, (val - vmin) / (vmax - vmin)))
current_pos = vmin
for entry in hu_color_map:
# Support both old format ([lower, upper], color) and new LookupTableEntry
if hasattr(entry, 'input_low'):
r_start, r_end = entry.input_low, entry.input_high
color = entry.output_values[0]
else:
r_start, r_end = entry[0]
color = entry[1]
# Skip ranges that don't overlap with the view
if r_end <= vmin or r_start >= vmax:
continue
# Clamp to view bounds
eff_start = max(r_start, vmin)
eff_end = min(r_end, vmax)
# 1. Handle Gap: If there is space between the last processed position and the new range start
if eff_start > current_pos:
norm_prev = get_norm(current_pos)
norm_curr = get_norm(eff_start)
segments.append((norm_prev, trans_white))
segments.append((norm_curr, trans_white))
# 2. Add Color Block
if use_transparency:
c = (color.r, color.g, color.b, color.a)
else:
c = (color.r, color.g, color.b, 1.0)
norm_start = get_norm(eff_start)
norm_end = get_norm(eff_end)
# Add two points to define a solid block of color (gradients handled by creating many small blocks externally)
segments.append((norm_start, c))
segments.append((norm_end, c))
current_pos = eff_end
# 3. Handle trailing Gap (between last range and view max)
if current_pos < vmax:
norm_prev = get_norm(current_pos)
norm_curr = get_norm(vmax)
segments.append((norm_prev, trans_white))
segments.append((norm_curr, trans_white))
# Fallback for completely empty overlap
if not segments:
segments = [(0.0, trans_white), (1.0, trans_white)]
# Create a LinearSegmentedColormap from the list.
# By doubling up coordinates (start, c) and (end, c), we create hard edges if needed,
# or smooth gradients if the input map was generated that way.
cmap = mcolors.LinearSegmentedColormap.from_list("hu_cmap", segments)
# Plot setup
fig, ax = plt.subplots(figsize=(6, 1))
fig.subplots_adjust(bottom=0.5)
norm = plt.Normalize(vmin=vmin, vmax=vmax)
# Use linspace for general readability, as specific map bounds might be irregular
cb = plt.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap),
cax=ax, orientation='horizontal',
ticks=np.linspace(vmin, vmax, 11))
cb.set_label('Hounsfield Units (HU)')
if use_transparency:
# Generate a high-res checkerboard pattern for background
ny, nx = 4, 200
checkers = ((np.indices((ny, nx)).sum(axis=0) % 2) * 0.3 + 0.6)
checkers_rgb = np.stack((checkers,) * 3, axis=-1)
# Display the checkerboard behind the colorbar data
# Extent matches the data range so x-axis aligns
ax.imshow(checkers_rgb, extent=(vmin, vmax, 0, 1),
interpolation='nearest', aspect='auto', zorder=-1)
plt.show()
