Dividir y conquistar

A veces, los programas se vuelven muy largos y difíciles de leer. Hablamos de código espagueti. Una forma de hacer que el código sea más fácil de leer y mantener es dividirlo en funciones más pequeñas y simplemente usarlas en flujos de trabajo más complejos. El principio de diseño de software se llama Dividir y conquistar.

from skimage.io import imread
from skimage.morphology import white_tophat, disk
from skimage.filters import gaussian, threshold_otsu
from skimage.measure import label, regionprops_table
import pandas as pd
import numpy as np

image = imread("../../data/blobs.tif")
footprint = disk(15)
background_subtracted = white_tophat(image, 
                                     footprint=footprint)
particle_radius = 5
denoised = gaussian(background_subtracted, 
                    sigma=particle_radius)
binary = denoised > threshold_otsu(denoised)
labels = label(binary)
requested_measurements = ["label", "area", "mean_intensity"]
regionprops = regionprops_table(image, 
                                labels, 
                                properties=requested_measurements)
table = pd.DataFrame(regionprops)
mean_total_intensity = np.mean(table["area"] * table["mean_intensity"])
mean_total_intensity

17136.90322580645

Una forma común y fácil de hacer que dicho código sea más fácil de leer es dividirlo en secciones que comiencen cada una con un comentario.

# configuration
file_to_analyze = "../../data/blobs.tif"
background_subtraction_radius = 15
particle_radius = 5
requested_measurements = ["area", "mean_intensity"]

# load data
image = imread(file_to_analyze)

# preprocess image
footprint = disk(background_subtraction_radius)
background_subtracted = white_tophat(image, 
                                     footprint=footprint)
denoised = gaussian(background_subtracted, 
                    sigma=particle_radius)

# segment image
binary = denoised > threshold_otsu(denoised)
labels = label(binary)

# extract features
regionprops = regionprops_table(image, 
                                labels, 
                                properties=requested_measurements)
table = pd.DataFrame(regionprops)

# descriptive statistics
mean_total_intensity = np.mean(table["area"] * table["mean_intensity"])
mean_total_intensity

17136.90322580645

Sería más profesional poner todo este código en subrutinas significativas y llamarlas desde una función central.

# reusable functions
def preprocess_image(image, background_subtraction_radius, particle_radius):
    """Apply background removal and denoising"""
    footprint = disk(background_subtraction_radius)
    background_subtracted = white_tophat(image, footprint=footprint)
    denoised = gaussian(background_subtracted, sigma=particle_radius)
    return denoised

def segment_image(image):
    """Apply thresholding and connected component analysis"""
    binary = image > threshold_otsu(image)
    labels = label(binary)
    return labels

def extract_features(image, labels, requested_measurements):
    """Measure specified properties"""
    regionprops = regionprops_table(image, 
                                    labels, 
                                    properties=requested_measurements)
    table = pd.DataFrame(regionprops)
    return table

Después de poner grupos de pasos de procesamiento en funciones, podemos llamarlos desde una función principal. Esta función puede reutilizarse más tarde para aplicar el mismo procedimiento a todas las imágenes en una carpeta. También sirve como índice, una visión general del flujo de trabajo de procesamiento de imágenes. Al leer solo esta función, conocemos todos los pasos de procesamiento y qué parámetros tienen.

def analyse_average_total_intensity(filename, 
                                    background_subtraction_radius = 15, 
                                    particle_radius = 5):
    """Load an image, segment objects and measure their mean total intensity."""
    image = imread(filename)
    denoised = preprocess_image(image, 
                                background_subtraction_radius, 
                                particle_radius)
    labels = segment_image(denoised)
    requested_measurements = ["area", "mean_intensity"]
    table = extract_features(image, 
                             labels, 
                             requested_measurements)

    # descriptive statistics
    mean_total_intensity = np.mean(table["area"] * table["mean_intensity"])
    
    return mean_total_intensity

# configuration
file_to_analyze = "../../data/blobs.tif"

Esta función central se puede leer fácilmente; tiene solo 6 líneas de código

analyse_average_total_intensity(file_to_analyze)

17136.90322580645

Esta función puede entonces reutilizarse también para otros archivos de imagen.

analyse_average_total_intensity("../../data/BBBC007_batch/20P1_POS0005_D_1UL.tif")

884.2620087336245