synthetic-sciences Clinical Imaging: Clinical & Physiological Imaging Analysis
Overview
Clinical Imaging provides computational tools for analyzing clinical and physiological imaging data. This skill covers diffusion MRI apparent diffusion coefficient (ADC) map computation, micro-CT bone morphometry (BV/TV, trabecular thickness), hemodynamic parameter analysis from blood pressure waveforms, circadian rhythm cosinor analysis, ciliary beat frequency measurement via FFT, tissue deformation analysis using optical flow, and amyloid plaque quantification from fluorescence microscopy.
When to Use This Skill
- Computing ADC maps from multi-b-value diffusion MRI data
- Analyzing bone microarchitecture from micro-CT volumes
- Processing blood pressure or hemodynamic waveform data
- Fitting circadian rhythm data with cosinor models
- Measuring ciliary beat frequency from high-speed video
- Quantifying tissue deformation from image sequences
- Counting and measuring amyloid plaques in fluorescence images
Related Skills: For DICOM file handling use pydicom. For biosignal processing (ECG, EMG, EDA) use neurokit2. For general image analysis use bioimage-analysis.
Installation
uv pip install nibabel SimpleITK scipy opencv-python scikit-image numpy pandas matplotlib
Quick Start
import numpy as np
from scipy.optimize import curve_fit
# Circadian rhythm cosinor analysis
def cosinor(t, mesor, amplitude, acrophase, period=24):
return mesor + amplitude * np.cos(2 * np.pi * t / period + acrophase)
# Example: body temperature over 48 hours
time_hours = np.arange(0, 48, 1)
temperature = 36.8 + 0.3 * np.cos(2 * np.pi * time_hours / 24 - 1.0) + \
np.random.normal(0, 0.1, len(time_hours))
popt, pcov = curve_fit(cosinor, time_hours, temperature,
p0=[36.8, 0.3, -1.0], maxfev=10000)
print(f"MESOR: {popt[0]:.2f} C")
print(f"Amplitude: {popt[1]:.3f} C")
print(f"Acrophase: {np.degrees(popt[2]):.1f} degrees ({popt[2]/(2*np.pi)*24:.1f} h)")
Core Capabilities
1. Diffusion MRI ADC Maps
Compute apparent diffusion coefficient maps from multi-b-value images.
import nibabel as nib
import numpy as np
def compute_adc_map(dwi_paths, b_values):
"""Compute ADC map from multi-b-value DWI images.
ADC = -ln(S/S0) / b, where S0 is the b=0 image.
Args:
dwi_paths: list of NIfTI file paths (one per b-value)
b_values: list of b-values in s/mm²
"""
# Load images
images = [nib.load(p).get_fdata() for p in dwi_paths]
b_values = np.array(b_values)
# Identify b=0 image
b0_idx = np.argmin(b_values)
S0 = images[b0_idx].astype(float)
S0[S0 == 0] = 1e-10 # Avoid division by zero
# Compute ADC using linear regression in log space
# ln(S) = ln(S0) - b * ADC
shape = S0.shape
adc_map = np.zeros(shape)
# Use non-zero b-values
nonzero_mask = b_values > 0
b_nonzero = b_values[nonzero_mask]
for idx, (img, b) in enumerate(zip(images, b_values)):
if b > 0:
ratio = img.astype(float) / S0
ratio = np.clip(ratio, 1e-10, None)
adc_contribution = -np.log(ratio) / b
adc_map += adc_contribution
adc_map /= nonzero_mask.sum()
# Mask out background (where S0 is very low)
tissue_mask = S0 > np.percentile(S0[S0 > 0], 10)
adc_map[~tissue_mask] = 0
# Clip to physiological range (0 - 3.5 x 10^-3 mm²/s)
adc_map = np.clip(adc_map, 0, 3.5e-3)
print(f"ADC map shape: {adc_map.shape}")
print(f"Mean ADC (tissue): {adc_map[tissue_mask].mean()*1e3:.3f} x10⁻³ mm²/s")
print(f"Median ADC (tissue): {np.median(adc_map[tissue_mask])*1e3:.3f} x10⁻³ mm²/s")
return adc_map, tissue_mask
def regional_adc_stats(adc_map, roi_mask):
"""Calculate ADC statistics within a region of interest."""
roi_values = adc_map[roi_mask & (adc_map > 0)]
return {
'mean': roi_values.mean() * 1e3,
'median': np.median(roi_values) * 1e3,
'std': roi_values.std() * 1e3,
'min': roi_values.min() * 1e3,
'max': roi_values.max() * 1e3,
'n_voxels': len(roi_values),
'unit': '10⁻³ mm²/s'
}
2. Micro-CT Bone Morphometry
Quantify bone microarchitecture from 3D micro-CT volumes.
import numpy as np
from scipy import ndimage
import skimage.filters
def bone_morphometry(volume, voxel_size_um, roi_mask=None):
"""Compute trabecular bone morphometric parameters.
Args:
volume: 3D numpy array (micro-CT volume)
voxel_size_um: voxel size in micrometers
roi_mask: optional ROI mask (binary 3D array)
"""
if roi_mask is not None:
vol = volume * roi_mask
else:
vol = volume.copy()
roi_mask = np.ones_like(volume, dtype=bool)
# Segment bone (Otsu thresholding)
threshold = skimage.filters.threshold_otsu(vol[roi_mask])
bone_mask = vol > threshold
# BV/TV: Bone Volume / Total Volume
total_voxels = roi_mask.sum()
bone_voxels = (bone_mask & roi_mask).sum()
bv_tv = bone_voxels / total_voxels
# Trabecular Thickness (Tb.Th) via distance transform
# Mean thickness = 2 * mean distance transform value within bone
bone_distance = ndimage.distance_transform_edt(bone_mask & roi_mask)
tb_th = 2 * bone_distance[bone_mask & roi_mask].mean() * voxel_size_um / 1000 # mm
# Trabecular Spacing (Tb.Sp) via distance transform of marrow
marrow_mask = ~bone_mask & roi_mask
marrow_distance = ndimage.distance_transform_edt(marrow_mask)
tb_sp = 2 * marrow_distance[marrow_mask].mean() * voxel_size_um / 1000 # mm
# Trabecular Number (Tb.N)
tb_n = bv_tv / tb_th if tb_th > 0 else 0 # 1/mm
# Cortical thickness (if applicable)
# Use morphological operations to identify cortex
from skimage.morphology import binary_erosion, ball
eroded = binary_erosion(bone_mask & roi_mask, ball(3))
cortex = (bone_mask & roi_mask) & ~eroded
cortex_dist = ndimage.distance_transform_edt(cortex)
ct_th = cortex_dist[cortex].mean() * voxel_size_um / 1000 if cortex.sum() > 0 else 0
results = {
'BV/TV': bv_tv,
'Tb.Th (mm)': tb_th,
'Tb.Sp (mm)': tb_sp,
'Tb.N (1/mm)': tb_n,
'Ct.Th (mm)': ct_th,
'bone_voxels': bone_voxels,
'total_voxels': total_voxels
}
for key, val in results.items():
if isinstance(val, float):
print(f"{key}: {val:.4f}")
return results
3. Hemodynamic Analysis
Process blood pressure waveforms.
import numpy as np
from scipy.signal import find_peaks
def analyze_blood_pressure(pressure_signal, sampling_rate_hz):
"""Analyze blood pressure waveform.
Args:
pressure_signal: 1D array of pressure values (mmHg)
sampling_rate_hz: sampling frequency
"""
# Detect systolic peaks
min_distance = int(0.5 * sampling_rate_hz) # Min 0.5s between beats
peaks, props = find_peaks(pressure_signal, distance=min_distance,
prominence=20, height=60)
# Detect diastolic troughs
troughs, _ = find_peaks(-pressure_signal, distance=min_distance)
# Systolic and diastolic pressures
systolic = pressure_signal[peaks]
# Match each peak to next trough
diastolic = []
for p in peaks:
next_troughs = troughs[troughs > p]
if len(next_troughs) > 0:
diastolic.append(pressure_signal[next_troughs[0]])
diastolic = np.array(diastolic[:len(systolic)])
# Pulse pressure
pulse_pressure = systolic[:len(diastolic)] - diastolic
# Mean arterial pressure (MAP)
map_pressure = diastolic + pulse_pressure / 3
# Heart rate
rr_intervals = np.diff(peaks) / sampling_rate_hz # seconds
heart_rate = 60 / rr_intervals # BPM
results = {
'systolic_mean': np.mean(systolic),
'systolic_std': np.std(systolic),
'diastolic_mean': np.mean(diastolic),
'diastolic_std': np.std(diastolic),
'pulse_pressure_mean': np.mean(pulse_pressure),
'MAP_mean': np.mean(map_pressure),
'heart_rate_mean': np.mean(heart_rate),
'heart_rate_std': np.std(heart_rate),
'n_beats': len(peaks)
}
print(f"BP: {results['systolic_mean']:.0f}/{results['diastolic_mean']:.0f} mmHg")
print(f"MAP: {results['MAP_mean']:.0f} mmHg")
print(f"HR: {results['heart_rate_mean']:.0f} ± {results['heart_rate_std']:.0f} BPM")
print(f"Pulse pressure: {results['pulse_pressure_mean']:.0f} mmHg")
return results
4. Circadian Rhythm Cosinor
Fit circadian data with cosinor model.
import numpy as np
from scipy.optimize import curve_fit
from scipy import stats
def cosinor_analysis(time_hours, measurements, period=24):
"""Cosinor analysis: Y = MESOR + Amplitude * cos(2*pi*t/T + acrophase).
Args:
time_hours: time points in hours
measurements: observed values
period: assumed period in hours (default: 24)
"""
def cosinor(t, mesor, amplitude, acrophase):
return mesor + amplitude * np.cos(2 * np.pi * t / period + acrophase)
# Initial estimates
mesor_init = np.mean(measurements)
amp_init = (np.max(measurements) - np.min(measurements)) / 2
acro_init = 0
popt, pcov = curve_fit(cosinor, time_hours, measurements,
p0=[mesor_init, amp_init, acro_init],
maxfev=10000)
perr = np.sqrt(np.diag(pcov))
mesor, amplitude, acrophase = popt
# Ensure amplitude is positive
if amplitude < 0:
amplitude = -amplitude
acrophase += np.pi
# Normalize acrophase to [0, 2*pi]
acrophase = acrophase % (2 * np.pi)
# R-squared
predicted = cosinor(time_hours, *popt)
ss_res = np.sum((measurements - predicted) ** 2)
ss_tot = np.sum((measurements - np.mean(measurements)) ** 2)
r_squared = 1 - ss_res / ss_tot
# F-test for rhythm significance
n = len(measurements)
f_stat = (ss_tot - ss_res) / 2 / (ss_res / (n - 3))
p_value = 1 - stats.f.cdf(f_stat, 2, n - 3)
# Convert acrophase to clock time
acrophase_hours = (-acrophase / (2 * np.pi) * period) % period
acro_h = int(acrophase_hours)
acro_m = int((acrophase_hours - acro_h) * 60)
results = {
'MESOR': mesor,
'amplitude': amplitude,
'acrophase_rad': acrophase,
'acrophase_hours': acrophase_hours,
'acrophase_clock': f"{acro_h:02d}:{acro_m:02d}",
'period': period,
'r_squared': r_squared,
'f_statistic': f_stat,
'p_value': p_value,
'significant': p_value < 0.05
}
print(f"MESOR: {mesor:.3f} ± {perr[0]:.3f}")
print(f"Amplitude: {amplitude:.3f} ± {perr[1]:.3f}")
print(f"Acrophase: {results['acrophase_clock']} ({acrophase_hours:.1f} h)")
print(f"R²: {r_squared:.4f}")
print(f"Rhythm p-value: {p_value:.2e} ({'significant' if p_value < 0.05 else 'not significant'})")
return results
5. Ciliary Beat Frequency
Measure CBF from high-speed video using FFT.
import numpy as np
import cv2
def measure_ciliary_beat_frequency(video_path, fps, roi=None):
"""Measure ciliary beat frequency from high-speed video via FFT.
Args:
video_path: path to video file
fps: frames per second of recording
roi: (x, y, w, h) region of interest tuple, or None for full frame
"""
cap = cv2.VideoCapture(video_path)
frames = []
while True:
ret, frame = cap.read()
if not ret:
break
gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
if roi:
x, y, w, h = roi
gray = gray[y:y+h, x:x+w]
frames.append(gray.astype(float))
cap.release()
frames = np.array(frames)
n_frames = len(frames)
print(f"Loaded {n_frames} frames at {fps} fps ({n_frames/fps:.1f} seconds)")
# Mean intensity time series per pixel
mean_intensity = frames.mean(axis=(1, 2))
# FFT of mean intensity
fft_vals = np.fft.rfft(mean_intensity - mean_intensity.mean())
freqs = np.fft.rfftfreq(n_frames, d=1/fps)
power = np.abs(fft_vals) ** 2
# Find dominant frequency (excluding DC and very low frequencies)
min_freq_idx = max(1, int(2 * n_frames / fps)) # Ignore < 2 Hz
max_freq_idx = len(freqs) - 1
dominant_idx = min_freq_idx + np.argmax(power[min_freq_idx:max_freq_idx])
cbf = freqs[dominant_idx]
print(f"Dominant ciliary beat frequency: {cbf:.1f} Hz")
print(f"Period: {1/cbf*1000:.1f} ms")
# Spatiotemporal frequency map
cbf_map = np.zeros(frames[0].shape)
for y in range(0, frames.shape[1], 4):
for x in range(0, frames.shape[2], 4):
pixel_ts = frames[:, y, x]
pixel_fft = np.fft.rfft(pixel_ts - pixel_ts.mean())
pixel_power = np.abs(pixel_fft) ** 2
peak_idx = min_freq_idx + np.argmax(pixel_power[min_freq_idx:max_freq_idx])
cbf_map[y, x] = freqs[peak_idx]
return cbf, cbf_map, freqs, power
6. Tissue Deformation Flow
Quantify tissue deformation using optical flow.
import cv2
import numpy as np
def analyze_tissue_deformation(frame1, frame2, method='lucas_kanade'):
"""Compute optical flow between two tissue images.
Args:
frame1, frame2: consecutive grayscale frames
method: 'lucas_kanade' or 'farneback'
"""
if method == 'farneback':
flow = cv2.calcOpticalFlowFarneback(
frame1, frame2, None,
pyr_scale=0.5, levels=3, winsize=15,
iterations=3, poly_n=5, poly_sigma=1.2, flags=0
)
else: # Lucas-Kanade (sparse)
pts = cv2.goodFeaturesToTrack(frame1, maxCorners=500,
qualityLevel=0.01, minDistance=10)
if pts is None:
return None
pts_new, status, err = cv2.calcOpticalFlowPyrLK(frame1, frame2, pts, None)
good_old = pts[status.ravel() == 1]
good_new = pts_new[status.ravel() == 1]
# Compute displacement field
dx = good_new[:, 0] - good_old[:, 0]
dy = good_new[:, 1] - good_old[:, 1]
magnitude = np.sqrt(dx**2 + dy**2)
print(f"Tracked {len(good_old)} points")
print(f"Mean displacement: {magnitude.mean():.2f} px")
print(f"Max displacement: {magnitude.max():.2f} px")
return {'dx': dx, 'dy': dy, 'magnitude': magnitude,
'old_pts': good_old, 'new_pts': good_new}
# For dense flow (Farneback)
u = flow[:, :, 0]
v = flow[:, :, 1]
magnitude = np.sqrt(u**2 + v**2)
# Divergence (expansion/contraction)
du_dx = np.gradient(u, axis=1)
dv_dy = np.gradient(v, axis=0)
divergence = du_dx + dv_dy
# Curl (rotation)
du_dy = np.gradient(u, axis=0)
dv_dx = np.gradient(v, axis=1)
curl = dv_dx - du_dy
# Strain tensor components
exx = du_dx
eyy = dv_dy
exy = 0.5 * (du_dy + dv_dx)
print(f"Mean displacement: {magnitude.mean():.3f} px")
print(f"Mean divergence: {divergence.mean():.6f}")
print(f"Mean curl: {curl.mean():.6f}")
return {
'flow': flow, 'magnitude': magnitude,
'divergence': divergence, 'curl': curl,
'strain_xx': exx, 'strain_yy': eyy, 'strain_xy': exy
}
7. Amyloid Plaque Quantification
Segment and measure amyloid plaques from fluorescence images.
import numpy as np
import skimage.io
import skimage.filters
import skimage.measure
from skimage.morphology import remove_small_objects
def quantify_amyloid_plaques(image_path, min_plaque_area=50):
"""Quantify amyloid plaques from fluorescence microscopy.
Args:
image_path: path to fluorescence image (ThT/Congo Red staining)
min_plaque_area: minimum plaque area in pixels
"""
image = skimage.io.imread(image_path)
if image.ndim == 3:
image = image[:, :, 0] # Take first channel
# Preprocessing
smoothed = skimage.filters.gaussian(image, sigma=2)
# Threshold
thresh = skimage.filters.threshold_otsu(smoothed)
binary = smoothed > thresh
# Remove small objects (noise)
cleaned = remove_small_objects(binary, min_size=min_plaque_area)
# Label and measure
labels = skimage.measure.label(cleaned)
regions = skimage.measure.regionprops(labels, intensity_image=image)
# Extract measurements
plaques = []
for r in regions:
plaques.append({
'label': r.label,
'area_px': r.area,
'perimeter': r.perimeter,
'eccentricity': r.eccentricity,
'mean_intensity': r.mean_intensity,
'max_intensity': r.max_intensity,
'centroid_y': r.centroid[0],
'centroid_x': r.centroid[1]
})
import pandas as pd
df = pd.DataFrame(plaques)
# Summary
total_area = image.shape[0] * image.shape[1]
plaque_area = df['area_px'].sum()
plaque_density = len(df) / (total_area / 1e6) # per million pixels
print(f"Plaque count: {len(df)}")
print(f"Total plaque area: {plaque_area} px ({100*plaque_area/total_area:.2f}%)")
print(f"Mean plaque area: {df['area_px'].mean():.0f} px")
print(f"Plaque density: {plaque_density:.1f} per Mpx")
return df
Typical Workflows
Workflow 1: Compute ADC Map from Multi-b-Value Diffusion MRI
adc_map, mask = compute_adc_map(
dwi_paths=['b0.nii.gz', 'b500.nii.gz', 'b1000.nii.gz'],
b_values=[0, 500, 1000]
)
stats = regional_adc_stats(adc_map, mask)
print(f"Mean ADC: {stats['mean']:.3f} x10⁻³ mm²/s")
Workflow 2: Analyze Bone Microarchitecture from Micro-CT
import tifffile
volume = tifffile.imread('microct_stack.tif')
results = bone_morphometry(volume, voxel_size_um=10)
print(f"BV/TV: {results['BV/TV']:.4f}")
print(f"Tb.Th: {results['Tb.Th (mm)']:.4f} mm")
Workflow 3: Fit Circadian Rhythm Data with Cosinor Analysis
import numpy as np
time = np.arange(0, 48, 2) # Every 2 hours for 48 hours
activity = 100 + 40 * np.cos(2 * np.pi * time / 24 - 0.5) + np.random.normal(0, 10, len(time))
results = cosinor_analysis(time, activity)
print(f"Peak activity at: {results['acrophase_clock']}")
Best Practices
- ADC maps — use at least 2 non-zero b-values for reliable fitting; b=0 and b=1000 s/mm² is standard for brain
- Bone morphometry — validate Otsu threshold against manual segmentation; report voxel size and scan parameters
- Cosinor analysis — collect data spanning at least 2 full cycles; report MESOR, amplitude, acrophase, and p-value
- Ciliary beat — recording must be at least 2x the expected CBF (Nyquist); typical CBF is 8-15 Hz requiring >30 fps
- Optical flow — use Farneback for dense deformation fields, Lucas-Kanade for sparse tracking of features
- Always report units — ADC in 10^-3 mm²/s, bone metrics in mm, blood pressure in mmHg, CBF in Hz
Troubleshooting
Problem: ADC map has noisy regions Solution: Apply Gaussian smoothing to DWI images before ADC calculation. Increase b-value range for better SNR. Mask out background using S0 threshold.
Problem: Bone segmentation includes soft tissue Solution: Apply manual ROI selection before thresholding. Use adaptive thresholding for inhomogeneous CT values. Verify Hounsfield unit calibration.
Problem: Cosinor analysis not significant Solution: Insufficient data points or too much noise. Collect more timepoints. Try different period values if 24h is not assumed. Check for ultradian rhythms.
Problem: CBF measurement gives wrong frequency Solution: Verify recording fps is correct. Check ROI contains actively beating cilia. Increase recording duration for better frequency resolution.