library.shoot_mask_cleaning

  1"""Shoot mask cleaning pipeline for plant root analysis.
  2
  3This module provides functions for cleaning and filtering shoot segmentation masks
  4from U-Net model predictions. The pipeline identifies exactly 5 shoot locations
  5using adaptive Y-zone detection and X-axis peak finding, handling edge cases like
  6failed germination and fallen shoots.
  7
  8Typical workflow:
  9    
 10    Step 1: Calculate global statistics across all masks (one-time setup)
 11    
 12        >>> from pathlib import Path
 13        >>> import cv2
 14        >>> from shoot_mask_cleaning import calculate_global_y_stats
 15        >>> 
 16        >>> mask_dir = Path('data/shoot_masks')
 17        >>> mask_paths = [str(f) for f in sorted(mask_dir.glob('*.png'))]
 18        >>> 
 19        >>> global_stats = calculate_global_y_stats(
 20        ...     mask_paths, 
 21        ...     kernel_size=5,
 22        ...     iterations=3,
 23        ...     y_min=200, 
 24        ...     y_max=750
 25        ... )
 26        >>> print(f"Global mean Y: {global_stats['global_mean']:.1f}")
 27        >>> print(f"Global std Y: {global_stats['global_std']:.1f}")
 28    
 29    Step 2: Process individual masks
 30    
 31        >>> from shoot_mask_cleaning import clean_shoot_mask_pipeline
 32        >>> 
 33        >>> for mask_path in mask_paths:
 34        ...     result = clean_shoot_mask_pipeline(mask_path, global_stats)
 35        ...     
 36        ...     # Save cleaned mask
 37        ...     output_path = f"cleaned/{result['filename']}"
 38        ...     cv2.imwrite(output_path, result['cleaned_mask'])
 39        ...     
 40        ...     # Check results
 41        ...     print(f"{result['filename']}: {result['method']}, "
 42        ...           f"{result['num_components']} components")
 43    
 44    Step 3: Debug problematic images
 45    
 46        >>> from shoot_mask_visualization import debug_peak_detection
 47        >>> 
 48        >>> # Visualize complete pipeline for troubleshooting
 49        >>> debug_peak_detection('data/shoot_masks/problem_image.png', global_stats)
 50
 51Recommended parameters (validated on 19 test images):
 52    
 53    Initial closing:
 54        - kernel_size: 7
 55        - iterations: 5
 56    
 57    Y-zone boundaries:
 58        - y_min: 200 (shoots never above this)
 59        - y_max: 750 (shoots never below this)
 60    
 61    Peak detection:
 62        - n_peaks: 5 (expected number of plants)
 63        - min_distance: 300 (minimum spacing between shoots)
 64        - x_min: 1000 (shoots never left of this)
 65        - initial_std: 2.0 (starting Y-zone width)
 66        - quality_check_threshold: 2.5 (when to validate peak quality)
 67        - min_peak_width: 20 (minimum width for valid peaks)
 68        - min_peak_height: 10 (minimum height for valid peaks)
 69        - min_area: 500 (for size-based fallback)
 70    
 71    Filtering:
 72        - band_width: 100 (X-distance around peaks)
 73
 74Dependencies:
 75    - numpy
 76    - opencv-python (cv2)
 77    - scipy (for peak detection)
 78
 79Author: Aaron Ciuffo
 80Date: December 2024
 81"""
 82
 83import numpy as np
 84import cv2
 85from pathlib import Path
 86from scipy.signal import find_peaks, peak_widths
 87
 88
 89def join_shoot_fragments(mask, kernel_size=7, iterations=5):
 90    """Join small shoot fragments using morphological closing.
 91    
 92    Applies morphological closing to connect nearby fragments and fill small gaps
 93    in the segmentation mask. This preprocessing step improves peak detection by
 94    creating more continuous shoot structures.
 95    
 96    Args:
 97        mask: Binary mask array with values 0 and 255.
 98        kernel_size: Size of square structuring element (default: 7).
 99        iterations: Number of closing iterations (default: 5).
100        
101    Returns:
102        Binary mask array with joined fragments (values 0 and 255).
103        
104    Example:
105        >>> mask = cv2.imread('shoot_mask.png', cv2.IMREAD_GRAYSCALE)
106        >>> joined = join_shoot_fragments(mask, kernel_size=5, iterations=3)
107        >>> cv2.imwrite('joined_mask.png', joined)
108    """
109    kernel = np.ones((kernel_size, kernel_size), np.uint8)
110    joined = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, kernel, iterations=iterations)
111    return joined
112
113
114def calculate_y_density(mask):
115    """Calculate pixel density per Y coordinate.
116    
117    Sums the number of foreground pixels along each row to create a 1D density
118    profile. Used for finding the vertical position of shoots.
119    
120    Args:
121        mask: Binary mask with values 0 and 255.
122        
123    Returns:
124        1D numpy array of pixel counts per row (length = mask height).
125        
126    Example:
127        >>> mask = cv2.imread('shoot_mask.png', cv2.IMREAD_GRAYSCALE)
128        >>> density = calculate_y_density(mask)
129        >>> print(f"Peak density at Y={np.argmax(density)}")
130    """
131    return np.sum(mask > 0, axis=1)
132
133
134def calculate_weighted_y_stats(density, y_min=200, y_max=750):
135    """Calculate weighted mean and standard deviation of Y positions within ROI.
136    
137    Computes the center of mass and spread of shoot pixels along the Y-axis,
138    restricted to the valid Y-range. Pixels outside the ROI are ignored to
139    prevent noise from biasing statistics.
140    
141    Args:
142        density: 1D array of pixel counts per row from calculate_y_density.
143        y_min: Minimum Y coordinate to consider (default: 200).
144        y_max: Maximum Y coordinate to consider (default: 750).
145        
146    Returns:
147        Dictionary with keys:
148            - 'mean': Weighted mean Y position (pixels).
149            - 'std': Weighted standard deviation (pixels).
150        
151    Example:
152        >>> density = calculate_y_density(mask)
153        >>> stats = calculate_weighted_y_stats(density, y_min=200, y_max=750)
154        >>> print(f"Center of mass: Y={stats['mean']:.1f} ± {stats['std']:.1f}")
155    """
156    # Clip density to ROI
157    roi_density = density.copy()
158    roi_density[:y_min] = 0
159    roi_density[y_max:] = 0
160    
161    y_positions = np.arange(len(roi_density))
162    total_pixels = np.sum(roi_density)
163    
164    if total_pixels == 0:
165        return {'mean': 0, 'std': 0}
166    
167    # Weighted mean (center of mass)
168    weighted_mean = np.sum(y_positions * roi_density) / total_pixels
169    
170    # Weighted standard deviation
171    weighted_var = np.sum(roi_density * (y_positions - weighted_mean)**2) / total_pixels
172    weighted_std = np.sqrt(weighted_var)
173    
174    return {
175        'mean': weighted_mean,
176        'std': weighted_std
177    }
178
179
180def calculate_global_y_stats(mask_paths, kernel_size=7, iterations=5, y_min=200, y_max=750):
181    """Calculate global weighted mean and standard deviation across all masks.
182    
183    Processes all masks to establish global shoot position statistics. These
184    statistics serve as priors for individual image processing, enabling adaptive
185    Y-zone detection that handles variation in shoot positions.
186    
187    Args:
188        mask_paths: List of paths (strings or Path objects) to mask files.
189        kernel_size: Kernel size for initial morphological closing (default: 7).
190        iterations: Number of closing iterations (default: 5).
191        y_min: Minimum Y coordinate to consider (default: 200).
192        y_max: Maximum Y coordinate to consider (default: 750).
193        
194    Returns:
195        Dictionary with keys:
196            - 'global_mean': Mean Y position across all images (pixels).
197            - 'global_std': Mean standard deviation across all images (pixels).
198            - 'all_stats': List of per-image statistics dictionaries.
199            - 'all_means': List of individual image means.
200            - 'all_stds': List of individual image standard deviations.
201            - 'y_min': Y minimum boundary used.
202            - 'y_max': Y maximum boundary used.
203        
204    Example:
205        >>> mask_files = [str(f) for f in Path('data/shoot').glob('*.png')]
206        >>> global_stats = calculate_global_y_stats(mask_files, y_min=200, y_max=750)
207        >>> print(f"Expected shoot zone: {global_stats['global_mean']:.0f} ± "
208        ...       f"{2 * global_stats['global_std']:.0f} pixels")
209    
210    Note:
211        This function should be called once per dataset to establish the global
212        statistics used for processing all individual images.
213    """
214    all_stats = []
215    
216    for path in mask_paths:
217        mask = cv2.imread(str(path), cv2.IMREAD_GRAYSCALE)
218        joined = join_shoot_fragments(mask, kernel_size, iterations)
219        density = calculate_y_density(joined)
220        stats = calculate_weighted_y_stats(density, y_min, y_max)
221        all_stats.append(stats)
222    
223    # Calculate mean of means and mean of stds
224    means = [s['mean'] for s in all_stats]
225    stds = [s['std'] for s in all_stats]
226    
227    global_mean = np.mean(means)
228    global_std = np.mean(stds)
229    
230    return {
231        'global_mean': global_mean,
232        'global_std': global_std,
233        'all_stats': all_stats,
234        'all_means': means,
235        'all_stds': stds,
236        'y_min': y_min,
237        'y_max': y_max
238    }
239
240
241def calculate_x_projection(mask, global_stats, std_multiplier=2.0):
242    """Calculate X-axis projection within the global Y-zone.
243    
244    Sums pixels along columns (Y-axis) within the vertical zone defined by global
245    statistics. Creates a 1D signal showing horizontal distribution of shoot mass.
246    
247    Args:
248        mask: Binary mask with values 0 and 255.
249        global_stats: Global statistics from calculate_global_y_stats.
250        std_multiplier: Standard deviation multiplier for Y-zone width (default: 2.0).
251        
252    Returns:
253        1D numpy array of pixel counts per column (length = mask width).
254        
255    Example:
256        >>> x_proj = calculate_x_projection(mask, global_stats, std_multiplier=2.0)
257        >>> import matplotlib.pyplot as plt
258        >>> plt.plot(x_proj)
259        >>> plt.xlabel('X coordinate')
260        >>> plt.ylabel('Pixel count')
261        >>> plt.show()
262    """
263    # Calculate Y bounds
264    global_upper = int(global_stats['global_mean'] - std_multiplier * global_stats['global_std'])
265    global_lower = int(global_stats['global_mean'] + std_multiplier * global_stats['global_std'])
266    
267    # Extract Y-zone region and sum along Y axis
268    zone_region = mask[global_upper:global_lower, :]
269    x_projection = np.sum(zone_region > 0, axis=0)
270    
271    return x_projection
272
273
274def find_shoot_peaks(x_projection, n_peaks=5, min_distance=300, x_min=1000):
275    """Find peak locations in X-projection using scipy peak detection.
276    
277    Identifies the N strongest peaks in the horizontal projection that satisfy
278    spacing and position constraints. Returns peaks sorted left to right.
279    
280    Args:
281        x_projection: 1D array of pixel counts from calculate_x_projection.
282        n_peaks: Number of peaks to find (default: 5).
283        min_distance: Minimum distance between peaks in pixels (default: 300).
284        x_min: Minimum X coordinate to consider (default: 1000).
285        
286    Returns:
287        Numpy array of peak X-coordinates, sorted left to right.
288        
289    Example:
290        >>> x_proj = calculate_x_projection(mask, global_stats)
291        >>> peaks = find_shoot_peaks(x_proj, n_peaks=5, min_distance=300)
292        >>> print(f"Found peaks at X positions: {peaks}")
293    """
294    # Mask out x < x_min
295    masked_projection = x_projection.copy()
296    masked_projection[:x_min] = 0
297    
298    # Find peaks with minimum distance constraint
299    peaks, properties = find_peaks(masked_projection, distance=min_distance)
300    
301    # Get peak heights and select top n_peaks
302    peak_heights = masked_projection[peaks]
303    top_indices = np.argsort(peak_heights)[-n_peaks:]
304    top_peaks = peaks[top_indices]
305    
306    # Sort by position (left to right)
307    top_peaks = np.sort(top_peaks)
308    
309    return top_peaks
310
311
312def validate_peak_quality(x_projection, peaks, min_peak_height=10, min_peak_width=20):
313    """Check if detected peaks are likely to be real shoots versus noise.
314    
315    Validates peak quality by checking both height (signal strength) and width
316    (spatial extent). Narrow spikes indicate noise, while wide peaks indicate
317    actual shoots. Requires at least 3 peaks to pass validation.
318    
319    Args:
320        x_projection: 1D array of pixel counts.
321        peaks: Array of peak X-coordinates from find_shoot_peaks.
322        min_peak_height: Minimum height for a valid peak (default: 10).
323        min_peak_width: Minimum width at half-height for valid peak (default: 20).
324        
325    Returns:
326        Boolean indicating if peaks are high quality (True) or likely noise (False).
327        
328    Example:
329        >>> peaks = find_shoot_peaks(x_proj, n_peaks=5)
330        >>> if validate_peak_quality(x_proj, peaks):
331        ...     print("High quality peaks detected")
332        ... else:
333        ...     print("Peaks may be noise, consider widening Y-zone")
334    """
335    peak_heights = x_projection[peaks]
336    
337    # Calculate peak widths at half prominence
338    widths, _, _, _ = peak_widths(x_projection, peaks, rel_height=0.5)
339    
340    # Check if at least 3 peaks are both tall enough AND wide enough
341    valid_peaks = np.sum((peak_heights >= min_peak_height) & (widths >= min_peak_width))
342    
343    return valid_peaks >= 3
344
345
346def find_shoot_peaks_with_size_fallback(mask, global_stats, n_peaks=5, min_distance=300,
347                                        x_min=1000, initial_std=2.0, max_std=3.5,
348                                        std_step=0.25, min_area=500, min_peak_width=20,
349                                        quality_check_threshold=2.5):
350    """Adaptively find shoot peaks with fallback to size-based detection.
351    
352    Three-stage detection strategy:
353    1. Normal cases (std 2.0-2.25): X-projection peaks without quality checks
354    2. Widened search (std 2.5-3.5): Progressive Y-zone widening with quality validation
355    3. Size fallback: Largest components by area when projection methods fail
356    
357    This handles both typical shoots and edge cases like failed germination or
358    shoots that have fallen outside the typical vertical zone.
359    
360    Args:
361        mask: Binary mask with values 0 and 255.
362        global_stats: Global statistics from calculate_global_y_stats.
363        n_peaks: Target number of peaks (default: 5).
364        min_distance: Minimum distance between peaks in pixels (default: 300).
365        x_min: Minimum X coordinate to consider (default: 1000).
366        initial_std: Starting std multiplier (default: 2.0).
367        max_std: Maximum std multiplier to try (default: 3.5).
368        std_step: Step size for widening (default: 0.25).
369        min_area: Minimum area for size fallback in pixels (default: 500).
370        min_peak_width: Minimum peak width for quality validation (default: 20).
371        quality_check_threshold: Std threshold to activate quality checks (default: 2.5).
372        
373    Returns:
374        Tuple of (peaks, x_projection, std_multiplier_used, method_used) where:
375            - peaks: Array of peak X-coordinates
376            - x_projection: The X-projection array used
377            - std_multiplier_used: The std multiplier that succeeded
378            - method_used: Either "projection" or "size_fallback"
379        
380    Example:
381        >>> peaks, x_proj, std, method = find_shoot_peaks_with_size_fallback(
382        ...     mask, global_stats, quality_check_threshold=2.5
383        ... )
384        >>> print(f"Method: {method}, Std: {std:.2f}")
385        >>> print(f"Peaks at: {peaks}")
386    """
387    std_multiplier = initial_std
388    
389    # Try standard X-projection approach with progressive widening
390    while std_multiplier <= max_std:
391        x_projection = calculate_x_projection(mask, global_stats, std_multiplier)
392        peaks = find_shoot_peaks(x_projection, n_peaks, min_distance, x_min)
393        
394        if len(peaks) >= n_peaks:
395            # Only check quality if we're in desperate territory (high std)
396            if std_multiplier >= quality_check_threshold:
397                if validate_peak_quality(x_projection, peaks, min_peak_height=10, 
398                                       min_peak_width=min_peak_width):
399                    print(f"Found {len(peaks)} quality peaks with std_multiplier={std_multiplier:.2f}")
400                    return peaks, x_projection, std_multiplier, "projection"
401                else:
402                    print(f"Found {len(peaks)} peaks but quality too low (std={std_multiplier:.2f})")
403            else:
404                # Low std values - just trust the peaks
405                print(f"Found {len(peaks)} peaks with std_multiplier={std_multiplier:.2f}")
406                return peaks, x_projection, std_multiplier, "projection"
407        
408        std_multiplier += std_step
409    
410    # Fallback: Use connected components filtered by size and X position
411    print(f"Projection method insufficient, trying size-based fallback...")
412    
413    num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(mask, connectivity=8)
414    
415    valid_components = []
416    for i in range(1, num_labels):
417        area = stats[i, cv2.CC_STAT_AREA]
418        x_left = stats[i, cv2.CC_STAT_LEFT]
419        width = stats[i, cv2.CC_STAT_WIDTH]
420        bbox_center = x_left + width // 2
421        
422        if area >= min_area and bbox_center >= x_min:
423            valid_components.append({
424                'label': i,
425                'bbox_center': bbox_center,
426                'area': area
427            })
428    
429    # Sort by area (largest first) and take top n_peaks
430    valid_components.sort(key=lambda x: x['area'], reverse=True)
431    selected = valid_components[:n_peaks]
432    
433    # Extract X positions and sort left to right
434    peaks = np.array([c['bbox_center'] for c in selected])
435    peaks = np.sort(peaks)
436    
437    # Create full projection for visualization
438    x_projection = np.sum(mask > 0, axis=0)
439    
440    print(f"Found {len(peaks)} peaks using size-based fallback (min_area={min_area})")
441    return peaks, x_projection, max_std, "size_fallback"
442
443
444def merge_with_narrow_bands(mask, peaks, global_stats, std_multiplier, band_width=100):
445    """Find components near peaks and keep entire components (not clipped to bands).
446    
447    Uses narrow X-bands around detected peaks to identify which components belong
448    to each shoot location. Keeps complete components rather than clipping them
449    to band boundaries, preserving shoot morphology.
450    
451    Args:
452        mask: Binary mask with values 0 and 255.
453        peaks: Array of peak X-coordinates from find_shoot_peaks_with_size_fallback.
454        global_stats: Global statistics from calculate_global_y_stats.
455        std_multiplier: Std multiplier used for peak detection.
456        band_width: Half-width of X-band for identifying components (default: 100).
457        
458    Returns:
459        Binary mask with complete components near peaks (values 0 and 255).
460        
461    Example:
462        >>> peaks, x_proj, std, method = find_shoot_peaks_with_size_fallback(
463        ...     mask, global_stats
464        ... )
465        >>> filtered = merge_with_narrow_bands(mask, peaks, global_stats, std, 
466        ...                                    band_width=100)
467    """
468    # Find all connected components in the full mask
469    num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(mask, connectivity=8)
470    
471    # Create empty output
472    output_mask = np.zeros_like(mask)
473    
474    # For each peak, find components whose center falls in the band
475    kept_labels = set()
476    
477    for i, peak_x in enumerate(peaks):
478        x_left = peak_x - band_width
479        x_right = peak_x + band_width
480        
481        # Check each component
482        for label in range(1, num_labels):
483            if label in kept_labels:
484                continue  # Already assigned to another peak
485            
486            # Get component's X-center
487            comp_x_left = stats[label, cv2.CC_STAT_LEFT]
488            comp_width = stats[label, cv2.CC_STAT_WIDTH]
489            comp_x_center = comp_x_left + comp_width // 2
490            
491            # If center is in this peak's band, keep the ENTIRE component
492            if x_left <= comp_x_center <= x_right:
493                output_mask[labels == label] = 255
494                kept_labels.add(label)
495                print(f"  Peak {i+1} at X={peak_x}: keeping component with center at X={comp_x_center}")
496    
497    return output_mask
498
499
500def filter_one_per_peak(mask, peaks, global_stats, std_multiplier, band_width=100,
501                       require_y_zone=True):
502    """Keep the largest component near each peak, optionally validating Y-zone.
503    
504    For each detected peak, identifies all components within the X-band and keeps
505    only the largest one. Optionally filters out components that don't touch the
506    Y-zone (disabled for size-fallback method to handle fallen shoots).
507    
508    Args:
509        mask: Binary mask with values 0 and 255.
510        peaks: Array of peak X-coordinates.
511        global_stats: Global statistics from calculate_global_y_stats.
512        std_multiplier: Std multiplier used for peak detection.
513        band_width: Half-width for assigning components to peaks (default: 100).
514        require_y_zone: If True, only keep components touching Y-zone (default: True).
515        
516    Returns:
517        Binary mask with exactly one component per peak (values 0 and 255).
518        
519    Example:
520        >>> # For projection method (require Y-zone)
521        >>> final_mask = filter_one_per_peak(mask, peaks, global_stats, std_used,
522        ...                                  band_width=100, require_y_zone=True)
523        >>> 
524        >>> # For size fallback (allow fallen shoots outside Y-zone)
525        >>> final_mask = filter_one_per_peak(mask, peaks, global_stats, std_used,
526        ...                                  band_width=100, require_y_zone=False)
527    
528    Note:
529        Set require_y_zone=False when using size-based fallback method to preserve
530        large shoots that have fallen outside the typical vertical zone.
531    """
532    # Calculate Y-zone bounds
533    global_upper = int(global_stats['global_mean'] - std_multiplier * global_stats['global_std'])
534    global_lower = int(global_stats['global_mean'] + std_multiplier * global_stats['global_std'])
535    
536    # Find all components
537    num_labels, labels, stats, _ = cv2.connectedComponentsWithStats(mask, connectivity=8)
538    
539    output_mask = np.zeros_like(mask)
540    
541    for peak_x in peaks:
542        x_left = peak_x - band_width
543        x_right = peak_x + band_width
544        
545        # Find all components in this peak's band
546        candidates = []
547        for i in range(1, num_labels):
548            x_center = stats[i, cv2.CC_STAT_LEFT] + stats[i, cv2.CC_STAT_WIDTH] // 2
549            
550            if x_left <= x_center <= x_right:
551                area = stats[i, cv2.CC_STAT_AREA]
552                
553                if require_y_zone:
554                    # Check if component touches Y-zone
555                    component_mask = (labels == i).astype(np.uint8)
556                    y_zone_pixels = np.sum(component_mask[global_upper:global_lower, :])
557                    
558                    if y_zone_pixels > 0:
559                        candidates.append((i, area))
560                else:
561                    # No Y-zone requirement - keep all
562                    candidates.append((i, area))
563        
564        # Keep the largest one from this peak
565        if candidates:
566            best_label = max(candidates, key=lambda x: x[1])[0]
567            output_mask[labels == best_label] = 255
568            print(f"  Peak at X={peak_x}: keeping component {best_label} "
569                  f"(area={stats[best_label, cv2.CC_STAT_AREA]})")
570        else:
571            print(f"  Peak at X={peak_x}: WARNING - no valid components found!")
572    
573    return output_mask
574
575
576def clean_shoot_mask_pipeline(mask_path, global_stats, closing_kernel_size=7,
577                              closing_iterations=5, quality_check_threshold=2.5,
578                              band_width=100):
579    """Complete pipeline to clean a shoot mask from raw predictions to final output.
580    
581    Runs the full processing workflow:
582    1. Load mask and apply initial morphological closing
583    2. Detect 5 shoot locations using adaptive peak finding
584    3. Filter to components near detected peaks
585    4. Keep exactly one component per peak (largest in each band)
586    
587    Args:
588        mask_path: Path to shoot mask file (string or Path object).
589        global_stats: Global statistics from calculate_global_y_stats.
590        closing_kernel_size: Initial closing kernel size (default: 7).
591        closing_iterations: Initial closing iterations (default: 5).
592        quality_check_threshold: Std threshold for quality checks (default: 2.5).
593        band_width: X-band width around peaks in pixels (default: 100).
594        
595    Returns:
596        Dictionary with keys:
597            - 'cleaned_mask': Final cleaned binary mask (0 and 255)
598            - 'peaks': Detected peak X-coordinates
599            - 'std_used': Std multiplier that succeeded
600            - 'method': Detection method ("projection" or "size_fallback")
601            - 'num_components': Number of components in final mask
602            - 'filename': Input filename
603        
604    Example:
605        >>> from pathlib import Path
606        >>> 
607        >>> # Process all masks
608        >>> mask_dir = Path('data/shoot_masks')
609        >>> for mask_file in mask_dir.glob('*.png'):
610        ...     result = clean_shoot_mask_pipeline(str(mask_file), global_stats)
611        ...     
612        ...     # Save cleaned mask
613        ...     output_path = f"cleaned/{result['filename']}"
614        ...     cv2.imwrite(output_path, result['cleaned_mask'])
615        ...     
616        ...     # Log results
617        ...     if result['num_components'] != 5:
618        ...         print(f"WARNING: {result['filename']} has "
619        ...               f"{result['num_components']} components")
620    
621    Note:
622        Parameters are optimized for typical plant imaging conditions:
623        - Shoots at X > 1000 pixels
624        - Shoots between Y = 200-750 pixels (adaptive per image)
625        - ~400-500 pixel spacing between plants
626    """
627    print(f"Processing: {Path(mask_path).name}")
628    
629    # Step 1: Load and initial joining
630    mask = cv2.imread(str(mask_path), cv2.IMREAD_GRAYSCALE)
631    joined_mask = join_shoot_fragments(mask, closing_kernel_size, closing_iterations)
632    
633    # Step 2: Find shoot locations
634    peaks, x_projection, std_used, method = find_shoot_peaks_with_size_fallback(
635        joined_mask, global_stats,
636        quality_check_threshold=quality_check_threshold
637    )
638    print(f"  Found {len(peaks)} peaks using {method} (std={std_used:.2f})")
639    
640    # Step 3: Filter to narrow bands around peaks
641    narrow_result = merge_with_narrow_bands(joined_mask, peaks, global_stats, 
642                                           std_used, band_width=band_width)
643    
644    # Step 4: Keep one component per peak
645    require_y_zone = (method != "size_fallback")
646    final_mask = filter_one_per_peak(narrow_result, peaks, global_stats, std_used,
647                                    band_width=band_width, require_y_zone=require_y_zone)
648    
649    # Count final components
650    num_labels, _, _, _ = cv2.connectedComponentsWithStats(final_mask, connectivity=8)
651    num_components = num_labels - 1
652    
653    return {
654        'cleaned_mask': final_mask,
655        'peaks': peaks,
656        'std_used': std_used,
657        'method': method,
658        'num_components': num_components,
659        'filename': Path(mask_path).name
660    }
661
662
663# Recommended parameters as a reference
664RECOMMENDED_PARAMS = {
665    'initial_closing': {
666        'kernel_size': 7,
667        'iterations': 5
668    },
669    'y_zone_boundaries': {
670        'y_min': 200,  # Shoots never above this
671        'y_max': 750   # Shoots never below this
672    },
673    'peak_detection': {
674        'n_peaks': 5,
675        'min_distance': 300,          # Minimum spacing between shoots
676        'x_min': 1000,                # Shoots never left of this
677        'initial_std': 2.0,           # Starting Y-zone width
678        'max_std': 3.5,               # Maximum Y-zone width
679        'quality_check_threshold': 2.5,  # When to validate peaks
680        'min_peak_width': 20,         # Minimum width for valid peaks
681        'min_peak_height': 10,        # Minimum height for valid peaks
682        'min_area': 500               # For size-based fallback
683    },
684    'filtering': {
685        'band_width': 100  # X-distance around peaks
686    }
687}