Describing Images with Feature Maps¶

A feature map turns an image into a new image that highlights one property, such as edges, thin lines, or texture. Feature maps are common inputs for clustering and machine learning.

This notebook computes seven of the 16 available feature maps for three SEM images and compares them.

JuSPICE Modules and Classes¶

juspice.io: load_data
SPICEData.features: built-in feature-extraction accessor (see juspice.feature_extract.FeatureExtractionAccessor)
juspice.feature_extract: FEATURE_MAP_NAMES
juspice.tracking: Tracker

Each spice.features.extract(feature_names=...) call returns a new SPICEData and is recorded in spice.history; the image itself is unchanged.

[1]:

from juspice.tracking import Tracker

tracker = Tracker(
    include_metadata=True, notes='Advanced feature-map extraction notebook'
)

tracker.recording_start()

/Users/amir/GIT_repositories/juspice_pre_release/venvs/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

[2]:

# --------- Block 0: Setup ----------

# Standard library
from pathlib import Path
import sys
import os
import glob
import importlib
import shutil

# Third-party
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# JuSPICE
import juspice.feature_extract
from juspice.io import load_data, save_data
from juspice.feature_extract import FEATURE_MAP_NAMES

# Paths
try:
    notebook_dir = Path(__file__).resolve().parent
except Exception:
    notebook_dir = Path.cwd()

cur = notebook_dir
repo_root = None
for _ in range(6):
    if (cur / 'juspice').exists() or (cur / 'pyproject.toml').exists():
        repo_root = cur
        break
    if cur.parent == cur:
        break
    cur = cur.parent
if repo_root is None:
    repo_root = notebook_dir
if str(repo_root) not in sys.path:
    sys.path.insert(0, str(repo_root))

sample_dir = os.path.join(repo_root, 'Sample_data', 'em')

[3]:

tracker.recording_stop()

Load the images¶

The three images share one history, so all extractions are recorded together.

[4]:

tracker.recording_start()

[5]:

# --------- Block 1: Load images ----------

image_paths = sorted(glob.glob(os.path.join(sample_dir, '*.png')))
if not image_paths:
    raise FileNotFoundError(f'No PNG images found in {sample_dir}')

print(f'Loading {len(image_paths)} images...')

# All images share a single DatasetHistory instance (the first-loaded
# image's history) so every feature extraction below — across every
# image and every loop — accumulates into one unified spice.history.
# `spice` / `img_spice` are just variable names for the SPICEData
# instances load_data() returns — any name would work; we use `spice`
# throughout these notebooks as an intuitive nod to JuSPICE / SPICEData.
spice = None
spice_dict = {}
for path in image_paths:
    img_spice = load_data(path)
    if spice is None:
        spice = img_spice
    else:
        img_spice.history = spice.history
    spice_dict[os.path.basename(path)] = img_spice

feature_results = {}
rows = []
print(f'Loaded {len(spice_dict)} images.')

Loading 3 images...
Loaded 3 images.

[6]:

tracker.recording_stop()

Reloading the module picks up recent code edits; it does not affect the results.

[7]:

tracker.recording_start()

[8]:

# --------- Block 2: Reload module ----------

importlib.reload(juspice.feature_extract)
print('Module reloaded')

Module reloaded

[9]:

tracker.recording_stop()

Smoothing at two scales¶

A Gaussian filter [6] averages each pixel with its neighbours. With sigma = 1, only fine noise is removed.

[10]:

tracker.recording_start()

[11]:

# --------- Block 3: gaussian_sigma1 ----------

feature_name = 'gaussian_sigma1'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='viridis')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting gaussian_sigma1 from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_14_1.png

[12]:

tracker.recording_stop()

With sigma = 3, small details fade and only larger structures remain.

[13]:

tracker.recording_start()

[14]:

# --------- Block 4: gaussian_sigma3 ----------

feature_name = 'gaussian_sigma3'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting gaussian_sigma3 from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_18_1.png

[15]:

tracker.recording_stop()

Edges¶

The Sobel magnitude [6] is large where brightness changes quickly, so boundaries appear bright.

[16]:

tracker.recording_start()

[17]:

# --------- Block 5: sobel_magnitude ----------

feature_name = 'sobel_magnitude'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting sobel_magnitude from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_22_1.png

[18]:

tracker.recording_stop()

Thin, line-like structures¶

The Frangi filter [24] highlights ridge- or tube-shaped structures such as thin lines. With default settings it looks for dark ridges.

[19]:

tracker.recording_start()

[20]:

# --------- Block 6: frangi_response ----------

feature_name = 'frangi_response'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting frangi_response from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_26_1.png

[21]:

tracker.recording_stop()

Oriented patterns¶

A Gabor filter [25] responds to stripe patterns of a given spacing and direction, here 45°.

[22]:

tracker.recording_start()

[23]:

# --------- Block 7: gabor_theta45 ----------

feature_name = 'gabor_theta45'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting gabor_theta45 from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_30_1.png

[24]:

tracker.recording_stop()

Texture¶

Local entropy measures how varied the brightness is within a disk of radius 5 pixels: low in smooth areas, high in textured ones.

[25]:

tracker.recording_start()

[26]:

# --------- Block 8: entropy_disk5 ----------

feature_name = 'entropy_disk5'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting entropy_disk5 from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_34_1.png

[27]:

tracker.recording_stop()

Directional order¶

Structure-tensor coherence [52] is near 1 where local texture runs in one direction (e.g. along an edge) and near 0 where it has none.

[28]:

tracker.recording_start()

[29]:

# --------- Block 9: tensor_coherence ----------

feature_name = 'tensor_coherence'
print(f'Extracting {feature_name} from all images...')

