Entropia in telerilevamento

Questa e' un risorsa a cui aggrapparsi quando non c'e' niente nell'informazione spettrale che possa aiutare

In pratica si usa l'entropia dell'immagine per distinguere nella tessitura dell'immagine

 

 

#!/usr/bin/env python3

"""

texture_features.py

Calcola mappe di texture (entropia locale + feature GLCM: contrast, homogeneity,

correlation, energy) da un'immagine grayscale o RGB (es. TIFF), per evidenziare

differenze di struttura spaziale tra oggetti di colore simile.

Pipeline:

    1. Legge l'immagine (TIFF/PNG/JPEG...). Se RGB, converte in luminanza

       (rgb2gray) per isolare la struttura spaziale dal colore.

    2. Normalizza a uint8.

    3. Calcola entropia locale (skimage.filters.rank.entropy) a chunk di righe.

    4. Calcola feature GLCM (contrast, homogeneity, correlation, energy) su

       finestra scorrevole, mediando su 4 angoli (0, 45, 90, 135 gradi) per

       essere invariante all'orientamento della texture.

    5. Salva ogni mappa come GeoTIFF (o TIFF semplice se l'input non ha georef).

Uso:

    python texture_features.py 264150.TIF --out_dir ./texture_out \

        --entropy_radius 5 --glcm_window 15 --glcm_step 4

    --glcm_step > 1 calcola il GLCM ogni N pixel e poi fa upsampling (nearest)

    per velocizzare — il GLCM è l'operazione più costosa dello script.

Dipendenze:

    pip install scikit-image numpy tifffile rasterio  # rasterio opzionale, solo per georef

"""

import argparse

import os

import numpy as np

from skimage.io import imread, imsave

from skimage.color import rgb2gray

from skimage.util import img_as_ubyte

from skimage.filters.rank import entropy as rank_entropy

from skimage.morphology import disk

from skimage.feature import graycomatrix, graycoprops

def load_grayscale(path, verbose=True):

    """Carica l'immagine e la converte in luminanza uint8 se RGB."""

    img = imread(path)

    if verbose:

        print(f"Immagine caricata: shape={img.shape}, dtype={img.dtype}")

    if img.ndim == 3:

        # scarta eventuale canale alpha

        rgb = img[:, :, :3]

        gray = rgb2gray(rgb)

        if verbose:

            print("Immagine RGB -> convertita in luminanza (rgb2gray) per isolare la struttura dal colore.")

    else:

        gray = img.astype(np.float64)

    gray = np.nan_to_num(gray, nan=0.0)

    gmin, gmax = gray.min(), gray.max()

    if gmax - gmin < 1e-12:

        return np.zeros_like(gray, dtype=np.uint8)

    gray_norm = (gray - gmin) / (gmax - gmin)

    return img_as_ubyte(gray_norm)

def compute_entropy_map(img_u8, radius, chunk_rows=512, verbose=True):

    """Entropia locale a chunk di righe, con overlap = radius."""

    n_rows, n_cols = img_u8.shape

    out = np.zeros((n_rows, n_cols), dtype=np.float32)

    footprint = disk(radius)

    for start in range(0, n_rows, chunk_rows):

        end = min(start + chunk_rows, n_rows)

        pad_top = min(radius, start)

        pad_bottom = min(radius, n_rows - end)

        chunk = img_u8[start - pad_top: end + pad_bottom, :]

        ent_chunk = rank_entropy(chunk, footprint)

        out[start:end, :] = ent_chunk[pad_top: pad_top + (end - start), :]

        if verbose:

            print(f"  entropia: righe {start}-{end}/{n_rows} completate")

    return out

