Image_Analyzer.py

import cv2
import numpy as np
from skimage import color
import math
import os
import csv
from skimage.feature import peak_local_max
from PIL import Image, ImageChops
import matplotlib.pyplot as plt

"""
THIS SCRIPT ANALYZES AN IMAGE OF PIG SKIN TO FIND THE MICRONEEDLE INSERTION SITE.
IT EXPECTS THAT THE MICRONEEDLE INSERTION SITE IS EITHER IN
1) A SINGLE RECTANGLE .... OR 
2) TWO RECTANGLES JOINED VIA A COMMON MIDDLE LINE.
FIRST IT FINDS THIS RECTANGLE, AND THEN WITHIN THIS RECTANGLE
IT FINDS THE MICRONEEDLE INSERTION SITE
"""

# --- CONFIGURATION ---
TARGET_AREA_MULTIPLIER = 3  # ROI2 target area = 3x ROI1 area
MIN_BLOB_AREA = 50          # Minimum pixel area for insertion site
MAX_BLOB_AREA = 15000       # Maximum pixel area for insertion site
SAFE_ZONE_MARGIN = 0.20     # Skips the outer 20% of each sub-rectangle to avoid marker lines


def extract_cmyk_magenta_component(image_path, opacity=100.0):
    """Replicates GIMP's 'Extract Component -> CMYK Magenta' with Opacity blending."""
    
    # 1. Handle file path vs array inputs correctly
    if isinstance(image_path, str):
        # Read with OpenCV (BGR format)
        bgr_img = cv2.imread(image_path)
        if bgr_img is None:
            raise FileNotFoundError(f"Could not load image: {image_path}")
    else:
        bgr_img = image_path

    # Convert OpenCV BGR to RGB array format
    rgb_array = cv2.cvtColor(bgr_img, cv2.COLOR_BGR2RGB)

    # 2. Extract the Magenta channel using standard CMYK conversion formulas
    # Scale array to float32 [0, 1.0] for math precision
    rgb_float = rgb_array.astype(np.float32) / 255.0
    r, g, b = rgb_float[..., 0], rgb_float[..., 1], rgb_float[..., 2]

    # Calculate Key (Black) channel
    k = 1.0 - np.maximum(np.maximum(r, g), b)

    # Calculate Magenta channel
    m = np.zeros_like(g)
    mask = (1.0 - k) > 0
    m[mask] = (1.0 - g[mask] - k[mask]) / (1.0 - k[mask])

    # Convert extracted data to a 0-255 uint8 grayscale format
    magenta_channel_data = np.uint8(m * 255)
    
    # 3. Apply the GIMP "Blending Options: Opacity" slider
    if opacity == 100.0:
        # Return directly as a single-channel grayscale NumPy array
        return magenta_channel_data 
    else:
        # Blend using Pillow, then convert the final result back to a NumPy array
        pil_original = Image.fromarray(rgb_array)
        pil_magenta_rgb = Image.fromarray(magenta_channel_data, mode="L").convert("RGB")
        
        final_output = ImageChops.blend(pil_original, pil_magenta_rgb, opacity / 100.0)
        return np.array(final_output)


def gimp_high_pass(image_input, std_dev=69, contrast=2):
   # std_dev 607.9, contrast=3.787 
   # This setting was found by experimenting with few images
   # in GIMP
    """Replicates GIMP's High Pass filter.

    :param image_input: PIL Image or NumPy array (grayscale/L channel)
    :param std_dev: GIMP's 'Std. Dev.' slider (corresponds to Gaussian blur
    radius)
    :param contrast: GIMP's 'Contrast' slider (scales edge intensity)
    """
    # 1. Ensure input is a floating-point NumPy array (0.0 to 255.0)
    img_array = np.asarray(image_input, dtype=np.float32)

    # 2. Calculate the kernel size for Gaussian Blur based on Std. Dev.
    # OpenCV requires an odd integer for ksize. GIMP uses a very wide radius.
    k_size = int(round(std_dev * 2)) * 2 + 1
    # Avoid zero or excessively massive kernels that crash memory
    k_size = max(3, min(k_size, 2001))

    # 3. Create the low-frequency background via Gaussian Blur
    low_freq = cv2.GaussianBlur(img_array, (0,0), sigmaX=std_dev)

    # 4. High Pass Formula: Original - Blurred + Neutral Gray (127.5)
    high_pass = img_array - low_freq + 127.5

    # 5. Apply the Contrast multiplier centered around neutral gray
    if contrast != 1.0:
        high_pass = (high_pass - 127.5) * contrast + 127.5

    # 6. Clip values to valid 0-255 bounds and convert back to 8-bit integers
    high_pass_clipped = np.clip(high_pass, 0, 255).astype(np.uint8)

    # 7. Return as a Pillow Image object
    return high_pass_clipped


