Object-Based Image Analysis

Practical Exam
Arunima Sen

---

Based on the existing tutorial notebooks, here's what I added:

1. I applied the Sobel filter to detect edges in the artificial landscape
2. I used thresholding to create masks that separate areas based on pixel intensity
3. I added k-means clustering to segment the image into different regions based on color, testing how the segmentation performed with different values of k

I worked with two images: Artifical_landscape_UTM_012025.tif with 3 bands (RGB) and ortho_subset_I.tif with 4 bands (RGB + NIR). The code was adapted to handle both cases by ensuring proper band extraction and normalization. Adjustments were made to visualization and processing steps to accommodate the additional NIR band where applicable.

Setup

# Check if running on Google Colab
if 'google.colab' in str(get_ipython()):
    import os
    repo_dir = "obia_tutorials"
    marker_file = os.path.join(repo_dir, ".setup_done")

    # Setup the environment only if it hasn't been done already
    if not os.path.exists(marker_file):
        # Clone the repository
        !git clone https://github.com/mariarodriguezn/obia_tutorials.git

        # Install the required packages
        !pip install -r obia_tutorials/requirements.txt

        # Create a marker file to avoid re-running the setup
        with open(marker_file, 'w') as f:
            f.write("Setup completed")

# Faster way to access the .tif file :)
!git clone https://github.com/arunima-sen/OBIA.git

# imports
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time
import rasterio

from scipy import stats
from skimage.color import label2rgb
from skimage.measure import regionprops, regionprops_table, perimeter
from skimage.segmentation import mark_boundaries, slic
from skimage.util import map_array
from sklearn.ensemble import RandomForestClassifier
from tqdm import tqdm

# extra imports
from scipy.ndimage import sobel
from skimage import io, filters, measure, morphology
from skimage.segmentation import mark_boundaries
from sklearn.cluster import KMeans

1. Image with 3 bands

Image Segmentation

# File path to the image
img_path1 = "OBIA/Artifical_landscape_UTM_012025.tif"

# Check the no. of bands
with rasterio.open(img_path1) as src:
    num_bands = src.count
    print(f"Number of bands: {num_bands}")

Number of bands: 3

# Read the image and extract the bands
with rasterio.open(img_path1) as src:
    # Extract RGB bands and normalize to [0, 1]
    bands = src.read([1, 2, 3]).astype(float) / 255
    red, green, blue = bands

# Stack RGB into an array
rgb = np.stack([red, green, blue], axis=-1)

# Stack all bands (only RGB in this case) into an array
image = np.stack([red, green, blue], axis=-1)

# Set the fixed compactness and number of superpixels
compactness = 0.3
n_segments = 1000

# Generate segments using SLIC
segments = slic(image, n_segments=n_segments, compactness=compactness, start_label=1)

# Create the figure with 3 subplots
fig, axs = plt.subplots(ncols=3, figsize=(15, 5), constrained_layout=True)

# Display the original RGB image
axs[0].imshow(rgb)
axs[0].set_title("Original RGB")

# Display the RGB image with segmentation boundaries for 3000 segments
axs[1].imshow(mark_boundaries(rgb, segments, color=(1, 0, 0), mode="thick"))
axs[1].set_title(f"SLIC Segmentation Boundaries ({n_segments} segments)")

# Display the mean RGB per segment
axs[2].imshow(label2rgb(segments, rgb, kind='avg'))
axs[2].set_title("Mean RGB per Segment")

# Remove axis for all subplots
for ax in axs:
    ax.set_axis_off()

# Display the combined figure
plt.show()

1.1 Sobel Edge Detection

# apply the filter to all bands
sobel_red_x = sobel(red, axis=0)
sobel_red_y = sobel(red, axis=1)
sobel_red = np.hypot(sobel_red_x, sobel_red_y)

sobel_green_x = sobel(green, axis=0)
sobel_green_y = sobel(green, axis=1)
sobel_green = np.hypot(sobel_green_x, sobel_green_y)

sobel_blue_x = sobel(blue, axis=0)
sobel_blue_y = sobel(blue, axis=1)
sobel_blue = np.hypot(sobel_blue_x, sobel_blue_y)

# combine from all channels
sobel_combined = (sobel_red + sobel_green + sobel_blue) / 3

