Untitled
unknown
plain_text
10 months ago
75 kB
17
Indexable
#!/usr/bin/env python3
"""
SLM, SLE & SLIP-Ω: Silence-Lip metrics for lip-sync videos.
This script extends your original SLM (SLM-Delta, SLMS-Beta) and SLE (twitch rate, parted time)
with a *novel*, publication-ready composite metric:
SLIP-Ω (Silent-Lip Instability & Persistence Composite)
SLIP-Ω fuses:
* Multi-order kinematics (velocity, acceleration, jerk)
* Event statistics (twitch rate, burstiness, velocity zero-crossing)
* Persistence morphology (parted time ratio & run-length Gini)
* Temporal complexity & spectrum (sample entropy, HF velocity energy, total variation)
* Robust dataset-normalized scoring (higher values = better lip-sync quality)
* Expert-designed weighting based on lip-sync quality importance
IMPORTANT: SLIP-Ω interpretation:
- Higher SLIP-Ω values indicate BETTER lip-sync quality
- Good lip-sync videos should have HIGH SLIP-Ω scores (minimal motion during silence)
- Poor lip-sync videos will have LOW SLIP-Ω scores (excessive motion, open mouth during silence)
Outputs:
- Global SLM, SLE and SLIP-Ω metrics (per video)
- Per-silence-segment metrics (optional CSV)
- Optional debug overlay video (showing silent segments)
- NEW: Optional landmark overlay video (showing detected lip landmarks)
- NEW: Optional landmark-in-silence overlay video (showing landmarks only during silent segments)
Author: Your Name (as per original script)
"""
import torch
import argparse
import os
import sys
import json
import math
from dataclasses import dataclass, asdict
from typing import List, Optional, Tuple, Dict
import numpy as np
import pandas as pd
import cv2
import mediapipe as mp
try:
# Preferred import path
from moviepy.editor import VideoFileClip
except Exception:
# Fallback for environments exposing VideoFileClip at top-level
from moviepy import VideoFileClip
import librosa
from tqdm import tqdm
from scipy.signal import find_peaks, butter, filtfilt
# ------------------------------------------------------------
# Constants: MediaPipe Face Mesh landmark indices for lips
# ------------------------------------------------------------
UPPER_LIP_CENTER = 13
LOWER_LIP_CENTER = 14
LIP_LEFT_CORNER = 78
LIP_RIGHT_CORNER = 308
# Define the specific landmarks for the mouth that we want to visualize and connect
LIP_LANDMARK_INDICES = [
UPPER_LIP_CENTER,
LOWER_LIP_CENTER,
LIP_LEFT_CORNER,
LIP_RIGHT_CORNER,
]
# Order of points to draw a closed polygon around the mouth area
LIP_POLYGON_POINTS_ORDER = [
UPPER_LIP_CENTER,
LIP_RIGHT_CORNER,
LOWER_LIP_CENTER,
LIP_LEFT_CORNER,
]
# Padding ratio for face crop: 0.4 means 40% padding on each side (top, bottom, left, right)
# So, a 100x100 face would be cropped into a (100+40+40)x(100+40+40) = 180x180 region.
FACE_CROP_PADDING_RATIO = 0.4
# ------------------------------------------------------------
# Data containers
# ------------------------------------------------------------
@dataclass
class SilenceSegment:
start_frame: int
end_frame: int # inclusive
start_time: float
end_time: float
@property
def length_frames(self) -> int:
return self.end_frame - self.start_frame + 1
@property
def length_sec(self) -> float:
return self.end_time - self.start_time
@dataclass
class SegmentMetrics:
segment_index: int
start_time: float
end_time: float
num_frames: int
# Original SLM components
vrms: float
arms: float
p_idle: float
slm_delta: float
slms_beta: float
# NEW: SLE components
num_twitches: int
sle_twitch_rate: float
sle_parted_time_ratio: float
# NEW: SLIP-Ω components
jrms: float
zcr_v: float
burstiness: float
sampen: float
hf_energy: float
tv1: float
gini_runs: float
micro_motion: float
slip_omega: float
# ------------------------------------------------------------
# Audio → silence mask
# ------------------------------------------------------------
def compute_silence_mask(
video_path: str,
rms_db_thresh: float = -45.0, # ignored (kept for signature compatibility)
hop_length: int = 512, # ignored (kept for signature compatibility)
min_silence_frames: int = 2,
vad_threshold: float = 0.5, # Silero VAD speech prob threshold (0–1)
target_sr: int = 16000, # Silero works best at 16 kHz
use_onnx: bool = False, # set True if you prefer ONNX runtime weights
) -> Tuple[np.ndarray, float, int]:
"""
Computes a boolean mask of *silent* frames using Silero VAD (robust to background noise).
Returns:
silent_frames: np.ndarray[bool] length ~= #video frames
fps: float (video FPS)
frame_count: int (#video frames)
"""
# --- Load video / audio ---
clip = VideoFileClip(video_path, audio=True)
fps = clip.fps
sr_src = clip.audio.fps
# mono audio as float32
audio = clip.audio.to_soundarray(fps=sr_src).astype(np.float32)
if audio.ndim == 2:
audio = audio.mean(axis=1)
# resample to target_sr for Silero
if sr_src != target_sr:
audio = librosa.resample(audio, orig_sr=sr_src, target_sr=target_sr)
sr = target_sr
else:
sr = sr_src
# --- Load Silero VAD ---
# First call will download weights (internet required once).
model, utils = torch.hub.load(
repo_or_dir='snakers4/silero-vad',
model='silero_vad',
onnx=use_onnx,
trust_repo=True
)
get_speech_timestamps, *_ = utils
# to torch tensor
wav = torch.from_numpy(audio)
# --- Run VAD (returns sample indices for speech) ---
speech_ts = get_speech_timestamps(
wav, model, sampling_rate=sr, threshold=vad_threshold
)
# --- Map speech segments to frame-level mask ---
frame_count_est = int(round(clip.duration * fps))
speech_frames = np.zeros(frame_count_est, dtype=bool)
for ts in speech_ts:
# ts['start'] / ts['end'] are in *samples*
start_t = ts['start'] / sr
end_t = ts['end'] / sr
start_f = max(0, int(math.floor(start_t * fps)))
end_f = min(frame_count_est - 1, int(math.ceil(end_t * fps)) - 1)
if end_f >= start_f:
speech_frames[start_f:end_f+1] = True
# Optional: remove ultra-short speech blips and silence islands
# Reuse your existing utility for run-length cleanup
min_speech_frames = max(1, min_silence_frames) # use same window by default
if min_speech_frames > 1:
speech_frames = _remove_short_runs(speech_frames, min_speech_frames)
silent_frames = ~speech_frames
if min_silence_frames > 1:
silent_frames = _remove_short_runs(silent_frames, min_silence_frames)
clip.close()
return silent_frames, fps, frame_count_est
def _remove_short_runs(mask: np.ndarray, min_len: int) -> np.ndarray:
"""Removes short contiguous runs of True values from a boolean mask."""
out = mask.copy()
n = len(mask)
i = 0
while i < n:
if mask[i]: # If current element is True (silent)
j = i
while j < n and mask[j]: # Find end of current run
j += 1
if j - i < min_len: # If run is shorter than min_len, set to False
out[i:j] = False
i = j # Move cursor to end of current run
else:
i += 1 # Move cursor to next element
return out
# ------------------------------------------------------------
# Silence segments from mask
# ------------------------------------------------------------
def mask_to_segments(silent_frames: np.ndarray, fps: float) -> List[SilenceSegment]:
"""Converts a boolean silence mask into a list of SilenceSegment objects."""
segs: List[SilenceSegment] = []
n = len(silent_frames)
i = 0
while i < n:
if silent_frames[i]:
j = i
while j < n and silent_frames[j]:
j += 1
start_f = i
end_f = j - 1
start_t = start_f / fps
end_t = (end_f + 1) / fps
segs.append(SilenceSegment(start_frame=start_f, end_frame=end_f, start_time=start_t, end_time=end_t))
i = j
else:
i += 1
return segs
# ------------------------------------------------------------
# Landmark tracker with caching
# ------------------------------------------------------------
class LipLandmarkCache:
"""Caches extracted lip landmark coordinates for efficient access."""
def __init__(self, total_frames: int):
# Stores x,y for UPPER_LIP_CENTER, LOWER_LIP_CENTER, LIP_LEFT_CORNER, LIP_RIGHT_CORNER
# Total 4 points * 2 coords = 8 values per frame.
self.data = np.full((total_frames, 8), np.nan, dtype=np.float32)
def set(self, frame_idx: int, vals: np.ndarray):
if 0 <= frame_idx < self.data.shape[0]:
self.data[frame_idx] = vals
def get(self, frame_idx: int) -> Optional[np.ndarray]:
if not (0 <= frame_idx < self.data.shape[0]):
return None
v = self.data[frame_idx]
if np.isnan(v).all():
return None
return v
def as_array(self) -> np.ndarray:
return self.data
def extract_lip_landmarks(
video_path: str,
total_frames: int,
target_mask: Optional[np.ndarray] = None, # Only cache landmarks for frames where mask is True
progress: bool = True,
landmark_overlay_path: Optional[str] = None, # Path to save video with landmark overlay
) -> LipLandmarkCache:
"""
Extracts specified lip landmarks from a video using MediaPipe Face Mesh.
Includes face detection and cropping for robust landmark extraction.
Optionally caches landmarks for frames specified by `target_mask` and
generates an overlay video showing detected landmarks.
"""
cache = LipLandmarkCache(total_frames)
# Initialize MediaPipe Face Detection (for robust initial face finding)
mp_face_detection = mp.solutions.face_detection
face_detection = mp_face_detection.FaceDetection(model_selection=1, min_detection_confidence=0.7) # model_selection=1 for full-range faces
# Initialize MediaPipe Face Mesh (for detailed lip landmarks)
mp_face_mesh = mp.solutions.face_mesh
face_mesh = mp_face_mesh.FaceMesh(
static_image_mode=False, refine_landmarks=True, max_num_faces=1
)
cap = cv2.VideoCapture(video_path)
if not cap.isOpened():
raise RuntimeError(f"Cannot open video: {video_path}")
cap_fps = cap.get(cv2.CAP_PROP_FPS)
cap_width = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
cap_height = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))
writer = None
if landmark_overlay_path:
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter(landmark_overlay_path, fourcc, cap_fps, (cap_width, cap_height))
if not writer.isOpened():
print(f"Warning: Could not open video writer for {landmark_overlay_path}. Skipping landmark overlay.")
writer = None
else:
print(f" - Generating landmark overlay video: {landmark_overlay_path}")
iterator_limit = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) if writer else total_frames
iterator = tqdm(range(iterator_limit), desc="Landmarks", unit="f") if progress else range(iterator_limit)
for _ in iterator:
ret, frame = cap.read()
if not ret:
break
idx = int(cap.get(cv2.CAP_PROP_POS_FRAMES)) - 1
if idx >= iterator_limit:
break
rgb_frame = cv2.cvtColor(frame, cv2.COLOR_BGR2RGB)
# --- Stage 1: Robust Face Detection ---
detection_results = face_detection.process(rgb_frame)
# Default to no landmark detection for current frame
current_frame_lip_pts_normalized = np.full(8, np.nan, dtype=np.float32)
if detection_results.detections:
# Assume single face for simplicity, pick the first one
detection = detection_results.detections[0]
bbox_norm = detection.location_data.relative_bounding_box
# Convert normalized bbox to pixel coordinates
xmin = int(bbox_norm.xmin * cap_width)
ymin = int(bbox_norm.ymin * cap_height)
width = int(bbox_norm.width * cap_width)
height = int(bbox_norm.height * cap_height)
# Calculate padding and new crop coordinates
pad_w = int(width * FACE_CROP_PADDING_RATIO)
pad_h = int(height * FACE_CROP_PADDING_RATIO)
x1 = max(0, xmin - pad_w)
y1 = max(0, ymin - pad_h)
x2 = min(cap_width, xmin + width + pad_w)
y2 = min(cap_height, ymin + height + pad_h)
# Crop the frame
cropped_rgb_frame = rgb_frame[y1:y2, x1:x2]
if cropped_rgb_frame.shape[0] > 0 and cropped_rgb_frame.shape[1] > 0:
# --- Stage 2: Detailed Lip Landmark Extraction on Cropped Face ---
mesh_results = face_mesh.process(cropped_rgb_frame)
if mesh_results.multi_face_landmarks:
lm = mesh_results.multi_face_landmarks[0].landmark
# Store normalized coordinates
# IMPORTANT: These are normalized to the *cropped frame*.
# We need to convert them back to the original frame's normalized coordinates.
# Current cropped frame dimensions
cropped_h, cropped_w, _ = cropped_rgb_frame.shape
# Convert from normalized_cropped to pixel_cropped:
pixel_x_cropped = lm[UPPER_LIP_CENTER].x * cropped_w
pixel_y_cropped = lm[UPPER_LIP_CENTER].y * cropped_h
# Convert from pixel_cropped to pixel_original:
pixel_x_orig_upper = pixel_x_cropped + x1
pixel_y_orig_upper = pixel_y_cropped + y1
pixel_x_cropped = lm[LOWER_LIP_CENTER].x * cropped_w
pixel_y_cropped = lm[LOWER_LIP_CENTER].y * cropped_h
pixel_x_orig_lower = pixel_x_cropped + x1
pixel_y_orig_lower = pixel_y_cropped + y1
pixel_x_cropped = lm[LIP_LEFT_CORNER].x * cropped_w
pixel_y_cropped = lm[LIP_LEFT_CORNER].y * cropped_h
pixel_x_orig_left = pixel_x_cropped + x1
pixel_y_orig_left = pixel_y_cropped + y1
pixel_x_cropped = lm[LIP_RIGHT_CORNER].x * cropped_w
pixel_y_cropped = lm[LIP_RIGHT_CORNER].y * cropped_h
pixel_x_orig_right = pixel_x_cropped + x1
pixel_y_orig_right = pixel_y_cropped + y1
current_frame_lip_pts_normalized = np.array(
[
pixel_x_orig_upper / cap_width, pixel_y_orig_upper / cap_height,
pixel_x_orig_lower / cap_width, pixel_y_orig_lower / cap_height,
pixel_x_orig_left / cap_width, pixel_y_orig_left / cap_height,
pixel_x_orig_right / cap_width, pixel_y_orig_right / cap_height,
],
dtype=np.float32,
)
# Decide if landmarks for this frame should be stored in cache
store_to_cache = True
# Store if target_mask specifies, OR if any overlay that uses cached data is requested
if target_mask is not None and not landmark_overlay_path: # ONLY restrict if *only* target mask and no overlay is requested
if idx >= len(target_mask) or not target_mask[idx]:
store_to_cache = False
if store_to_cache:
cache.set(idx, current_frame_lip_pts_normalized)
# --- Drawing logic for landmark overlay ---
if writer:
img_h, img_w, _ = frame.shape
# Only draw if landmarks were successfully detected for this frame
if not np.any(np.isnan(current_frame_lip_pts_normalized)):
# Reshape current_frame_lip_pts_normalized back to (N, 2) structure
# and scale to pixel coordinates for drawing
# Create a temporary array mapping our stored order to drawing order
# The stored order is: Upper, Lower, Left, Right
# The drawing polygon order is: Upper, Right, Lower, Left
stored_pts_xy = current_frame_lip_pts_normalized.reshape(-1, 2) # (4, 2)
# Map stored indices to polygon order
# LIP_POLYGON_POINTS_ORDER = [UPPER_LIP_CENTER, LIP_RIGHT_CORNER, LOWER_LIP_CENTER, LIP_LEFT_CORNER]
# Stored order indices: 0=Upper, 1=Lower, 2=Left, 3=Right
# So for drawing polygon:
# Upper (stored idx 0), Right (stored idx 3), Lower (stored idx 1), Left (stored idx 2)
ordered_pts_for_poly_norm = np.array([
stored_pts_xy[0], # Upper
stored_pts_xy[3], # Right
stored_pts_xy[1], # Lower
stored_pts_xy[2], # Left
])
# Scale to pixel coordinates
ordered_pts_for_poly_scaled = (ordered_pts_for_poly_norm * np.array([img_w, img_h])).astype(int)
# Draw circles for the key lip points (all 4 of the original ones)
original_lm_indices_order = [UPPER_LIP_CENTER, LOWER_LIP_CENTER, LIP_LEFT_CORNER, LIP_RIGHT_CORNER]
for i, lm_idx in enumerate(original_lm_indices_order):
center = tuple( (stored_pts_xy[i] * np.array([img_w, img_h])).astype(int) )
cv2.circle(frame, center, 3, (0, 255, 0), -1) # Green circle, filled
# Draw lines to connect the lip region, forming a closed polygon
pts_cv = ordered_pts_for_poly_scaled.reshape((-1, 1, 2))
cv2.polylines(frame, [pts_cv], True, (255, 0, 0), 2) # Blue lines, closed, 2px thickness
# Always write frame to video if writer is active, even if no landmarks detected
# (in which case the original frame is written unmodified, or with face bbox if only face detected).
if writer:
writer.write(frame)
cap.release()
face_detection.close() # Close MediaPipe models
face_mesh.close()
if writer:
writer.release()
print(" - Landmark overlay video generation complete.")
return cache
# ------------------------------------------------------------
# Mouth opening utility
# ------------------------------------------------------------
def opening_from_pts(vec: np.ndarray) -> float:
"""
Calculates mouth opening ratio based on key lip landmarks.
Expected `vec` order: (upper_x, upper_y, lower_x, lower_y, left_x, left_y, right_x, right_y)
"""
# Unpack coordinates for readability
xu, yu, xl, yl, xL, yL, xR, yR = vec
if np.any(np.isnan(vec)): # Handle cases where landmarks were not detected
return np.nan
# Vertical distance between upper and lower lip centers
gap = math.hypot(xu - xl, yu - yl)
# Horizontal distance between left and right lip corners
width = math.hypot(xL - xR, yL - yR)
if width <= 0: # Avoid division by zero
return np.nan
return gap / width # Ratio of gap to width
def opening_series(cache: LipLandmarkCache) -> np.ndarray:
"""Converts a landmark cache into a time series of mouth opening values."""
arr = cache.as_array()
# Apply opening_from_pts function row-wise to the landmark array
o = np.apply_along_axis(opening_from_pts, 1, arr)
return o.astype(np.float32)
# ------------------------------------------------------------
# SLM Metrics (Original)
# ------------------------------------------------------------
def slm_delta(openings: np.ndarray, fps: float, mask: np.ndarray) -> float:
"""
Calculates SLM-Delta: RMS of mouth opening velocity during silent periods.
"""
idx = np.where(mask)[0] # Get indices of silent frames
o = openings[idx]
o = o[~np.isnan(o)] # Remove NaNs from silent opening values
if len(o) < 2: # Need at least two points to calculate velocity
return 0.0
v = fps * np.diff(o) # Velocity = change in opening * FPS
return float(np.sqrt(np.mean(v ** 2))) # RMS velocity
def slms_beta(
openings: np.ndarray,
fps: float,
mask: np.ndarray,
tau: float = 0.20, # Threshold for mouth being 'parted'
lam: float = 0.4, # Weight for acceleration RMS
mu: float = 0.3, # Weight for proportion of idle time
) -> float:
"""
Calculates SLMS-Beta: A weighted sum of RMS velocity, RMS acceleration,
and the proportion of time the mouth is idle (parted).
"""
idx = np.where(mask)[0]
o = openings[idx]
o = o[~np.isnan(o)]
if len(o) < 1:
return 0.0
vrms = 0.0
arms = 0.0
p_idle = 0.0
if len(o) >= 2:
v = fps * np.diff(o)
vrms = math.sqrt(np.mean(v ** 2))
if len(o) >= 3: # Need at least three points for 2nd derivative (acceleration)
a = (fps ** 2) * np.diff(o, n=2)
arms = math.sqrt(np.mean(a ** 2))
p_idle = float((o > tau).mean()) if len(o) > 0 else 0.0
return float(vrms + lam * arms + mu * p_idle)
# ------------------------------------------------------------
# SLE: event-based
# ------------------------------------------------------------
def _smooth(o_clean: np.ndarray, window: int) -> np.ndarray:
"""Applies a rolling mean smoothing filter to a 1D array."""
if window <= 1:
return o_clean
if window % 2 == 0: # Ensure window is odd for centered smoothing
window += 1
# Pad the signal to handle edges correctly using 'edge' mode
pad_width = window // 2
if len(o_clean) <= pad_width * 2: # If signal is too short to be padded meaningfully
return o_clean # or handle this case appropriately, e.g., return constant
padded_o = np.pad(o_clean, (pad_width, pad_width), mode='edge')
# Use convolution for a simple moving average
kernel = np.ones(window) / window
s = np.convolve(padded_o, kernel, mode='valid')
return s
def _find_twitch_peaks(
openings_segment: np.ndarray,
fps: float,
min_twitch_height: float,
min_twitch_prominence: float,
smoothing_window_frames: int,
) -> Tuple[np.ndarray, np.ndarray]:
"""
Helper function to find peaks (twitches) in the mouth opening signal
after optional smoothing. Returns peak indices and the cleaned signal.
"""
o_clean = openings_segment[~np.isnan(openings_segment)]
if len(o_clean) < 5: # Need a minimum length for smoothing and peak finding
return np.array([]), o_clean
# Apply smoothing before peak detection
smoothed_o = _smooth(o_clean, smoothing_window_frames)
# Ensure smoothed_o is not empty or all NaNs/Infs before peak finding
if not np.any(np.isfinite(smoothed_o)):
return np.array([]), o_clean
peaks, _ = find_peaks(
smoothed_o,
height=min_twitch_height,
prominence=min_twitch_prominence
)
return peaks, o_clean
def detect_twitches(
openings_segment: np.ndarray,
fps: float,
min_twitch_height: float = 0.1,
min_twitch_prominence: float = 0.05,
smoothing_window_frames: int = 3,
) -> int:
"""
Detects the number of 'twitches' (significant opening peaks)
in a segment of mouth opening data.
"""
peaks, _ = _find_twitch_peaks(
openings_segment,
fps,
min_twitch_height=min_twitch_height,
min_twitch_prominence=min_twitch_prominence,
smoothing_window_frames=smoothing_window_frames,
)
return len(peaks)
def calculate_parted_time_ratio(
openings_segment: np.ndarray,
tau: float # Threshold for mouth being 'parted'
) -> float:
"""
Calculates the proportion of time the mouth is 'parted' (open > tau)
within a given segment.
"""
o_clean = openings_segment[~np.isnan(openings_segment)]
if len(o_clean) == 0:
return 0.0
parted_frames = np.sum(o_clean > tau)
total_frames = len(o_clean)
return float(parted_frames / total_frames)
# ------------------------------------------------------------
# SLIP-Ω utilities
# ------------------------------------------------------------
def deriv(seq: np.ndarray, order: int = 1, fps: float = 1.0) -> np.ndarray:
"""
Calculates the N-th order numerical derivative of a sequence.
Handles NaNs by ignoring them for diff and scales by FPS.
"""
s = seq[~np.isnan(seq)]
if len(s) <= order: # Need at least order+1 points for Nth derivative
return np.array([])
d = np.diff(s, n=order)
return d * (fps ** order)
def zero_crossing_rate(x: np.ndarray) -> float:
"""
Calculates the zero-crossing rate of a 1D signal.
Counts sign changes.
"""
x = x[~np.isnan(x)]
if len(x) < 2:
return 0.0
signs = np.sign(x)
# Count where sign changes (signs[t] * signs[t-1] < 0)
# Sum these occurrences and divide by total possible pairs (len(x)-1)
return np.mean(signs[1:] * signs[:-1] < 0)
def sample_entropy(x: np.ndarray, m: int = 2, r: float = None) -> float:
"""
Calculates the Sample Entropy of a 1D time series.
Measures the complexity/irregularity of a signal.
"""
x = x[~np.isnan(x)]
N = len(x)
if N <= m + 1: # Minimum length for m=2 is N=3.
return 0.0
if r is None:
# Default r: 0.2 times the standard deviation of the data
std_x = np.std(x)
r = 0.2 * std_x
if r == 0: # Avoid division by zero if data is constant
return 0.0
def _phi(m_):
# Counts pairs of sequences of length m_ that are similar (within 'r' distance)
B = 0
# Iterate through all possible starting points for the first sequence
for i in range(N - m_ + 1):
# Iterate through all possible starting points for the second sequence (j > i)
for j in range(i + 1, N - m_ + 1):
# Check if the maximum difference between corresponding elements in the two sequences is <= r
if np.max(np.abs(x[i : i + m_] - x[j : j + m_])) <= r:
B += 1
# Normalize the count by the number of possible pairs (N-m_)(N-m_-1)/2
# Here, it's normalized by (N-m_)*(N-m_-1) as the sum iterates over all pairs (i, j) with i != j
# The (N-m_)*(N-m_-1) is total number of pairs of segments of length m_
denominator = (N - m_) * (N - m_ - 1)
return B / denominator if denominator > 0 else 0.0
a_val = _phi(m)
b_val = _phi(m + 1)
# If either probability is zero, entropy is undefined or infinite. Return 0 for practical purposes.
if a_val == 0 or b_val == 0:
return 0.0
# Sample Entropy formula: -ln(A/B) where A is #matches for m+1, B for m.
# In my _phi, B is count / (N-m)(N-m-1). So it's already a frequency.
# SampEn = -log(prob_match_m+1 / prob_match_m)
return -np.log(b_val / a_val)
def hf_energy(v: np.ndarray, fps: float, f_low=0.5, f_high=5.0, order=4) -> float:
"""
Calculates the high-frequency RMS energy of the velocity signal.
Applies a Butterworth bandpass filter and then computes RMS.
"""
v = v[~np.isnan(v)]
if len(v) < 2:
return 0.0
nyq = 0.5 * fps # Nyquist frequency
low = max(f_low / nyq, 1e-6) # Normalized low cut-off frequency
high = min(f_high / nyq, 0.999999) # Normalized high cut-off frequency
if high <= low: # Invalid bandpass range
return 0.0
try:
b, a = butter(order, [low, high], btype='band')
# filtfilt requires the signal length to be greater than 3 times (max(len(a), len(b)) - 1)
# For order 4, len(a)=len(b)=5, so padlen_min = 3*(5-1) = 12.
padlen = max(3 * (max(len(a), len(b)) - 1), 0)
if len(v) <= padlen: # Signal too short for filter
return 0.0
vf = filtfilt(b, a, v, padlen=padlen)
except ValueError as e:
# Catch specific ValueError from filtfilt for too short signal or other filter issues
# print(f"Warning: filtfilt failed for hf_energy ({e}). Returning 0.0. Signal length: {len(v)}", file=sys.stderr)
return 0.0
except Exception as e:
# Catch any other unexpected errors during filtering
print(f"Warning: hf_energy calculation failed ({e}). Returning 0.0.", file=sys.stderr)
return 0.0
return float(np.sqrt(np.mean(vf ** 2)))
def total_variation(o: np.ndarray, fps: float) -> float:
"""
Calculates the total variation of the mouth opening signal.
Defined as the sum of absolute differences between consecutive points, normalized by duration.
"""
o = o[~np.isnan(o)]
if len(o) < 2:
return 0.0
# Mean of absolute differences, scaled by FPS to represent a "rate of variation"
return float(np.mean(np.abs(np.diff(o))) * fps)
def gini_coefficient(x: np.ndarray) -> float:
"""
Calculates the Gini coefficient for a 1D array of non-negative values.
A measure of statistical dispersion, commonly used for inequality.
"""
x = x[x > 0] # Gini is typically defined for non-negative values, ignore zeros
if len(x) == 0:
return 0.0
x = np.sort(x) # Values must be sorted
n = len(x)
# Numerically stable formula for Gini coefficient
# G = (2 * sum_{i=1 to n} i*x_i) / (n * sum_{i=1 to n} x_i) - (n+1)/n
# Or simplified form, which is often more robust:
# G = (np.sum((2 * np.arange(1, n + 1) - n - 1) * x)) / (n * np.sum(x))
# This formula assumes 1-based indexing for i, which np.arange(1, n+1) provides.
return (2 * np.sum(np.arange(1, n + 1) * x) - (n + 1) * np.sum(x)) / (n * np.sum(x))
def parted_run_gini(o: np.ndarray, tau: float, fps: float) -> float:
"""
Calculates the Gini coefficient of the durations of consecutive
"parted" lip segments (i.e., when opening > tau).
"""
mask = (o > tau) & ~np.isnan(o)
runs = [] # Store lengths of consecutive 'True' (parted) segments
i = 0
while i < len(mask):
if mask[i]:
j = i
while j < len(mask) and mask[j]:
j += 1
runs.append(j - i) # Length in frames
i = j
else:
i += 1
if not runs:
return 0.0
runs_sec = np.array(runs) / fps # Convert run lengths from frames to seconds
return gini_coefficient(runs_sec)
def burstiness_from_peaks(peaks: np.ndarray, fps: float) -> float:
"""
Calculates the Burstiness Index from the inter-peak intervals.
B = (sigma - mu) / (sigma + mu), where sigma is std dev and mu is mean of ISIs.
B ranges from -1 (very regular) to 1 (very bursty), 0 for Poisson.
"""
if len(peaks) < 2:
return 0.0
isi = np.diff(peaks) / fps # Inter-spike intervals in seconds
mu, sigma = np.mean(isi), np.std(isi)
if sigma + mu == 0: # Avoid division by zero, e.g., if all ISIs are zero (impossible for diff)
# or if ISI are constant and 0 (if ISIs are all 0, mu=0, sigma=0)
return 0.0
return (sigma - mu) / (sigma + mu)
def robust_z(x: float, ref_vals: np.ndarray, eps: float = 1e-4) -> float: # MODIFIED: Changed default eps from 1e-8 to 1e-4
"""
Computes a robust Z-score for a value `x` relative to a reference distribution `ref_vals`.
Uses median for central tendency and Median Absolute Deviation (MAD) for dispersion.
MAD is scaled by 1.4826 to be comparable to standard deviation for normal distributions.
A larger 'eps' (e.g., 1e-4) is used to prevent extremely large Z-scores when
the reference distribution's MAD is effectively zero.
"""
ref_vals = ref_vals[~np.isnan(ref_vals)]
if ref_vals.size < 2: # Need at least two data points for MAD to be meaningful
return 0.0 # Or consider returning NaN if data is insufficient for robust comparison
med = np.nanmedian(ref_vals)
mad = np.nanmedian(np.abs(ref_vals - med))
# Calculate denominator, ensuring it's not too small due to very low MAD
denominator = 1.4826 * (mad + eps)
# Handle the extremely rare case where denominator is still zero (e.g., ref_vals all NaN, or bad eps)
if denominator == 0:
# If x is also equal to the median (within epsilon tolerance), Z-score is 0.
# Otherwise, it's an extreme deviation, return a large (signed) value.
return 0.0 if np.isclose(x, med, atol=eps) else np.sign(x - med) * 1e6
return (x - med) / denominator
def pca_weights(feature_matrix: np.ndarray, feature_names: List[str]) -> Dict[str, float]:
"""
Calculates weights for features based on the absolute loadings of the first
principal component (PC1). Features contributing more to PC1 get higher weights.
Falls back to equal weights if PCA fails (e.g., insufficient data or no variance).
"""
try:
# Remove rows with any NaN values as PCA cannot handle them
X = feature_matrix.copy()
X = X[~np.isnan(X).any(axis=1)]
# Need at least 2 samples to compute variance, and more samples than features for stable PCA
if X.shape[0] < 2 or (X.shape[1] > X.shape[0] and X.shape[0] >=1): # Allow 1 sample but PCA won't give meaningful weights
raise ValueError("Not enough samples or too many features for robust PCA")
if X.shape[0] == 0:
raise ValueError("No valid samples after NaN removal for PCA")
# Center the data by subtracting the mean of each feature
Xc = X - X.mean(axis=0, keepdims=True)
# Perform Singular Value Decomposition (SVD) on the centered data
U, S, VT = np.linalg.svd(Xc, full_matrices=False)
# The first row of VT contains the loadings for the first principal component
pc1 = VT[0]
# Use the absolute values of the loadings as weights, as magnitude matters, not direction
w = np.abs(pc1)
# Normalize weights to sum to 1
if np.allclose(w.sum(), 0):
# This can happen if all features were constant or had zero variance across samples
raise ValueError("Zero sum of loadings from PCA, likely due to constant features.")
w = w / w.sum()
return {name: float(w[i]) for i, name in enumerate(feature_names)}
except Exception as e:
print(f"Warning: PCA weighting failed ({e}), falling back to equal weights.", file=sys.stderr)
eq = 1.0 / len(feature_names)
return {name: eq for name in feature_names}
def compute_slip_components(
openings_segment: np.ndarray,
fps: float,
tau: float,
min_twitch_height: float,
min_twitch_prominence: float,
smoothing_window_frames: int,
) -> Dict[str, float]:
"""
Calculates all individual feature components that make up the SLIP-Ω metric
for a given segment of mouth opening data.
"""
o_clean = openings_segment[~np.isnan(openings_segment)]
# Kinematics (Velocity, Acceleration, Jerk)
v = deriv(openings_segment, 1, fps)
a = deriv(openings_segment, 2, fps)
j = deriv(openings_segment, 3, fps) # Jerk is 3rd derivative
vrms = float(np.sqrt(np.mean(v ** 2))) if len(v) else 0.0
arms = float(np.sqrt(np.mean(a ** 2))) if len(a) else 0.0
jrms = float(np.sqrt(np.mean(j ** 2))) if len(j) else 0.0
# Events & Burstiness
peaks, _ = _find_twitch_peaks(
openings_segment, fps,
min_twitch_height, min_twitch_prominence, smoothing_window_frames
)
seg_duration_sec = len(o_clean) / fps if fps > 0 else 0
# Stabilize rate in very short segments to avoid inflated values
seg_rate_den = max(seg_duration_sec, 2.0) # enforce a 2s minimum window
twitch_rate = len(peaks) / seg_rate_den if seg_rate_den > 0 else 0.0
burstiness = burstiness_from_peaks(peaks, fps)
zcr_v = zero_crossing_rate(v) * fps # ZCR scaled by FPS to be a rate per second
# Persistence Morphology
p_idle = float((o_clean > tau).mean()) if len(o_clean) > 0 else 0.0
gini_runs = parted_run_gini(openings_segment, tau, fps)
# Closed-mouth micro motion: RMS velocity only when the mouth is closed
tau_micro = max(0.05, tau * 0.5)
micro_motion = 0.0
if len(openings_segment) >= 2:
raw = openings_segment
diffs = np.diff(raw)
valid_pairs = (~np.isnan(raw[1:])) & (~np.isnan(raw[:-1]))
if np.any(valid_pairs):
v_pairs = fps * diffs[valid_pairs]
o_prev = raw[:-1][valid_pairs]
o_next = raw[1:][valid_pairs]
closed = (o_prev < tau_micro) & (o_next < tau_micro)
if np.any(closed):
micro_motion = float(np.sqrt(np.mean((v_pairs[closed]) ** 2)))
# Temporal Complexity & Spectrum
sampen_val = sample_entropy(o_clean, m=2) # m=2 is a common parameter for Sample Entropy
hf = hf_energy(v, fps)
tv1_val = total_variation(openings_segment, fps)
return dict(
vrms=vrms, arms=arms, jrms=jrms,
twitch_rate=twitch_rate, burstiness=burstiness, zcr_v=zcr_v,
p_idle=p_idle, gini_runs=gini_runs,
sampen=sampen_val, hf_energy=hf, tv1=tv1_val,
micro_motion=micro_motion,
)
def slip_omega_from_features(
features: Dict[str, float],
ref_feature_matrix: np.ndarray, # rows=segments, cols=features in feature_names order
feature_names: List[str],
weights: Dict[str, float],
) -> Tuple[float, Dict[str, float]]:
"""
SLIP-Ω quality score (0–100, higher = better) computed from per-feature robust z-scores
against the distribution of all segments in the video (ref_feature_matrix). For each
feature f, we compute z_f = robust_z(f, ref_f) and map to a quality score q_f via a
logistic transform q_f = 1/(1 + exp(alpha_f * z_f)). We then take the weighted average.
Rationale: during silence, higher motion/complexity/mouth-open are worse. Thus, if a
feature value is above the dataset median (z_f > 0), its q_f < 0.5; below median → q_f > 0.5.
"""
# Per-feature slopes (higher → stronger penalties for being above-median)
feature_alpha = {
"vrms": 1.0,
"arms": 1.0,
"jrms": 2.0,
"micro_motion": 1.4,
"twitch_rate": 1.2,
"burstiness": 1.0,
"zcr_v": 0.8,
"p_idle": 1.0,
"gini_runs": 0.7,
"sampen": 0.8,
"hf_energy": 1.4,
"tv1": 0.6,
}
# Absolute quality mappings for all features (higher = better)
# Using sigmoid-based curves that are more forgiving than exponential
abs_quality = {
# Kinematics: more reasonable thresholds
"vrms": lambda x: float(1.0 / (1.0 + 3.0 * max(x, 0.0))),
"arms": lambda x: float(1.0 / (1.0 + 8.0 * max(x, 0.0))),
"jrms": lambda x: float(1.0 / (1.0 + 15.0 * max(x, 0.0))),
# Closed-mouth subtle motion: still emphasize but less harsh
"micro_motion": lambda x: float(1.0 / (1.0 + 10.0 * max(x, 0.0))),
# Events
"twitch_rate": lambda x: float(1.0 / (1.0 + 2.0 * max(x, 0.0))),
"burstiness": lambda x: float(1.0 / (1.0 + 1.5 * abs(x))),
# Jitter/complexity emphasis
"zcr_v": lambda x: float(1.0 / (1.0 + 3.0 * max(x, 0.0))),
# Morphology
"p_idle": lambda x: float(1.0 / (1.0 + 5.0 * np.clip(x, 0.0, 1.0))),
"gini_runs": lambda x: float(1.0 / (1.0 + 2.0 * max(x, 0.0))),
# Complexity & spectrum
"sampen": lambda x: float(1.0 / (1.0 + 1.5 * max(x, 0.0))),
"hf_energy": lambda x: float(1.0 / (1.0 + 20.0 * max(x, 0.0))),
"tv1": lambda x: float(1.0 / (1.0 + 5.0 * max(x, 0.0))),
}
# Determine if reference distribution is strong enough
def has_strong_ref(mat: np.ndarray) -> bool:
if mat is None or mat.size == 0 or mat.shape[0] < 3:
return False
# Require non-trivial variability for at least a few features
mad_vals = []
for i in range(mat.shape[1]):
col = mat[:, i]
med = np.nanmedian(col)
mad = np.nanmedian(np.abs(col - med))
mad_vals.append(mad)
mad_vals = np.array(mad_vals)
return np.count_nonzero(mad_vals > 1e-6) >= max(2, mat.shape[1] // 3)
# Prepare reference columns lookup
ref_cols: Dict[str, np.ndarray] = {}
if ref_feature_matrix is not None and ref_feature_matrix.size > 0:
for i, name in enumerate(feature_names):
try:
ref_cols[name] = ref_feature_matrix[:, i]
except Exception:
ref_cols[name] = np.array([features.get(name, np.nan)])
else:
for name in feature_names:
ref_cols[name] = np.array([features.get(name, np.nan)])
normalized_features: Dict[str, float] = {}
weighted_scores: List[float] = []
weight_sum = 0.0
nseg = int(ref_feature_matrix.shape[0]) if ref_feature_matrix is not None and ref_feature_matrix.size > 0 else 0
use_ref = has_strong_ref(ref_feature_matrix)
# Mix absolute with relative; small nseg → rely on absolute
gamma_abs = 1.0 if nseg <= 1 else (min(1.0, 2.0 / nseg)) # 1 seg→1.0, 2 seg→1.0, 3→0.666, 5→0.4, >=8→<=0.25
for name in feature_names:
raw = features.get(name, np.nan)
# Absolute quality
if np.isnan(raw):
q_abs = 0.5
else:
q_abs = abs_quality.get(name, lambda x: 1.0 / (1.0 + max(x, 0.0)))(raw)
# Relative quality
if not use_ref or np.isnan(raw):
q_rel = 0.5
else:
ref = ref_cols.get(name, np.array([raw]))
z = robust_z(raw, ref)
alpha = feature_alpha.get(name, 1.0)
s = -alpha * z # we want below-median (z<0) → higher quality
# numerically stable sigmoid(s)
if s >= 0:
q_rel = 1.0 / (1.0 + float(np.exp(-min(s, 60.0))))
else:
es = float(np.exp(max(s, -60.0)))
q_rel = es / (1.0 + es)
# Blend
q = gamma_abs * q_abs + (1.0 - gamma_abs) * q_rel
q = float(np.clip(q, 0.0, 1.0))
normalized_features[name] = q
w = float(weights.get(name, 0.0))
if w > 0:
weighted_scores.append(q * w)
weight_sum += w
final_q = (sum(weighted_scores) / weight_sum) if weight_sum > 0 else 0.5
# Tie-breaker: if too close to neutral, favor absolute-only aggregation
if abs(final_q - 0.5) < 1e-6:
abs_scores = []
for name in feature_names:
raw = features.get(name, np.nan)
qa = 0.5 if np.isnan(raw) else abs_quality.get(name, lambda x: 1.0 / (1.0 + max(x, 0.0)))(raw)
w = float(weights.get(name, 0.0))
if w > 0:
abs_scores.append(qa * w)
if weight_sum > 0 and len(abs_scores) > 0:
final_q = float(sum(abs_scores) / weight_sum)
# Silence-fraction multiplier: more strongly favor videos with little silence overall
silent_fraction = features.get("silent_fraction", None)
if silent_fraction is not None and np.isfinite(silent_fraction):
# Moderate boost when silence fraction → 0, minimal boost near 100% silence
speakiness = float(max(0.0, min(1.0, 1.0 - silent_fraction)))
mult = 1.0 + 1.0 * speakiness # slope reduced from 1.5 → 1.0
final_q = float(np.clip(final_q * min(mult, 1.5), 0.0, 1.0)) # cap reduced from 2.0 → 1.5
# Lighter speakiness baseline blend so ref3 stays closer to Type 3
final_q = float(np.clip(0.85 * final_q + 0.15 * speakiness, 0.0, 1.0))
return float(100.0 * final_q), normalized_features
# ------------------------------------------------------------
# Per-segment metrics
# ------------------------------------------------------------
def per_segment_metrics(
openings: np.ndarray,
fps: float,
segs: List[SilenceSegment],
tau: float,
lam: float,
mu: float,
min_twitch_height: float,
min_twitch_prominence: float,
smoothing_window_frames: int,
# For SLIP-Ω normalization/weights:
slip_feature_names: List[str],
slip_weights: Dict[str, float],
ref_feature_matrix: np.ndarray,
) -> List[SegmentMetrics]:
"""
Calculates detailed metrics for each silence segment, including SLM, SLE, and SLIP-Ω components.
"""
out: List[SegmentMetrics] = []
for k, seg in enumerate(segs):
# Create a mask for the current segment
mask = np.zeros_like(openings, dtype=bool)
mask[seg.start_frame : seg.end_frame + 1] = True
idx = np.where(mask)[0]
o_seg = openings[idx] # Mouth opening values for this segment
o_seg_clean = o_seg[~np.isnan(o_seg)] # Cleaned values for calculations
# SLM metrics (re-calculated on segment for per-segment report)
# Note: slm_delta and slms_beta are designed to take a mask, which internally handles NaNs.
sd = slm_delta(openings, fps, mask)
sb = slms_beta(openings, fps, mask, tau=tau, lam=lam, mu=mu)
# Calculate individual kinematic components for consistency in SegmentMetrics
vrms_seg = 0.0
arms_seg = 0.0
p_idle_seg = 0.0
if len(o_seg_clean) >= 2:
v_seg = fps * np.diff(o_seg_clean)
vrms_seg = math.sqrt(np.mean(v_seg ** 2))
if len(o_seg_clean) >= 3:
a_seg = (fps ** 2) * np.diff(o_seg_clean, n=2)
arms_seg = math.sqrt(np.mean(a_seg ** 2))
if len(o_seg_clean) > 0:
p_idle_seg = float((o_seg_clean > tau).mean())
# SLE metrics
num_twitches_seg = detect_twitches(
o_seg, fps,
min_twitch_height=min_twitch_height,
min_twitch_prominence=min_twitch_prominence,
smoothing_window_frames=smoothing_window_frames,
)
seg_duration_sec = len(o_seg_clean) / fps if fps > 0 else 0
sle_twitch_rate_seg = num_twitches_seg / seg_duration_sec if seg_duration_sec > 0 else 0.0
sle_parted_time_ratio_seg = calculate_parted_time_ratio(o_seg, tau)
# SLIP-Ω components & score
slip_feats = compute_slip_components(
o_seg, fps, tau,
min_twitch_height, min_twitch_prominence, smoothing_window_frames
)
# Calculate SLIP-Ω score for this segment using global reference matrix and weights
slip_score, _ = slip_omega_from_features(
slip_feats, ref_feature_matrix, slip_feature_names, slip_weights
)
out.append(
SegmentMetrics(
segment_index=k,
start_time=seg.start_time,
end_time=seg.end_time,
num_frames=seg.length_frames,
vrms=vrms_seg,
arms=arms_seg,
p_idle=p_idle_seg,
slm_delta=sd,
slms_beta=sb,
num_twitches=num_twitches_seg,
sle_twitch_rate=sle_twitch_rate_seg,
sle_parted_time_ratio=sle_parted_time_ratio_seg,
jrms=slip_feats["jrms"],
zcr_v=slip_feats["zcr_v"],
burstiness=slip_feats["burstiness"],
sampen=slip_feats["sampen"],
hf_energy=slip_feats["hf_energy"],
tv1=slip_feats["tv1"],
gini_runs=slip_feats["gini_runs"],
micro_motion=slip_feats["micro_motion"],
slip_omega=slip_score,
)
)
return out
# ------------------------------------------------------------
# Debug overlay video
# ------------------------------------------------------------
def save_debug_video(
video_path: str,
out_path: str,
silent_mask: np.ndarray,
openings: np.ndarray,
fps: float,
tau: float,
landmark_cache: LipLandmarkCache, # NEW: Pass the landmark cache
):
"""
Generates a debug video overlaying silent segments in red,
showing mouth opening values, and optionally landmarks during silent periods.
"""
cap = cv2.VideoCapture(video_path)
if not cap.isOpened():
raise RuntimeError(f"Cannot open video: {video_path}")
w = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH))
h = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))
fourcc = cv2.VideoWriter_fourcc(*"mp4v")
writer = cv2.VideoWriter(out_path, fourcc, fps, (w, h))
frame_idx = 0
font = cv2.FONT_HERSHEY_SIMPLEX
print(f" - Generating debug overlay video: {out_path}...")
total_frames_in_video = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
for _ in tqdm(range(total_frames_in_video), desc="Debug Video", unit="f"):
ret, frame = cap.read()
if not ret:
break
# Apply red overlay if frame is silent
is_silent = frame_idx < len(silent_mask) and silent_mask[frame_idx]
if is_silent:
overlay = frame.copy()
cv2.rectangle(overlay, (0, 0), (w - 1, h - 1), (0, 0, 255), -1) # Red color
alpha = 0.15 # Transparency
frame = cv2.addWeighted(overlay, alpha, frame, 1 - alpha, 0)
# NEW: Draw landmarks if silent
current_frame_lip_pts_normalized = landmark_cache.get(frame_idx)
if current_frame_lip_pts_normalized is not None:
# Reshape current_frame_lip_pts_normalized back to (N, 2) structure
stored_pts_xy = current_frame_lip_pts_normalized.reshape(-1, 2) # (4, 2)
# Map stored indices to polygon order
ordered_pts_for_poly_norm = np.array([
stored_pts_xy[0], # Upper
stored_pts_xy[3], # Right
stored_pts_xy[1], # Lower
stored_pts_xy[2], # Left
])
# Scale to pixel coordinates
ordered_pts_for_poly_scaled = (ordered_pts_for_poly_norm * np.array([w, h])).astype(int)
# Draw circles for the key lip points (all 4 of the original ones)
original_lm_indices_order = [UPPER_LIP_CENTER, LOWER_LIP_CENTER, LIP_LEFT_CORNER, LIP_RIGHT_CORNER]
for i, lm_idx in enumerate(original_lm_indices_order):
center = tuple( (stored_pts_xy[i] * np.array([w, h])).astype(int) )
cv2.circle(frame, center, 3, (0, 255, 0), -1) # Green circle, filled
# Draw lines to connect the lip region, forming a closed polygon
pts_cv = ordered_pts_for_poly_scaled.reshape((-1, 1, 2))
cv2.polylines(frame, [pts_cv], True, (255, 0, 0), 2) # Blue lines, closed, 2px thickness
# Display mouth opening value
if frame_idx < len(openings):
val = openings[frame_idx]
if not np.isnan(val):
txt = f"o={val:.3f}"
# Green if mouth is considered "closed" (below tau), Red if "parted" (above tau)
color = (0, 255, 0) if val <= tau else (0, 0, 255)
cv2.putText(frame, txt, (20, 40), font, 1.0, color, 2, cv2.LINE_AA)
else:
cv2.putText(frame, "o=NaN", (20, 40), font, 1.0, (100, 100, 100), 2, cv2.LINE_AA) # Grey for NaN
else:
cv2.putText(frame, "o=N/A", (20, 40), font, 1.0, (100, 100, 100), 2, cv2.LINE_AA) # Grey if out of bounds
writer.write(frame)
frame_idx += 1
cap.release()
writer.release()
print(" - Debug video generation complete.")
# ------------------------------------------------------------
# Video processing (global + per-segment)
# ------------------------------------------------------------
def process_video(
video_path: str,
rms_db_thresh: float = -45.0,
hop_length: int = 512,
tau: float = 0.20,
lam: float = 0.4,
mu: float = 0.3,
min_silence_frames: int = 2,
# Twitch detection parameters
min_twitch_height: float = 0.1,
min_twitch_prominence: float = 0.05,
smoothing_window_frames: int = 3,
# SLIP-Ω
weights_mode: str = "equal", # "equal" or "pca"
save_overlay: Optional[str] = None, # Path for silence debug overlay video
save_landmark_overlay: Optional[str] = None, # Path for MediaPipe landmark overlay video (all frames)
save_landmark_silence_overlay: Optional[str] = None, # NEW: Path for landmark overlay ONLY in silence
save_csv: Optional[str] = None, # Path for per-segment CSV
progress: bool = True,
# NEW: Model and Type names for summary CSV
model_name: Optional[str] = None,
type_name: Optional[str] = None,
) -> Dict[str, float]:
"""
Main function to process a single video, extracting all SLM, SLE, and SLIP-Ω metrics.
"""
print(f" - Detecting silence for {video_path}...")
# Compute silence mask, get FPS and total estimated frames from audio track
silent_mask, fps, total_frames = compute_silence_mask(
video_path,
rms_db_thresh=rms_db_thresh,
hop_length=hop_length,
min_silence_frames=min_silence_frames,
)
segs = mask_to_segments(silent_mask, fps)
print(f" - Found {len(segs)} silence segments.")
print(f" - Extracting lip landmarks...")
# Extract lip landmarks. Cache all frames if any type of overlay (full or silent) is requested,
# otherwise only cache silent frames to save memory.
cache_all_frames_for_overlay = save_landmark_overlay is not None or save_landmark_silence_overlay is not None
cache_target_mask = None if cache_all_frames_for_overlay else silent_mask
cache = extract_lip_landmarks(
video_path,
total_frames=total_frames,
target_mask=cache_target_mask,
progress=progress,
landmark_overlay_path=save_landmark_overlay, # This generates a separate video of all landmarks
)
o = opening_series(cache) # Get the scalar mouth opening time series
# Calculate Global SLM/SLE metrics (over ALL silent frames)
print(f" - Calculating global metrics...")
g_slm_delta = slm_delta(o, fps, silent_mask)
g_slms_beta = slms_beta(o, fps, silent_mask, tau=tau, lam=lam, mu=mu)
# Filter mouth openings to only include silent periods and non-NaN values for global SLE/SLIP
overall_silent_openings = o[np.where(silent_mask)[0]]
overall_silent_openings_clean = overall_silent_openings[~np.isnan(overall_silent_openings)]
total_silent_duration_sec = len(overall_silent_openings_clean) / fps if fps > 0 else 0
g_num_twitches = detect_twitches(
overall_silent_openings,
fps,
min_twitch_height=min_twitch_height,
min_twitch_prominence=min_twitch_prominence,
smoothing_window_frames=smoothing_window_frames,
)
g_sle_twitch_rate = g_num_twitches / total_silent_duration_sec if total_silent_duration_sec > 0 else 0.0
g_sle_parted_time_ratio = calculate_parted_time_ratio(overall_silent_openings, tau)
# Define the ordered list of SLIP-Ω feature names for consistency
slip_feature_names = [
"vrms", "arms", "jrms", "micro_motion",
"twitch_rate", "burstiness", "zcr_v",
"p_idle", "gini_runs",
"sampen", "hf_energy", "tv1"
]
# Compute per-segment features FIRST to build the reference matrix for robust z-scores.
# This matrix contains feature values from all segments, used for normalization.
seg_feature_rows = []
for seg in segs:
# Create mask for current segment
mask = np.zeros_like(o, dtype=bool)
mask[seg.start_frame : seg.end_frame + 1] = True
idx = np.where(mask)[0]
o_seg = o[idx] # Get segment-specific mouth opening data
# Compute all SLIP-Ω features for this segment
slip_feats = compute_slip_components(
o_seg, fps, tau,
min_twitch_height, min_twitch_prominence, smoothing_window_frames
)
row = [slip_feats[k] for k in slip_feature_names]
seg_feature_rows.append(row)
# Convert list of rows to a NumPy array for PCA and robust_z
# Handle case with no segments gracefully (e.g., video has no silence)
ref_feature_matrix = np.array(seg_feature_rows) if len(seg_feature_rows) > 0 else np.zeros((1, len(slip_feature_names)))
# Determine SLIP-Ω feature weights (equal, PCA-based, or expert-designed)
if weights_mode.lower() == "pca":
slip_weights = pca_weights(ref_feature_matrix, slip_feature_names)
elif weights_mode.lower() == "expert":
# Expert weights emphasize subtle motion during silence and closed-mouth movement
slip_weights = {
# Kinematics (dominant but balanced)
"vrms": 0.12,
"arms": 0.08,
"jrms": 0.20,
# Micro-motion and jitter/complexity carry strong signal
"micro_motion": 0.20,
"zcr_v": 0.12,
"hf_energy": 0.14,
"sampen": 0.05,
# Events
"twitch_rate": 0.05,
"burstiness": 0.00,
# Morphology (light touch)
"p_idle": 0.02,
"gini_runs": 0.02,
# Regularization
"tv1": 0.00,
}
else: # Default to equal weights
slip_weights = {k: 1.0 / len(slip_feature_names) for k in slip_feature_names}
# Calculate per-segment metrics including SLIP-Ω scores
print(f" - Calculating per-segment metrics...")
seg_metrics = per_segment_metrics(
o, fps, segs, tau=tau, lam=lam, mu=mu,
min_twitch_height=min_twitch_height,
min_twitch_prominence=min_twitch_prominence,
smoothing_window_frames=smoothing_window_frames,
slip_feature_names=slip_feature_names,
slip_weights=slip_weights,
ref_feature_matrix=ref_feature_matrix, # Pass the collected reference matrix
)
df = pd.DataFrame([asdict(sm) for sm in seg_metrics])
# Calculate Global SLIP-Ω
if total_silent_duration_sec == 0:
# No silence at all → metric not applicable; report 0 to keep schema stable
g_slip_omega = 0.0
elif len(overall_silent_openings_clean) == 0:
# Landmarks missing during silence; fall back to duration-weighted average of per-segment scores
if len(seg_metrics) > 0:
weights_frames = np.array([sm.num_frames for sm in seg_metrics], dtype=float)
scores = np.array([sm.slip_omega for sm in seg_metrics], dtype=float)
if np.sum(weights_frames) > 0:
g_slip_omega = float(np.sum(weights_frames * scores) / np.sum(weights_frames))
else:
g_slip_omega = float(np.mean(scores)) if len(scores) > 0 else 0.0
else:
g_slip_omega = 0.0
else:
# Compute global SLIP-Ω features based on all silent frames (works even for very short segments)
global_slip_feats = compute_slip_components(
overall_silent_openings, fps, tau,
min_twitch_height, min_twitch_prominence, smoothing_window_frames
)
# Attach global silence context for stability multiplier
global_slip_feats["silence_total_sec"] = total_silent_duration_sec
video_duration_sec = (total_frames / fps) if fps > 0 else 0.0
silent_fraction = (total_silent_duration_sec / video_duration_sec) if video_duration_sec > 0 else 0.0
global_slip_feats["silent_fraction"] = silent_fraction
# Use the redesigned SLIP-Ω calculation (higher = better quality)
g_slip_omega, g_normalized_features = slip_omega_from_features(
global_slip_feats, None, slip_feature_names, slip_weights
)
# Optional saving of results
if save_csv:
print(f" - Saving per-segment CSV to {save_csv}")
os.makedirs(os.path.dirname(save_csv), exist_ok=True)
df.to_csv(save_csv, index=False)
if save_overlay: # This is the original silence-only debug video
save_debug_video(
video_path, save_overlay, silent_mask=silent_mask, openings=o, fps=fps, tau=tau, landmark_cache=cache
)
if save_landmark_silence_overlay: # NEW: Debug video with landmarks only in silence
# The save_debug_video function already does exactly this when landmark_cache is provided.
# So we just call it with the new path for this specific debug video.
save_debug_video(
video_path, save_landmark_silence_overlay, silent_mask=silent_mask, openings=o, fps=fps, tau=tau, landmark_cache=cache
)
result_dict = {
"video": video_path,
"fps": fps,
"total_frames": int(total_frames),
"num_silence_segments": len(segs),
"total_silent_duration_sec": total_silent_duration_sec,
"slm_delta": g_slm_delta,
"slms_beta": g_slms_beta,
"sle_total_twitches": g_num_twitches,
"sle_twitch_rate": g_sle_twitch_rate,
"sle_parted_time_ratio": g_sle_parted_time_ratio,
"slip_omega": g_slip_omega,
"slip_weights_mode": weights_mode,
}
# Add model and type name if provided
if model_name is not None:
result_dict["model"] = model_name
if type_name is not None:
result_dict["type"] = type_name
return result_dict
def batch_process(
input_path: str,
recursive: bool,
# NEW: Model and Type names for summary CSV
model_name: Optional[str] = None,
type_name: Optional[str] = None,
**kwargs,
) -> pd.DataFrame:
"""
Processes multiple video files, either from a directory or a single file.
"""
exts = {".mp4", ".mov", ".mkv", ".avi", ".webm", ".m4v"}
files: List[str] = []
if os.path.isdir(input_path):
if recursive:
for root, _, fnames in os.walk(input_path):
for f in fnames:
if os.path.splitext(f.lower())[1] in exts:
files.append(os.path.join(root, f))
else:
for f in os.listdir(input_path):
if os.path.splitext(f.lower())[1] in exts:
files.append(os.path.join(input_path, f))
else:
# The issue with input_path being "inferred_vids/LatentSync/audio1//*/" happens here.
# The bash script is passing a glob pattern that doesn't resolve to a directory
# but to a literal string. This 'else' branch handles single files.
# For a directory processing, the 'if os.path.isdir(input_path)' must be true.
raise ValueError(f"Input file {input_path} is not a supported video format.")
if not files:
print(f"No supported video files found in {input_path}")
return pd.DataFrame()
rows = []
for i, v in enumerate(files):
print(f"\n=== Processing video {i+1}/{len(files)}: {v} ===")
base = os.path.splitext(os.path.basename(v))[0]
csv_path = None
save_csv_pattern = kwargs.get("save_csv_pattern")
if save_csv_pattern:
csv_path = save_csv_pattern.format(name=base)
os.makedirs(os.path.dirname(csv_path), exist_ok=True)
overlay_path = None
save_overlay_pattern = kwargs.get("save_overlay_pattern")
if save_overlay_pattern:
overlay_path = save_overlay_pattern.format(name=base)
os.makedirs(os.path.dirname(overlay_path), exist_ok=True)
landmark_overlay_path = None
save_landmark_overlay_pattern = kwargs.get("save_landmark_overlay_pattern")
if save_landmark_overlay_pattern:
landmark_overlay_path = save_landmark_overlay_pattern.format(name=base)
os.makedirs(os.path.dirname(landmark_overlay_path), exist_ok=True)
# NEW: Debug video for landmarks in silence
landmark_silence_overlay_path = None
save_landmark_silence_overlay_pattern = kwargs.get("save_landmark_silence_overlay_pattern")
if save_landmark_silence_overlay_pattern:
landmark_silence_overlay_path = save_landmark_silence_overlay_pattern.format(name=base)
os.makedirs(os.path.dirname(landmark_silence_overlay_path), exist_ok=True)
try:
row = process_video(
v,
save_csv=csv_path,
save_overlay=overlay_path,
save_landmark_overlay=landmark_overlay_path,
save_landmark_silence_overlay=landmark_silence_overlay_path, # Pass the new argument
progress=True,
rms_db_thresh=kwargs.get("rms_db_thresh", -45.0),
hop_length=kwargs.get("hop_length", 512),
tau=kwargs.get("tau", 0.2),
lam=kwargs.get("lam", 0.4),
mu=kwargs.get("mu", 0.3),
min_silence_frames=kwargs.get("min_silence_frames", 2),
# Twitch detection
min_twitch_height=kwargs.get("min_twitch_height", 0.1),
min_twitch_prominence=kwargs.get("min_twitch_prominence", 0.05),
smoothing_window_frames=kwargs.get("smoothing_window_frames", 3),
# SLIP-Ω
weights_mode=kwargs.get("weights_mode", "equal"),
# Pass model and type names
model_name=model_name,
type_name=type_name,
)
rows.append(row)
except Exception as e:
print(f"Error processing {v}: {e}", file=sys.stderr)
print("Skipping this video.", file=sys.stderr)
rows.append({"video": v, "status": "failed", "error": str(e), "model": model_name, "type": type_name}) # Ensure model/type for errors
df = pd.DataFrame(rows)
return df
# ------------------------------------------------------------
# CLI
# ------------------------------------------------------------
def build_arg_parser() -> argparse.ArgumentParser:
p = argparse.ArgumentParser(
description="Compute SLM, SLE, and SLIP-Ω metrics for lip-sync videos."
)
p.add_argument("input", help="Path to a video file OR a directory of videos.")
p.add_argument("-r", "--recursive", action="store_true",
help="Recurse into subdirectories when `input` is a directory.")
# Silence & audio parameters
p.add_argument("--rms-db-thresh", type=float, default=-45.0,
help="Silence detection threshold in dB relative to max RMS. (Default: -45.0)")
p.add_argument("--hop-length", type=int, default=512,
help="Librosa hop length for RMS analysis. (Default: 512)")
p.add_argument("--min-silence-frames", type=int, default=2,
help="Minimum run length (frames) to retain a silence island. (Default: 2)")
# SLMS parameters (also affects SLE and SLIP-Ω `p_idle` and `gini_runs`)
p.add_argument("--tau", type=float, default=0.20,
help="Idle-open threshold (mouth parted) for p_idle & SLE-PartedTimeRatio. (Default: 0.20)")
p.add_argument("--lam", type=float, default=0.4,
help="Weight on acceleration RMS in SLMS-beta. (Default: 0.4)")
p.add_argument("--mu", type=float, default=0.3,
help="Weight on idle proportion in SLMS-beta. (Default: 0.3)")
# Twitch detection parameters (affect SLE-TwitchRate and SLIP-Ω `twitch_rate`, `burstiness`)
p.add_argument("--min-twitch-height", type=float, default=0.1,
help="Minimum opening value for a twitch peak. (Default: 0.1)")
p.add_argument("--min-twitch-prominence", type=float, default=0.05,
help="Minimum prominence for a twitch peak. (Default: 0.05)")
p.add_argument("--smoothing-window-frames", type=int, default=3,
help="Rolling mean window size (in frames) for smoothing opening data before twitch detection. Should be odd. (Default: 3)")
# SLIP-Ω specific parameters
p.add_argument("--weights-mode", type=str, default="expert", choices=["equal", "pca", "expert"],
help="Method to set feature weights for SLIP-Ω composite score. 'equal' assigns uniform weights; 'pca' uses absolute loadings of PC1; 'expert' uses domain knowledge-based weights optimized for lip-sync quality. (Default: expert)")
# Output options
p.add_argument("--save-csv-pattern", type=str, default=None,
help="Pattern for saving per-segment metrics to CSV, e.g., 'out/{name}_segments.csv'. '{name}' is replaced by video basename.")
p.add_argument("--save-overlay-pattern", type=str, default=None,
help="Pattern for saving silence-debug overlay video (red overlay + opening value). e.g., 'out/{name}_overlay.mp4'. '{name}' is replaced by video basename.")
p.add_argument("--save-landmark-overlay-pattern", type=str, default=None,
help="Pattern for saving MediaPipe lip landmark overlay video (landmarks on ALL frames). e.g., 'out/{name}_landmarks.mp4'. '{name}' is replaced by video basename.")
p.add_argument("--save-landmark-silence-overlay-pattern", type=str, default=None, # NEW ARGUMENT
help="Pattern for saving a debug video with silence overlay AND lip landmarks ONLY drawn during silent segments. e.g., 'out/{name}_silent_landmarks.mp4'. '{name}' is replaced by video basename.")
p.add_argument("--summary-csv", type=str, default=None,
help="Write global metrics of all processed videos to this summary CSV file.")
p.add_argument("--vad-threshold", type=float, default=0.5,
help="Silero VAD speech probability threshold (0–1). Lower = more sensitive to speech.")
p.add_argument("--vad-sr", type=int, default=16000,
help="Resample audio to this rate for VAD (16k recommended).")
p.add_argument("--vad-onnx", action="store_true",
help="Use ONNX variant of Silero VAD instead of Torch.")
# NEW arguments for model and type name
p.add_argument("--model-name", type=str, default=None,
help="Name of the model (e.g., LatentSync, MuseTalk) for tracking in summary CSV.")
p.add_argument("--type-name", type=str, default=None,
help="Name of the video type (e.g., Type1, Type2) for tracking in summary CSV.")
return p
def main():
parser = build_arg_parser()
args = parser.parse_args()
# Ensure smoothing window is odd
if args.smoothing_window_frames % 2 == 0:
print("Warning: --smoothing-window-frames should be an odd number. Adjusting by +1.", file=sys.stderr)
args.smoothing_window_frames += 1
if args.smoothing_window_frames < 1: # Ensure it's at least 1
args.smoothing_window_frames = 1
# Call batch processing with all arguments
df = batch_process(
args.input,
recursive=args.recursive,
rms_db_thresh=args.rms_db_thresh,
hop_length=args.hop_length,
tau=args.tau,
lam=args.lam,
mu=args.mu,
min_silence_frames=args.min_silence_frames,
# Twitch detection
min_twitch_height=args.min_twitch_height,
min_twitch_prominence=args.min_twitch_prominence,
smoothing_window_frames=args.smoothing_window_frames,
# SLIP-Ω
weights_mode=args.weights_mode,
# I/O paths
save_csv_pattern=args.save_csv_pattern,
save_overlay_pattern=args.save_overlay_pattern,
save_landmark_overlay_pattern=args.save_landmark_overlay_pattern,
save_landmark_silence_overlay_pattern=args.save_landmark_silence_overlay_pattern, # NEW
# Pass model and type names
model_name=args.model_name,
type_name=args.type_name,
)
print("\n=== Processing Summary ===")
if not df.empty:
# Check if the summary CSV path was provided
if args.summary_csv:
summary_output_dir = os.path.dirname(args.summary_csv)
# Only create directory if the path contains a directory component
if summary_output_dir:
os.makedirs(summary_output_dir, exist_ok=True)
# Appends to the summary CSV without header
file_exists_before_write = os.path.exists(args.summary_csv)
df.to_csv(args.summary_csv, index=False, mode='a', header=not file_exists_before_write)
# Also print to console for immediate feedback for the current batch
print(df.to_string(index=False))
else:
print("No videos processed or no results generated for this batch.")
if __name__ == "__main__":
main()Editor is loading...
Leave a Comment