def load_side_mapping(csv_path):
    """
    An image has insertion sites for cMN and bMN
    cMN == coated MN == peanut_protein_coated_MN
    bMN == blank MN == excipient coated MN
    This loads the CSV file that maps filenames to cMN or bMN.
    """
    mapping = {}
    if not os.path.exists(csv_path):
        print(f"Warning: Mapping CSV {csv_path} not found. Cannot assign cMN/bMN properly.")
        return mapping
        
    with open(csv_path, mode='r') as f:
        reader = csv.reader(f)
        next(reader, None)  # Skip header
        for row in reader:
            if len(row) >= 2:
                mapping[row[0].strip()] = row[1].strip()
    return mapping

def calculate_erythema_index(bgr_img):
    img_float = bgr_img.astype(float)
    b, g, r = cv2.split(img_float)
    ei = np.log10(r + 1.0) - np.log10(g + 1.0)
    return ei

def get_a_star_channel(bgr_img):
    """
    Represents the green-to-red color axis within the CIELAB ($L^*a^*b^*$) color space.
    lab_img[:, :, 1]) isolates pure redness while stripping away all brightness and shadow data.
    In dermatological research, tracking the a* channel is desirable
    over tracking the raw "Red" channel from a standard RGB camera 
    for two major reasons:
        1. It is "Immune" to Light and Shadows
        2. It Directly Measures Erythema (Blood Flow)
        since a* channel is highly sensitive to
        subtle shifts in blood perfusion
    """
    rgb_img = cv2.cvtColor(bgr_img, cv2.COLOR_BGR2RGB)
    lab_img = color.rgb2lab(rgb_img)
    return lab_img[:, :, 1]

def get_magenta(bgr_img):
    """
    Extracts CMYK Magenta channel for metric calculations.
    
    We experimented with images in GIMP by changing filters or extracting color components.
    We found that magenta channel from CMYK space was especially
    good at decreasing background and showing the microneedle
    insertion site as a fluorescent area.
    Therefore, we used this channel for image analysis.
    """
    bgr = bgr_img.astype(np.float32) / 255.0
    b, g, r = cv2.split(bgr)
    k = 1.0 - np.maximum(np.maximum(r, g), b)
    m = (1.0 - g - k) / (1.0 - k + 1e-7)
    return (m * 255.0).astype(np.uint8)