# plot
fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# original
axs[0].imshow(np.dstack([red, green, blue]))
axs[0].set_title("Original RGB Image")
axs[0].axis("off")

# sobel
axs[1].imshow(sobel_combined, cmap="Reds")
axs[1].set_title("Sobel Edge Detection")
axs[1].axis("off")

plt.show()

1.2 Thresholding

# RGB channels
red, green, blue = image[:, :, 0], image[:, :, 1], image[:, :, 2]

# mean thresholding
threshold_red = np.mean(red)
threshold_green = np.mean(green)
threshold_blue = np.mean(blue)

# create binary masks
mask_red = red > threshold_red
mask_green = green > threshold_green
mask_blue = blue > threshold_blue

# combine
combined_mask = mask_red & mask_green & mask_blue

# median filter
structuring_element = morphology.square(10)
filtered_mask = filters.median(combined_mask, structuring_element)


fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
masks = [combined_mask, filtered_mask]
titles = ["Initial RGB mask", "RGB mask after median filter"]

for i in range(2):

    axs[i, 0].imshow(masks[i], cmap="gray")
    axs[i, 0].set_title(titles[i])
    axs[i, 0].axis("off")

    # yellow boundary for objects
    rgb_with_boundaries = mark_boundaries(image, masks[i], color=(255, 255, 0), outline_color=(255, 255, 0), mode="thick")
    axs[i, 1].imshow(rgb_with_boundaries)
    axs[i, 1].set_title(f"RGB with boundaries of {titles[i]}")
    axs[i, 1].axis("off")


    labeled_segments = measure.label(masks[i].astype("int"), background=-1)
    n_segments = len(np.unique(labeled_segments))
    axs[i, 2].imshow(labeled_segments, cmap="Spectral", interpolation="nearest")
    axs[i, 2].set_title(f"Resulting segments: {n_segments}")
    axs[i, 2].axis("off")


plt.tight_layout()
plt.show()

WARNING:matplotlib.image:Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
WARNING:matplotlib.image:Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

1.3 K-Means Clustering

# reshape the image into a 2D array
pixels = image.reshape((-1, 3))

# My list of k-values to test
k_values = [4, 6, 10]

fig, axs = plt.subplots(1, 4, figsize=(20, 7))
axs[0].imshow(image)
axs[0].set_title("Original RGB Image")
axs[0].axis('off')

# Iterate over k_values
for i, k in enumerate(k_values):

    kmeans = KMeans(n_clusters=k)
    kmeans.fit(pixels)

    segmentation = kmeans.labels_

    segmentation_image = segmentation.reshape(image.shape[:2])

    axs[i + 1].imshow(segmentation_image, cmap='tab20b', interpolation='nearest')
    axs[i + 1].set_title(f"Segmentation with K={k}")
    axs[i + 1].axis('off')

plt.tight_layout()
plt.show()