feature_results[feature_name] = {}
for path_name, img_spice in spice_dict.items():
    fmap_spice = img_spice.features.extract(feature_names=feature_name)
    feature_results[feature_name][path_name] = fmap_spice.data
    rows.append({
        'image': path_name, 'feature': feature_name,
        'mean': fmap_spice.data[..., 0].mean(),
    })

fig, axes = plt.subplots(1, len(spice_dict), figsize=(4 * len(spice_dict), 4))
if len(spice_dict) == 1:
    axes = [axes]
for i, path_name in enumerate(spice_dict.keys()):
    feature_data = feature_results[feature_name][path_name][..., 0]
    axes[i].imshow(feature_data, cmap='magma')
    axes[i].set_title(path_name[:15])
    axes[i].axis('off')
plt.suptitle(feature_name)
plt.tight_layout()
plt.show()

Extracting tensor_coherence from all images...

../../_images/notebooks_feature_extraction_feature_extract_basics_38_1.png

[30]:

tracker.recording_stop()

Compare the images¶

The table lists the mean of each feature map per image. Means allow a quick comparison but hide where a feature occurs.

[31]:

tracker.recording_start()

[32]:

# --------- Block 10: Summary table ----------

summary_df = pd.DataFrame(rows)
summary_pivot = summary_df.pivot(index='image', columns='feature', values='mean')
print(f'Feature Summary Statistics ({len(feature_results)} features extracted)')
print('=' * 60)
summary_pivot

Feature Summary Statistics (7 features extracted)
============================================================

[32]:

feature	entropy_disk5	frangi_response	gabor_theta45	gaussian_sigma1	gaussian_sigma3	sobel_magnitude	tensor_coherence
image
1b935635dd.png	4.446555	0.012061	0.000411	0.492312	0.492282	0.035964	0.624449
1ff313314c.png	4.815572	0.041037	0.000442	0.529878	0.529876	0.048803	0.523054
2fb8fe344c.png	4.344232	0.004808	0.000372	0.445533	0.445683	0.032921	0.465044

[33]:

tracker.recording_stop()

Review and save¶

The shared history holds one call per image and feature (3 × 7 = 21). We then save the data with its metadata and history script.

[34]:

tracker.recording_start()

[35]:

# --------- Block 11: Inspect unified human-readable history ----------

# All images loaded in Block 1 share the same DatasetHistory instance, so
# spice.history.to_lines() reconstructs every feature extraction performed
# across every image and every feature-map loop into ONE unified, readable
# pipeline: a single load_data(...) call followed by
# (n_images x n_feature_blocks) feature_spice = spice.features.extract(...)
# calls, each using the actual runtime feature_names value.
readable_lines = spice.history.to_lines()
n_extract_calls = sum(
    1 for line in readable_lines
    if line.startswith('feature_spice = spice.features.extract(')
)

expected_calls = len(spice_dict) * len(feature_results)
print(
    f'Unified history for spice: {len(spice_dict)} images x '
    f'{len(feature_results)} feature blocks = {n_extract_calls} extraction calls '
    f'(expected {expected_calls})'
)
assert n_extract_calls == expected_calls

print()
print('Reconstructed feature-extraction pipeline for this SPICEData object:')
print('\n'.join(readable_lines))

Unified history for spice: 3 images x 7 feature blocks = 21 extraction calls (expected 21)

Reconstructed feature-extraction pipeline for this SPICEData object:
import juspice
spice = juspice.io.load_data('/Users/amir/GIT_repositories/juspice_pre_release/Sample_data/em/1b935635dd.png')
feature_spice = spice.features.extract(feature_names=['gaussian_sigma1'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gaussian_sigma1'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gaussian_sigma1'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gaussian_sigma3'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gaussian_sigma3'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gaussian_sigma3'], n_features=1)
feature_spice = spice.features.extract(feature_names=['sobel_magnitude'], n_features=1)
feature_spice = spice.features.extract(feature_names=['sobel_magnitude'], n_features=1)
feature_spice = spice.features.extract(feature_names=['sobel_magnitude'], n_features=1)
feature_spice = spice.features.extract(feature_names=['frangi_response'], n_features=1)
feature_spice = spice.features.extract(feature_names=['frangi_response'], n_features=1)
feature_spice = spice.features.extract(feature_names=['frangi_response'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gabor_theta45'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gabor_theta45'], n_features=1)
feature_spice = spice.features.extract(feature_names=['gabor_theta45'], n_features=1)
feature_spice = spice.features.extract(feature_names=['entropy_disk5'], n_features=1)
feature_spice = spice.features.extract(feature_names=['entropy_disk5'], n_features=1)
feature_spice = spice.features.extract(feature_names=['entropy_disk5'], n_features=1)
feature_spice = spice.features.extract(feature_names=['tensor_coherence'], n_features=1)
feature_spice = spice.features.extract(feature_names=['tensor_coherence'], n_features=1)
feature_spice = spice.features.extract(feature_names=['tensor_coherence'], n_features=1)

[36]:

tracker.recording_stop()

[37]:

# --------- Block 12: Save ----------
save_data(spice)

stem = 'feature_extract_basics'
print(f'Data:           {os.path.join(notebook_dir, stem + ".npy")}')
print(f'JSON sidecar:   {os.path.join(notebook_dir, stem + ".json")}')
print(f'History script: {os.path.join(notebook_dir, stem + "_history.py")}')

Data:           /Users/amir/GIT_repositories/juspice_pre_release/notebooks/feature_extraction/feature_extract_basics.npy
JSON sidecar:   /Users/amir/GIT_repositories/juspice_pre_release/notebooks/feature_extraction/feature_extract_basics.json
History script: /Users/amir/GIT_repositories/juspice_pre_release/notebooks/feature_extraction/feature_extract_basics_history.py