def compute_glcm_maps(img_u8, window, step=4, levels=32, verbose=True):

    """

    Calcola contrast/homogeneity/correlation/energy con GLCM su finestra

    scorrevole (mediata su 4 angoli), campionando ogni `step` pixel e poi

    facendo upsampling nearest per velocizzare. Restituisce un dict di mappe

    2D della stessa dimensione dell'input.

    """

    n_rows, n_cols = img_u8.shape

    half = window // 2

    angles = [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4]

    props = ['contrast', 'homogeneity', 'correlation', 'energy']

    rows_idx = list(range(half, n_rows - half, step))

    cols_idx = list(range(half, n_cols - half, step))

    small_maps = {p: np.zeros((len(rows_idx), len(cols_idx)), dtype=np.float32) for p in props}

    total = len(rows_idx)

    for ri, r in enumerate(rows_idx):

        for ci, c in enumerate(cols_idx):

            patch = img_u8[r - half:r + half, c - half:c + half]

            patch_q = (patch.astype(np.float32) / 255 * (levels - 1)).astype(np.uint8)

            glcm = graycomatrix(patch_q, distances=[1], angles=angles,

                                 levels=levels, symmetric=True, normed=True)

            for p in props:

                small_maps[p][ri, ci] = graycoprops(glcm, p).mean()

        if verbose and ri % max(1, total // 10) == 0:

            print(f"  GLCM: riga {ri}/{total}")

    # upsampling nearest alla dimensione originale

    from skimage.transform import resize

    full_maps = {}

    for p in props:

        full_maps[p] = resize(small_maps[p], (n_rows, n_cols), order=0,

                               preserve_range=True, anti_aliasing=False).astype(np.float32)

    return full_maps

def save_map(arr, path, verbose=True):

    """Salva una mappa float32 come TIFF a 32 bit."""

    imsave(path, arr.astype(np.float32))

    if verbose:

        print(f"Salvato: {path}")

def main():

    parser = argparse.ArgumentParser(description="Feature di texture (entropia + GLCM) da immagine grayscale/RGB.")

    parser.add_argument("image_path", help="Path all'immagine (es. 264150.TIF)")

    parser.add_argument("--out_dir", default="./texture_out", help="Cartella di output")

    parser.add_argument("--entropy_radius", type=int, default=5, help="Raggio disk() per entropia locale")

    parser.add_argument("--glcm_window", type=int, default=15, help="Dimensione finestra quadrata per GLCM")

    parser.add_argument("--glcm_step", type=int, default=4,

                         help="Passo di campionamento per il GLCM (>1 = piu' veloce, poi upsampling)")

    parser.add_argument("--glcm_levels", type=int, default=32, help="Numero di livelli di quantizzazione per GLCM")

    parser.add_argument("--chunk_rows", type=int, default=512, help="Righe per chunk nel calcolo dell'entropia")

    args = parser.parse_args()

    os.makedirs(args.out_dir, exist_ok=True)

    print(f"Caricamento immagine: {args.image_path}")

    img_u8 = load_grayscale(args.image_path)

    gray_path = os.path.join(args.out_dir, "grayscale_input.tif")

    save_map(img_u8.astype(np.float32), gray_path)

    print(f"\nCalcolo entropia locale (raggio={args.entropy_radius}) ...")

    ent_map = compute_entropy_map(img_u8, args.entropy_radius, chunk_rows=args.chunk_rows)

    save_map(ent_map, os.path.join(args.out_dir, "entropy.tif"))

    print(f"\nCalcolo feature GLCM (finestra={args.glcm_window}, step={args.glcm_step}) ...")

    glcm_maps = compute_glcm_maps(img_u8, args.glcm_window, step=args.glcm_step, levels=args.glcm_levels)

    for name, arr in glcm_maps.items():

        save_map(arr, os.path.join(args.out_dir, f"glcm_{name}.tif"))

    print("\nCompletato. Mappe generate:")

    print("  - entropy.tif           (disordine generale)")

    print("  - glcm_contrast.tif     (differenze locali di intensita', bordi netti)")

    print("  - glcm_homogeneity.tif  (uniformita' locale)")

    print("  - glcm_correlation.tif  (pattern direzionali/periodici)")

    print("  - glcm_energy.tif       (regolarita'/ripetitivita')")

    print("\nConfronta queste mappe negli stessi punti per separare oggetti con colore simile ma struttura diversa.")

if __name__ == "__main__":

    main()