/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(
/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(
/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(
<ipython-input-42-263c294a7b18>:16: ConvergenceWarning: Number of distinct clusters (9) found smaller than n_clusters (10). Possibly due to duplicate points in X.
  kmeans.fit(pixels)

2. Image with 4 bands

Image Segmentation

# File path to the image
img_path2 = "obia_tutorials/sample_data/ortho_subset_I.tif"

# Check the no. of bands
with rasterio.open(img_path) as src:
    num_bands = src.count
    print(f"Number of bands: {num_bands}")

Number of bands: 3

# File path to the image
img_path2 = "obia_tutorials/sample_data/ortho_subset_I.tif"

# Read the image and extract the bands
with rasterio.open(img_path2) as src:
    # Extract red, green, blue, and NIR bands, normalize to [0, 1]
    bands = src.read([1, 2, 3, 4]).astype(float) / 255
    red, green, blue, nir = bands

# Stack RGB into an array
rgb = np.stack([red, green, blue], axis=-1)

# Calculate NDVI
ndvi = (nir - red) / (nir + red + 1e-6)  # small value added to avoid division by zero

# Stack all bands into an array
image = np.stack([red, green, blue, nir], axis=-1)

# Set the fixed compactness and number of superpixels
compactness = 0.3
n_segments = 3000

# Generate segments using SLIC
segments = slic(image, n_segments=n_segments, compactness=compactness, start_label=1)

# Create the figure with 3 subplots
fig, axs = plt.subplots(ncols=3, figsize=(15, 5), constrained_layout=True)

# Display the original RGB image
axs[0].imshow(rgb)
axs[0].set_title("Original RGB")

# Display the RGB image with segmentation boundaries for 3000 segments
axs[1].imshow(mark_boundaries(rgb, segments, color=(1, 0, 0), mode="thick"))
axs[1].set_title(f"SLIC Segmentation Boundaries ({n_segments} segments)")

# Display the mean RGB per segment
axs[2].imshow(label2rgb(segments, rgb, kind='avg'))
axs[2].set_title("Mean RGB per Segment")

# Remove axis for all subplots
for ax in axs:
    ax.set_axis_off()

# Display the combined figure
plt.show()

2.1 Sobel Edge Detection

# apply the filter to all bands
sobel_red = np.hypot(sobel(red, axis=0), sobel(red, axis=1))
sobel_green = np.hypot(sobel(green, axis=0), sobel(green, axis=1))
sobel_blue = np.hypot(sobel(blue, axis=0), sobel(blue, axis=1))
sobel_nir = np.hypot(sobel(nir, axis=0), sobel(nir, axis=1))

# combine from all channels
sobel_combined = (sobel_red + sobel_green + sobel_blue + sobel_nir) / 4

# plot
fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# original
axs[0].imshow(np.dstack([red, green, blue]))
axs[0].set_title("Original RGB Image")
axs[0].axis("off")

# sobel
axs[1].imshow(sobel_combined, cmap="Reds")
axs[1].set_title("Sobel Edge Detection")
axs[1].axis("off")

plt.show()

2.2 Thresholding

# extract bands
red, green, blue, nir = image[:, :, 0], image[:, :, 1], image[:, :, 2], image[:, :, 3]

# mean thresholding for each channel
threshold_red = np.mean(red)
threshold_green = np.mean(green)
threshold_blue = np.mean(blue)
threshold_nir = np.mean(nir)

# create binary masks for each channel
mask_red = red > threshold_red
mask_green = green > threshold_green
mask_blue = blue > threshold_blue
mask_nir = nir > threshold_nir

combined_mask = mask_red & mask_green & mask_blue & mask_nir

# apply a median filter
structuring_element = morphology.square(10)
filtered_mask = filters.median(combined_mask, structuring_element)

# plot
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
masks = [combined_mask, filtered_mask]
titles = ["Initial RGB+NIR mask", "RGB+NIR mask after median filter"]

for i in range(2):
    axs[i, 0].imshow(masks[i], cmap="gray")
    axs[i, 0].set_title(titles[i])
    axs[i, 0].axis("off")

    rgb_with_boundaries = mark_boundaries(rgb, masks[i], color=(1, 0, 0), outline_color=(1, 0, 0), mode="thicker")
    axs[i, 1].imshow(rgb_with_boundaries)
    axs[i, 1].set_title(f"RGB with boundaries of {titles[i]}")
    axs[i, 1].axis("off")

    labeled_segments = measure.label(masks[i].astype("int"), background=-1)
    n_segments = len(np.unique(labeled_segments))
    axs[i, 2].imshow(labeled_segments, cmap="Spectral", interpolation="nearest")
    axs[i, 2].set_title(f"Resulting segments: {n_segments}")
    axs[i, 2].axis("off")

plt.tight_layout()
plt.show()

2.3 K-Means Clustering

# reshape the image into a 2D array
pixels = image.reshape((-1, 4))  # now using 4 bands (R, G, B, NIR)

# My list of k-values to test
k_values = [4, 6, 10]

fig, axs = plt.subplots(1, 4, figsize=(20, 7))
axs[0].imshow(rgb)
axs[0].set_title("Original RGB Image")
axs[0].axis('off')

# Iterate over k-values
for i, k in enumerate(k_values):
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(pixels)

    segmentation = kmeans.labels_

    segmentation_image = segmentation.reshape(image.shape[:2])

    axs[i + 1].imshow(segmentation_image, cmap='tab20b', interpolation='nearest')
    axs[i + 1].set_title(f"Segmentation with K={k}")
    axs[i + 1].axis('off')

plt.tight_layout()
plt.show()

/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(
/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(
/usr/local/lib/python3.11/dist-packages/sklearn/cluster/_kmeans.py:870: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  warnings.warn(