Confocal_40x_ROIs_Pdgfrb_Morpho_BatchScript.qmd

---
title-block-banner: true
title: "Morphological analysis of PDGFR-β+ cells in defined ROIs"
subtitle: "Batch processing Python script"
date: today
date-format: full
author: 
  - name: "Daniel Manrique-Castano"
    orcid: 0000-0002-1912-1764
    degrees:
      - PhD
    affiliation: 
      - name: Univerisity Laval 
        department: Psychiatry and Neuroscience
        group: Laboratory of neurovascular interactions 
note: "GitHub: https://daniel-manrique.github.io/"
keywords: 
  - Ki67
  - PDGFR-β
  - Brain injury
  - Cell proliferation
  - Bayesian modeling 
   
license: "CC BY"

format:
   pdf: 
    toc: true
    number-sections: true
    colorlinks: true
   html:
    code-fold: true
    embed-resources: true
    toc: true
    toc-depth: 2
    toc-location: left
    number-sections: true
    theme: spacelab

knitr:
  opts_chunk: 
    warning: false
    message: false
    
csl: science.csl
bibliography: references.bib
---

```{python}
import os
import pandas as pd
import numpy as np
import skimage
from skimage.io import imread, imsave
from skimage.filters import gaussian
from skimage import segmentation, morphology
from skimage.measure import label, regionprops_table
from skimage.feature import peak_local_max
from scipy import ndimage as ndi
from skan import draw, Skeleton, summarize
from skan.csr import skeleton_to_csgraph
from skimage.segmentation import watershed
from scipy import ndimage as ndi
import matplotlib.pyplot as plt
from skimage.util import invert
from skimage.filters import threshold_otsu
from skimage.filters import unsharp_mask
from skimage.exposure import equalize_adapthist
from skimage import color
from skimage import util
from skimage import transform

# Specify your directory
dir_path = "D:/Research/Stroke_PDGFR-B_Reactivity/Images_Raw/Confocal_40x_ROIs_CD31-Pdgfrb-CD13/Images_Cells"
save_dir = "D:/Research/Stroke_PDGFR-B_Reactivity/Images_Raw/Confocal_40x_ROIs_CD31-Pdgfrb-CD13/Images_Morphology"

# Initialize an empty DataFrame for the results
all_results = pd.DataFrame()

# Initialize an empty DataFrame for the branch data
all_branch_data = pd.DataFrame()

# Loop over all .tif files in the directory
for filename in os.listdir(dir_path):
    if filename.endswith(".tif"):
        # Construct the full file path
        image_path = os.path.join(dir_path, filename)

        # Load the image
        Raw = imread(image_path)
        
        # Create a subdirectory for this image's results
        image_save_dir = os.path.join(save_dir, filename.replace('.tif', ''))
        os.makedirs(image_save_dir, exist_ok=True)

        # Perform the processing steps...

        Smooth = gaussian(Raw, sigma=5)
        Unsharp = unsharp_mask(Smooth , radius=10, amount=2)
        #Sigmoid = adjust_sigmoid(Smooth, gain=2)
        #Segmentation = segmentation.morphological_chan_vese(Sigmoid, num_iter=30, smoothing=1)
        
        # Applying Clahe.
        Clahe = equalize_adapthist(Unsharp , clip_limit = 0.02)
        # Rescaling img2 from 0 to 255.
        Clahe = Clahe*255.0

        # Apply otsu
        Otsu = threshold_otsu(Clahe)
        # Pixels with intensity greater than the "threshold" are kept.
        Otsu = 255*(Clahe  > Otsu)
      
        # Closing  
        Closing = skimage.morphology.isotropic_closing(Otsu, radius=2)
        
        # Reduce Image
        scale_factor = 0.5  # reduce the size by 50%
        inverted_closing = util.invert(Closing)
        Closing_size = (np.array(inverted_closing.shape) * scale_factor).astype(int)
        # Resize the image
        Closing_PNG = transform.resize(inverted_closing, Closing_size)
      
        #Holes = morphology.remove_small_holes(Segmentation, 5 ** 3)
        Objects = morphology.remove_small_objects(Closing, min_size=500)
        inverted_Objects = util.invert(Objects)
        Objects_size = (np.array(inverted_Objects.shape) * scale_factor).astype(int)
        # Resize the image
        Objects_PNG = transform.resize(inverted_Objects, Objects_size)
        
        
        #Objects_Inv = np.invert(Objects)
        
        #Labels = label(Objects)
        
        # Compute the distance transform of the binary image
        distance = ndi.distance_transform_edt(Objects)

       # Find the local maxima of the distance transform
        coordinates = peak_local_max(distance, min_distance=400, labels=Objects)

       # Create an image with these local maxima as seeds
        seeds = np.zeros(distance.shape, dtype=bool)
        seeds[tuple(coordinates.T)] = True
        seeds = ndi.label(seeds)[0]

       # Apply the watershed algorithm
        Labels = watershed(-distance, seeds, mask=Objects)
        
        imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Labels.tif')), Labels)

        # Initialize a list to hold the filenames
        object_filenames = []

        # Loop over each label
        for i in range(1, Labels.max() + 1):
            # Create a new image containing only the current label
            single_object = (Labels == i)

            # Construct the filename for this object
            object_filename = f"{filename.replace('.tif', '')}_object{i}.tif"

            # Save the image
            imsave(os.path.join(image_save_dir, object_filename), single_object)

            # Append the filename to the list
            object_filenames.append(object_filename)
       

    
        # Save the processed images and individual objects
        #imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Raw.tif')), Raw)
        #imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Smooth.tif')), Smooth)
        #imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Undharp.tif')), Unsharp)
        imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Closing.png')), Closing_PNG)
        imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Objects.png')), Objects_PNG)
        #imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Labels.tif')), Labels)
        #imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Labels.tif')), Watershed_Labels)


        # Compute the properties for this image
        props = regionprops_table(Labels, properties=('image', 'perimeter', 'solidity', 'centroid', 'area', 'convex_area', 'eccentricity', 'euler_number', 'feret_diameter_max', 'axis_major_length', 'axis_minor_length')) 

        # Create a DataFrame and add the image name as the first column
        Results = pd.DataFrame(props)
        Results.insert(0, 'Image_Name', filename)
        
        # Add the filenames as a new column in the DataFrame
        Results['Object_Filename'] = object_filenames

        # Append the results for this image to the overall results
        all_results = pd.concat([all_results, Results])

        # Perform the skeleton analysis
        Img_Skeleton = morphology.skeletonize(Labels)
        # Invert the binary image to make the skeleton lines black on a white background
        inverted_skeleton = util.invert(Img_Skeleton)
        # Convert the binary image to an RGB image
        rgb_skeleton = color.gray2rgb(inverted_skeleton)
        
        imsave(os.path.join(image_save_dir, filename.replace('.tif', '_Img_Skeleton.tif')), rgb_skeleton)

        pixel_graph, coordinates = skeleton_to_csgraph(Img_Skeleton)

        Branch_data = summarize(Skeleton(Img_Skeleton))
        Branch_data.insert(0, 'Image_Name', filename)
        all_branch_data = pd.concat([all_branch_data, Branch_data])

       
        Branch_types = draw.overlay_euclidean_skeleton_2d(Unsharp, Branch_data, skeleton_color_source='branch-type')

        # Save the figure
        plt.axis('off') # Optional: remove axis
        plt.savefig(os.path.join(image_save_dir, filename.replace('.tif', '_BranchTypes.png')), bbox_inches='tight', pad_inches=0)


# Save the overall results to a .csv file
all_results.to_csv(os.path.join(save_dir, 'Pdgfrb_Morphology.csv'), index=False)
all_branch_data.to_csv(os.path.join(save_dir, 'Pdgfrb_Skeleton.csv'), index=False)

```

```{r}
sessionInfo()
```