def find_master_bounding_box(bgr_img):
    """
    Finds the overall rectangular boundary containing the marker drawing.
    Uses a multi-stage thresholding and morphology pipeline to isolate the hand-drawn grid line.
    """
    # Convert the input BGR color image to a 1-channel grayscale canvas for thresholding operations
    gray = cv2.cvtColor(bgr_img, cv2.COLOR_BGR2GRAY)
    # Extract high-level shape parameters (height and width) of the current image canvas
    img_h, img_w = bgr_img.shape[:2]

    # 1. CREATE THE STENCIL
    # Apply a heavy median blur with a large 25x25 kernel window. This smooths out high-frequency 
    # structural noise like skin textures, body hair, wrinkles, and micro-needling points, 
    # leaving behind only the broad, thick pen strokes of the marker boundaries.
    blurred = cv2.medianBlur(gray, 25)
    
    # Use local adaptive thresholding to compute localized thresholds across varying illumination zones.
    # The cv2.THRESH_BINARY_INV flag flips the array so that dark ink strokes turn into bright white 
    # foreground tracks (255) against an absolute black baseline (0).
    pristine_mask = cv2.adaptiveThreshold(blurred, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, 
                                          cv2.THRESH_BINARY_INV, 81, 5)

    # 2. CREATE THE HIDDEN TRACKER
    # Construct a solid 25x25 rectangular structuring element kernel block for morphological bridging
    kernel_search = cv2.getStructuringElement(cv2.MORPH_RECT, (25, 25))
    
    # Run a Morphological Closing step (dilation followed by erosion). This forces nearby disjointed line 
    # fragments, faint pen strokes, or small structural breaks in the hand-drawn grid box to fuse together, 
    # locking the outer bounding frame into a single continuous, solid object channel.
    search_mask = cv2.morphologyEx(pristine_mask, cv2.MORPH_CLOSE, kernel_search)

    # 3. EDGE REJECTION + LARGEST AREA
    contours, _ = cv2.findContours(search_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    best_box = None
    max_area = 0
    margin = 15 

    for cnt in contours:
        x, y, w, h = cv2.boundingRect(cnt)
        area = w * h
        
        # RULE 1: The shape CANNOT touch the extreme edges of the photo.
        if x > margin and y > margin and (x + w) < (img_w - margin) and (y + h) < (img_h - margin):
            # RULE 2: Out of the valid, strictly internal shapes, pick the largest one.
            if area > max_area:
                max_area = area
                best_box = (x, y, w, h)

    return best_box

def get_safe_zone(rect, margin_pct):
    """Shrinks a rectangle by a percentage margin to avoid edge artifacts."""
    x, y, w, h = rect
    sx = int(x + w * margin_pct)
    sy = int(y + h * margin_pct)
    sw = int(w * (1 - 2 * margin_pct))
    sh = int(h * (1 - 2 * margin_pct))
    return (sx, sy, sw, sh)



def show_image(img, title="Default"):
    """
    Helper function used while debugging and making the script.
    It helped visualize various intermediate images.
    """
    fig, ax = plt.subplots(figsize=(8, 6))
    plt.title(title, fontsize=20) 
    ax.imshow(img, cmap='gray')
    
    
def find_best_site_in_zone(img, zone_rect, label, debug_dir):
    """Finds microneedle insertion site and saves
    intermediate debug frames.
    """
    
    zx, zy, zw, zh = zone_rect
    
    crop_pct = 0.02 # crops a region inside the sharpie rectangle so that edge artifacts are removed
    cx, cy = int(zx + (zw * crop_pct)), int(zy + (zh * crop_pct))
    cw, ch = int(zw * (1 - 2 * crop_pct)), int(zh * (1 - 2 * crop_pct))

    roi_bgr = img[cy:cy+ch, cx:cx+cw]
    
    
    # --- REDNESS METHOD  ---
    b, g, r = cv2.split(roi_bgr.astype(np.float32))
    redness = np.clip(r - g, 0, 255).astype(np.uint8)
    # high pass filter accentuated the insertion site helping in its identification
    redness_high_pass = gimp_high_pass(redness)
    
    # --- MAGENTA HIGH-PASS METHOD ---
    magenta_img = extract_cmyk_magenta_component(roi_bgr)
    # show_image(magenta_img, "Magenta")

    # redness = gimp_high_pass(magenta_img)
    # show_image(redness, "magenta High Pass")

    # We found that both red and magenta channels were useful
    # in finding the microneedle insertion sites.
    # In some images red channel was better and in some magenta was better
    # So we created a fused channel to help identify sites in ALL images
    fused_redness = cv2.addWeighted(redness_high_pass, 0.75, magenta_img, 0.25, 0)
    
    # --- HAIR REMOVAL (Median Blur) ---
    # Erases thin hairs by replacing them with the median color of the surrounding skin
    median_k = max(5, int(cw * 0.02))
    if median_k % 2 == 0: median_k += 1
    redness = cv2.medianBlur(fused_redness, median_k)
    # show_image(redness, "Redness_Image")
    
    
    # STEP 1: CRISP MASK
    p85 = np.percentile(redness, 85)
    _, crisp_mask = cv2.threshold(redness, p85, 255, cv2.THRESH_BINARY)
    
    
    # --- PIXEL CLIPPER ---
    # Snips off any thin squiggly lines (surviving hairs) protruding from the main blob
    clipper_size = max(3, int(cw * 0.015))
    clipper_kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (clipper_size, clipper_size))
    crisp_mask = cv2.morphologyEx(crisp_mask, cv2.MORPH_OPEN, clipper_kernel)
    
    
    cv2.imwrite(os.path.join(debug_dir, f"debug_1_CRISP_MASK_{label}.jpg"), crisp_mask)
    
    # STEP 2: HEATMAP CORE & SPATIAL WEIGHTING
    kernel_size = int(cw // 2) 
    if kernel_size % 2 == 0: kernel_size += 1
    
    # Convert heatmap to float for accurate multiplication
    heatmap = cv2.GaussianBlur(crisp_mask, (kernel_size, kernel_size), 0).astype(np.float32)  
   
    # --- Center-Bias Gaussian Mask ---
    # We found that some images had too much noise on skin
    # Since the microneedle insertion site was largely positioned in the
    # middle of the hand-drawn rectangle, we created a 'center-bias'
    # This creates a mathematical "hill" that fades out at the edges
    yy, xx = np.mgrid[0:ch, 0:cw]
    center_x, center_y = cw // 2, ch // 2
    # The 0.4 multiplier sets how wide the "hill" is. 
    sigma_x, sigma_y = cw * 0.4, ch * 0.4 
    spatial_weight = np.exp(-(((xx - center_x)**2) / (2 * sigma_x**2) + ((yy - center_y)**2) / (2 * sigma_y**2)))
    
    # Multiply the raw redness heatmap by the center-bias mask
    weighted_heatmap = heatmap * spatial_weight

    # Find the hot zone using the weighted map
    max_heat = np.max(weighted_heatmap)
    _, hot_zone = cv2.threshold(weighted_heatmap, max_heat * 0.75, 255, cv2.THRESH_BINARY)
    hot_zone = hot_zone.astype(np.uint8) # Convert back to 8-bit image
    cv2.imwrite(os.path.join(debug_dir, f"debug_2_HOTZONE_{label}.jpg"), hot_zone)
    
    # STEP 3: RECONSTRUCTION
    bubble_radius = max(int(cw * 0.175), 1)
    bubble_kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (bubble_radius, bubble_radius))
    protective_bubble = cv2.dilate(hot_zone, bubble_kernel, iterations=1)
    cv2.imwrite(os.path.join(debug_dir, f"debug_3_BUBBLE_{label}.jpg"), protective_bubble)
    
    final_isolated_mask = cv2.bitwise_and(crisp_mask, protective_bubble)
    
    # This melts the individual microneedle puncture dots into one solid contiguous blob
    bridge_size = max(5, int(cw * 0.06)) # Sized to bridge gaps between microneedle holes
    bridge_kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (bridge_size, bridge_size))
    final_isolated_mask = cv2.morphologyEx(final_isolated_mask, cv2.MORPH_CLOSE, bridge_kernel)

    cv2.imwrite(os.path.join(debug_dir, f"debug_4_FINAL_ISOLATED_{label}.jpg"), final_isolated_mask)
    
    # STEP 4: DRAW BOUNDARY
    contours, _ = cv2.findContours(final_isolated_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    
    if contours:
        best_cnt = max(contours, key=cv2.contourArea)
        
        # # CIRCUMSCRIBED CIRCLE
        # Old version
        # Sometimes this gave a large bounding box around the
        # microneedle insertion site.
        # So we moved to the 'equivalent area' technique
        # xcc, ycc, wcc, hcc = cv2.boundingRect(best_cnt)
        
        
        # Equivalent Area Bounding Box
        # Area of rect = Area of contour marking microneedle insertion site.
        # This contour can be irregular rectangularish in shape
        
        # 1. Get the standard circumscribed outer box for the aspect ratio
        x_old, y_old, w_old, h_old = cv2.boundingRect(best_cnt)
        aspect_ratio = w_old / float(h_old) if h_old != 0 else 1.0
        
        # 2. Get the true biological area of the red/magenta blob
        actual_area = cv2.contourArea(best_cnt)
        
        # 3. Calculate the center of mass (Centroid) of the blob
        M = cv2.moments(best_cnt)
        if M['m00'] != 0:
            cx_blob = int(M['m10'] / M['m00'])
            cy_blob = int(M['m01'] / M['m00'])
        else:
            cx_blob, cy_blob = x_old + w_old // 2, y_old + h_old // 2
            
        # 4. Calculate the mathematically equivalent "mean" dimensions
        w = int(np.sqrt(actual_area * aspect_ratio))
        h = int(np.sqrt(actual_area / aspect_ratio))
        
        # 5. Define the new top-left corner so the box is perfectly centered
        x = cx_blob - (w // 2)
        y = cy_blob - (h // 2)
        
        debug_final = roi_bgr.copy()
        cv2.drawContours(debug_final, [best_cnt], -1, (255, 0, 255), 2) 
        cv2.rectangle(debug_final, (x,y), (x+w, y+h), (0,255,0), 2)
        cv2.imwrite(os.path.join(debug_dir, f"debug_5_CONTOUR_AND_BOX_{label}.jpg"), debug_final)
        
        return (cx + x, cy + y, w, h)
        
    return (cx + cw//4, cy + ch//4, cw//2, ch//2)

def create_roi_masks(img_shape, rect):
    """
    Returns the ROI_1, ROI_2 AND ROI_3
    """
    h_img, w_img = img_shape[:2]
    x1, y1, w1, h1 = rect
    area1 = w1 * h1
    
    mask_roi1 = np.zeros((h_img, w_img), dtype=bool)
    mask_roi1[y1:y1+h1, x1:x1+w1] = True
    
    aspect = float(w1) / float(h1)
    target_area2 = TARGET_AREA_MULTIPLIER * area1
    
    w2 = int(round(math.sqrt(target_area2 * aspect)))
    h2 = int(round(math.sqrt(target_area2 / aspect)))
    
    cx, cy = x1 + (w1 / 2.0), y1 + (h1 / 2.0)
    x2, y2 = int(round(cx - w2 / 2.0)), int(round(cy - h2 / 2.0))
    
    x2, y2 = max(0, x2), max(0, y2)
    if x2 + w2 > w_img: w2 = w_img - x2
    if y2 + h2 > h_img: h2 = h_img - y2
    
    mask_roi2 = np.zeros((h_img, w_img), dtype=bool)
    mask_roi2[y2:y2+h2, x2:x2+w2] = True
    
    mask_roi3 = mask_roi2 & ~mask_roi1 # This is bezel area
    return mask_roi1, mask_roi3, (x2, y2, w2, h2)

def process_image(geom_img_path, metric_img_path, actual_filename, left_label, right_label, processed_dir, debug_dir, bezel_dir, debug=False):
    """
    Uses geom_img_path (whitewashed or raw) to find boundaries.
    In some cases there was high background in images.
    As a result the handmade rectangles on skin
    could not be found. So we made a copy of these images, opened them in GIMP and painted around the
    hand drawn sharpie rectangle with white color. These were called 'whitewashed'images.
    These whitewashed images were opened for finding ROIs.
    But for makig all measurements, the corresponding 'RAW'/'ORIGINAL' images were opened and used.
    
    Uses metric_img_path (strictly raw/original) to extract data statistics.
    """
    img_geom = cv2.imread(geom_img_path)
    print(geom_img_path)
    img_metrics = cv2.imread(metric_img_path)

    if img_geom is None or img_metrics is None:
        print(f"Could not load images for: {actual_filename}")
        return None
        
    print(f"Processing Geometry: {os.path.basename(geom_img_path)} | Extracting Metrics From: {os.path.basename(metric_img_path)}")
    debug_img = img_metrics.copy() # non whitewashed original image
    
    # Calculate maps on the METRIC image (the raw image)
    
    # WE USED THREE METRICS FOR ERYTHEMA   
    ei_map = calculate_erythema_index(img_metrics)
    a_star_map = get_a_star_channel(img_metrics)
    magenta_map = get_magenta(img_metrics)
    
    # Find master box on the GEOMETRY image (the whitewashed/raw image)
    # This is hand drawn rectangle(s) using sharpie
    master_box = find_master_bounding_box(img_geom)
    if not master_box:
        print("  -> Could not detect the master marker boundary.")
        return None
    
    
    bx, by, bw, bh = master_box
    cv2.rectangle(debug_img, (bx, by), (bx+bw, by+bh), (0, 255, 255), 2)
    cv2.putText(debug_img, "Master Bounds", (bx, by-10), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 255), 2)
    metrics = {"Filename": actual_filename, "Status": "Success"}

    
    # --- SINGLE SITE/BOX DETECTION ---
    # For week 1 IP the bMN insertion was not recorded.
    # So this part of script analyzes those images
    is_single_site = "single_" in geom_img_path.lower()
    
    if is_single_site:
        print(f"  -> Single site detected for {actual_filename}.")
        single_safe = get_safe_zone((bx, by, bw, bh), SAFE_ZONE_MARGIN)
        if debug:
            cv2.rectangle(debug_img, (single_safe[0], single_safe[1]), 
                          (single_safe[0]+single_safe[2], single_safe[1]+single_safe[3]), (255, 0, 255), 1)
        
        # Process only the one region, dictating its label from the CSV
        # label means cMN or bMN
        regions = [(right_label, single_safe)]
        
        # Pre-fill the missing side's data with zeros to keep the CSV aligned
        metrics[f"{left_label}_Delta_EI"] = 0.0
        metrics[f"{left_label}_Delta_a"] = 0.0
        metrics[f"{left_label}_Delta_magenta"] = 0.0
        metrics[f"{left_label}_Footprint_Area(px^2)"] = 0
        metrics[f"{left_label}_Aspect_Ratio(px/px)"] = 0.0
        metrics[f"{left_label}_Total_Holes"] = 0
        
    else:
        # --- NORMAL DUAL-SITE INSERTION SPLIT ---
        # This is the case for majority of the images
        # where there are two rectangles hand-drawn with sharpie
        # and there is a bMN or cMN site in each rectangle
        half_w = bw // 2
        left_rect = (bx, by, half_w, bh)
        right_rect = (bx + half_w, by, bw - half_w, bh)
        
        cv2.line(debug_img, (bx + half_w, by), (bx + half_w, by + bh), (0, 255, 255), 2)
        
        left_safe = get_safe_zone(left_rect, SAFE_ZONE_MARGIN)
        right_safe = get_safe_zone(right_rect, SAFE_ZONE_MARGIN)
        
        if debug:
            for zx, zy, zw, zh in [left_safe, right_safe]:
                cv2.rectangle(debug_img, (zx, zy), (zx+zw, zy+zh), (255, 0, 255), 1)
            
        regions = [(left_label, left_safe), (right_label, right_safe)]
        
        

    # Create a fresh, clean copy just for the presentation graphic
    presentation_img = img_metrics.copy()
    
    
    for label, safe_zone in regions:
        # Find microneedle site strictly on the geometry image
        site_rect = find_best_site_in_zone(img_geom, safe_zone, label, debug_dir)
        
        if site_rect:
            x1, y1, w1, h1 = site_rect
            mask_roi1, mask_roi3, rect2 = create_roi_masks(img_geom.shape, site_rect)
            x2, y2, w2, h2 = rect2
            
            # --- 1) CHROMATIC METRICS (From Actual Image) ---
            # Three EI metrics were measured
            # Informal analysis showed all three gave similar conclusion
            # For the manuscript we used the delta_ei metric
            delta_ei = np.mean(ei_map[mask_roi1]) - np.mean(ei_map[mask_roi3])
            delta_a = np.mean(a_star_map[mask_roi1]) - np.mean(a_star_map[mask_roi3])
            delta_magenta = np.mean(magenta_map[mask_roi1]) - np.mean(magenta_map[mask_roi3])
            
            # --- 2 & 3) GEOMETRIC & SHEAR METRICS ---
            # We tried to measure geometric footprint
            # But it was not clear what the relationship was with skin inflammation
            # because skin stretches when microneedle is pushed and this skews the result
            # Therefore, although computed and recorded, we did not use them for analysis.
            footprint_area = w1 * h1
            aspect_ratio = float(w1) / float(h1) if h1 > 0 else 1.0
            
            # --- 4) TEXTURAL ANALYSIS (From Actual Image) ---
            # This section was for computing the number of holes 
            # generated from individual microneedles.
            # However, we found that the results did not match 
            # the actual holes counted manually as a spot check for few images.
            # So we did not use this data for analysis
            green_roi = img_metrics[y1:y1+h1, x1:x1+w1]
            green_gray = cv2.cvtColor(green_roi, cv2.COLOR_BGR2GRAY)
            
            adaptive_holes = cv2.adaptiveThreshold(green_gray, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, 
                                                   cv2.THRESH_BINARY_INV, 15, 3)
            adaptive_holes = cv2.morphologyEx(adaptive_holes, cv2.MORPH_OPEN, np.ones((2, 2), np.uint8))
            
            blur_size = max(3, int(w1 // 8))
            if blur_size % 2 == 0: blur_size += 1
            density_map = cv2.GaussianBlur(adaptive_holes.astype(float), (blur_size, blur_size), 0)
            
            if np.max(density_map) > 0:
                density_map = (density_map / np.max(density_map) * 255).astype(np.uint8)
            else:
                density_map = density_map.astype(np.uint8)
                
            min_needle_spacing = max(3, int(w1 // 16))
            peaks = peak_local_max(density_map, min_distance=min_needle_spacing, threshold_abs=40)
            total_holes = len(peaks)
            
            # heatmap_color = cv2.applyColorMap(density_map, cv2.COLORMAP_JET)
            heatmap_color = cv2.applyColorMap(density_map, cv2.COLORMAP_PINK)

            visual_heatmap = cv2.addWeighted(green_roi, 0.5, heatmap_color, 0.5, 0)
            
            cv2.imwrite(os.path.join(debug_dir, f"textural_hole_heatmap_{label}.jpg"), visual_heatmap)
            
            
            # --- MAGENTA INFLAMMATION HEATMAP ---
            # 1. Crop the already-calculated magenta map to ROI footprint
            magenta_roi = magenta_map[y1:y1+h1, x1:x1+w1]
            magenta_smooth = magenta_roi
            # 2. Apply COOL colormap
            magenta_color_map = cv2.applyColorMap(magenta_smooth, cv2.COLORMAP_COOL)
            
            # 3. Blend it transparently over the raw green_roi image
            visual_magenta_heat = cv2.addWeighted(green_roi, 0.5, magenta_color_map, 0.5, 0)
            
            # 4. Save it to the debug folder
            cv2.imwrite(os.path.join(debug_dir, f"magenta_inflammation_heatmap_{label}.jpg"), visual_magenta_heat)
            
            magenta_color = (255, 0, 255)
            blue_color = (255, 0, 0)
            black_color = (0, 0, 0)
            
            # --- Draw the boxes on the main full-size output image ---
            cv2.rectangle(debug_img, (x1, y1), (x1+w1, y1+h1), blue_color, 2)
            cv2.rectangle(debug_img, (x2, y2), (x2+w2, y2+h2), magenta_color, 2)
            
            # --- Draw Cropped Site with Blue and magenta boxes ---
            
            # Create a fresh copy of the raw image for this specific crop
            presentation_img = img_metrics.copy()
            
            # 1. Draw the exact Green Footprint (ROI 1)
            cv2.rectangle(presentation_img, (x1, y1), (x1+w1, y1+h1), blue_color, 2) 
            
            # 2. Draw the exact Blue Bezel (ROI 2) at its true mathematical location
            cv2.rectangle(presentation_img, (x2, y2), (x2+w2, y2+h2), magenta_color, 2)
            
            # 3. Calculate padded boundaries FOR THE CROP ONLY (captures extra skin outside the box)
            pad = 50  # Increased to 50 pixels so the margin is easily visible
            crop_y1 = max(0, y2 - pad)
            crop_y2 = min(presentation_img.shape[0], y2 + h2 + pad)
            crop_x1 = max(0, x2 - pad)
            crop_x2 = min(presentation_img.shape[1], x2 + w2 + pad)
            
            # 4. Crop the image cleanly using the padded boundary
            cropped_bezel = presentation_img[crop_y1:crop_y2, crop_x1:crop_x2]
            
            # 5. Save the cropped ROI to the bezel folder
            # These are the images reported in manuscript and supplementary data
            # It removes identificable information of pig subjects and humans
            bezel_output_path = os.path.join(bezel_dir, f"bezelled_{label}_{actual_filename}")
            cv2.imwrite(bezel_output_path, cropped_bezel)        
            

            print(f"  {label}:")
            print(f"    Delta EI: {delta_ei:.4f} | Delta a*: {delta_a:.4f} | Delta Magenta: {delta_magenta:.4f}")
            print(f"    Size Footprint: {footprint_area} px^2 | Application Aspect Ratio: {aspect_ratio:.2f} px/px")
            
            metrics[f"{label}_Delta_EI"] = round(delta_ei, 4)
            metrics[f"{label}_Delta_a"] = round(delta_a, 4)
            metrics[f"{label}_Delta_magenta"] = round(delta_magenta, 4)
            metrics[f"{label}_Footprint_Area(px^2)"] = footprint_area
            metrics[f"{label}_Aspect_Ratio(px/px)"] = round(aspect_ratio, 4)
            metrics[f"{label}_Total_Holes"] = total_holes
            
        else:
            # Microneedle insertion site not found
            print(f"  {label}: No insertion site found in safe zone.")
            metrics[f"{label}_Delta_EI"] = 0.0
            metrics[f"{label}_Delta_a"] = 0.0
            metrics[f"{label}_Delta_magenta"] = 0.0
            metrics[f"{label}_Footprint_Area(px^2)"] = 0
            metrics[f"{label}_Aspect_Ratio(px/px)"] = 1.0
            metrics[f"{label}_Total_Holes"] = 0

    final_output_path = os.path.join(processed_dir, f"processed_{actual_filename}")
    cv2.imwrite(final_output_path, debug_img)
    
    print(f"  -> Saved final ROI map: {final_output_path}\n")
    return metrics

def batch_process_folder(input_folder, mapping_file, master_csv_path):
    """Handles mapping dynamic filenames and dual image pipeline, saving local and master CSVs."""
    processed_image_dir = os.path.join(input_folder, "processed_image")
    os.makedirs(processed_image_dir, exist_ok=True)
    
    # --- Create Bezel folder ---
    bezel_dir = os.path.join(input_folder, "bezelled_images")
    os.makedirs(bezel_dir, exist_ok=True)

    mapping = load_side_mapping(mapping_file)

    valid_exts = ('.jpg', '.jpeg', '.png', '.tif', '.tiff')
    all_files = [f for f in os.listdir(input_folder) if f.lower().endswith(valid_exts) and not f.startswith("processed_")]
    
    if not all_files:
        print(f"No images found in {input_folder}")
        return

    files_to_process = []
    for f in all_files:
        if f.startswith("whitewashed_"):
            actual = f.replace("whitewashed_", "", 1)
            files_to_process.append((f, actual))
            
        elif f.startswith("single_"):
            # for week 1 IP case
            true_name = f.replace("single_", "", 1) # Safer than f[7:]
            if true_name in all_files:
                files_to_process.append((f, true_name))
            else:
                raise FileNotFoundError(f"No paired raw file found for ...> {f}")
                
        else:
            if f"whitewashed_{f}" not in all_files and f"single_{f}" not in all_files:
                files_to_process.append((f, f))

    print(f"\n--- Starting Folder: {input_folder} ---")
    print(f"Found {len(files_to_process)} target processing combinations. Initializing tree paths...")

    # Name the local folder CSV based on the folder name (e.g., Week_3_results.csv)
    folder_name = os.path.basename(os.path.normpath(input_folder))
    local_csv_path = os.path.join(input_folder, f"{folder_name}_results.csv")

    fieldnames = ['Week', 'Pig_Name', 'Filename', 'cMN_Delta_EI', 'cMN_Delta_a', 'cMN_Delta_magenta',
                  'cMN_Footprint_Area(px^2)', 'cMN_Aspect_Ratio(px/px)', 'cMN_Total_Holes',
                  'bMN_Delta_EI', 'bMN_Delta_a', 'bMN_Delta_magenta',
                  'bMN_Footprint_Area(px^2)', 'bMN_Aspect_Ratio(px/px)', 'bMN_Total_Holes', 'Status']

    # Open both the local CSV (write mode to overwrite just this folder's old data) 
    # and the Master CSV (append mode to add to the running list)
    with open(local_csv_path, mode='w', newline='') as local_csv, \
         open(master_csv_path, mode='a', newline='') as master_csv:
        
        local_writer = csv.DictWriter(local_csv, fieldnames=fieldnames)
        master_writer = csv.DictWriter(master_csv, fieldnames=fieldnames)
        
        local_writer.writeheader()
        
        # Check if the master file needs headers (only if it's completely empty)
        if os.stat(master_csv_path).st_size == 0:
            master_writer.writeheader()
        
        for geom_filename, metric_filename in files_to_process:
            left_label = mapping.get(metric_filename)
            if not left_label:
                print(f"Skipping {metric_filename} -> Could not find file mapping in the CSV.")
                continue
                
            right_label = "bMN" if left_label == "cMN" else "cMN"
            # In manuscript bMN == ecMN and cMN == ppcMN
            geom_path = os.path.join(input_folder, geom_filename)
            metric_path = os.path.join(input_folder, metric_filename)
            
            base_filename = os.path.splitext(metric_filename)[0]
            specific_debug_dir = os.path.join(input_folder, "debug_images", base_filename)
            os.makedirs(specific_debug_dir, exist_ok=True)
            

            result_row = process_image(geom_path, metric_path, metric_filename, left_label, right_label, 
                                       processed_dir=processed_image_dir, debug_dir=specific_debug_dir, 
                                       bezel_dir=bezel_dir, debug=True) 
            
            if result_row:
                # Assign the folder name to the 'Week' column
                week_number = folder_name.split("_")
                result_row['Week'] = week_number[1]
                pig_name = metric_filename.split("_")[1]
                result_row["Pig_Name"] = pig_name

                local_writer.writerow(result_row)
                master_writer.writerow(result_row)
                
    print(f"Finished {folder_name}. Local summary saved to: {local_csv_path}.")


#########################################################
#########################################################
# --- RUN THE SCRIPT ---
folders = ['Week_1', 'Week_2', 'Week_3', 'Week_4', 'Week_5', 'Week_6']
csv_mapping_file = "./Image_bMN_cMN_MAPPING.csv" 
master_output_file = "final_experiment_results.csv"

# Clear out the old master file at the start of a fresh run so data doesn't duplicate
open(master_output_file, 'w').close() 

for folder in folders:
    input_folder = f"./{folder}"
    # Check if the folder exists before trying to process it
    if os.path.exists(input_folder):
        batch_process_folder(input_folder, mapping_file=csv_mapping_file, master_csv_path=master_output_file)
    else:
        print(f"\nSkipping {folder} - Directory not found.")

print(f"\nALL FOLDERS COMPLETE. Master cumulative log saved to: {master_output_file}")