Venus SAR Exploratory Lab

Processing Modes

Click "Run Processing" to generate the image

RAW

Select a mode and click Run to process the image.

Choose a processing mode from the left sidebar.
Click the "Run Processing" button to apply the selected algorithm.
Results will appear in the main panel with explanations here.

Geological Feature Candidates +

Processing Log +

ML Pipeline

Current stage: Preprocessing

Click diagram to view larger version

About This Demo

Exploratory AI/ML pipeline for Venus Magellan SAR data.

Method: Unsupervised texture clustering in Haralick feature space.

Data: NASA PDS / USGS Astrogeology

No geological labels assigned - this is unsupervised surface-unit exploration.

Overview & Architecture

This developer guide provides copy-paste ready Python code examples for implementing the image processing and machine learning methods used in the Venus SAR Exploratory Lab.

Pipeline Architecture

SAR Image (GeoTIFF)
    │
    ├── Preprocessing
    │   ├── Normalization (0-1 range)
    │   └── Speckle reduction (median filter)
    │
    ├── Feature Extraction
    │   ├── GLCM Haralick textures
    │   ├── Edge detection (Sobel, Canny)
    │   └── Structure tensor (orientation)
    │
    └── ML Analysis
        ├── K-Means clustering
        ├── PCA visualization
        └── UMAP embedding

Required Libraries

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, filters, feature, exposure
from skimage.filters import median, sobel_h, sobel_v
from skimage.morphology import disk
from skimage.feature import graycomatrix, graycoprops, canny
from sklearn.cluster import MiniBatchKMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import umap

Data Source

The Venus Magellan C3-MDIR Global Mosaic is publicly available from NASA/USGS:

# Data URL
url = "https://planetarymaps.usgs.gov/mosaic/Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif"

# Download and load using requests + PIL/tifffile
import requests
from PIL import Image
import io

response = requests.get(url)
image = Image.open(io.BytesIO(response.content))
data = np.array(image)

# Or use tifffile for GeoTIFF
import tifffile
data = tifffile.imread('Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif')

# Normalize to 0-1 range
data = data.astype(np.float32)
data = (data - data.min()) / (data.max() - data.min())

Note: The full mosaic is ~170 MB (9388 x 18775 pixels). Extract a Region of Interest (ROI) for faster processing.

RAW Mode - Image Normalization

Basic preprocessing to normalize SAR amplitude values:

def process_raw(image: np.ndarray) -> np.ndarray:
    """
    Normalize SAR image to 0-1 range.
    
    Parameters:
        image: 2D numpy array of SAR backscatter values
    
    Returns:
        Normalized image (float32, range 0-1)
    """
    img = image.astype(np.float32)
    img_min, img_max = img.min(), img.max()
    
    if img_max - img_min > 0:
        normalized = (img - img_min) / (img_max - img_min)
    else:
        normalized = np.zeros_like(img)
    
    return normalized

# Visualization
fig, ax = plt.subplots(figsize=(10, 10))
im = ax.imshow(normalized_image, cmap='gray', vmin=0, vmax=1)
plt.colorbar(im, label='Normalized Radar Backscatter')
ax.set_title('Venus Magellan SAR - Raw')
plt.show()

DESPECKLE - Speckle Reduction

SAR images contain coherent speckle noise. Median filtering reduces speckle while preserving edges:

from skimage.filters import median
from skimage.morphology import disk

def process_despeckle(image: np.ndarray, radius: int = 3) -> np.ndarray:
    """
    Apply median filter for speckle reduction.
    
    Parameters:
        image: Normalized SAR image (0-1 range)
        radius: Radius of circular structuring element (default: 3)
    
    Returns:
        Despeckled image
    """
    despeckled = median(image, disk(radius))
    return despeckled

# Compare original vs despeckled
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
axes[0].imshow(original, cmap='gray')
axes[0].set_title('Original SAR Image')
axes[1].imshow(despeckled, cmap='gray')
axes[1].set_title(f'Despeckled (radius={radius})')
plt.tight_layout()

Tip: Larger radius = more smoothing but may blur fine details. Typical values: 2-5 pixels.

TEXTURE - GLCM Haralick Features

Gray-Level Co-occurrence Matrix (GLCM) extracts texture descriptors:

from skimage.feature import graycomatrix, graycoprops

def compute_glcm_features(block: np.ndarray, levels: int = 16) -> dict:
    """
    Compute GLCM Haralick texture features for an image block.
    
    Parameters:
        block: 2D image block (uint8, 0-255 or scaled to levels)
        levels: Number of gray levels for GLCM
    
    Returns:
        Dictionary with contrast, homogeneity, energy, correlation
    """
    # Scale to desired levels
    block_scaled = ((block - block.min()) / 
                    (block.max() - block.min() + 1e-8) * (levels-1)).astype(np.uint8)
    
    # Compute GLCM
    glcm = graycomatrix(block_scaled, 
                        distances=[1], 
                        angles=[0, np.pi/4, np.pi/2, 3*np.pi/4],
                        levels=levels, 
                        symmetric=True, 
                        normed=True)
    
    # Extract Haralick features (average over angles)
    features = {
        'contrast': graycoprops(glcm, 'contrast').mean(),
        'homogeneity': graycoprops(glcm, 'homogeneity').mean(),
        'energy': graycoprops(glcm, 'energy').mean(),
        'correlation': graycoprops(glcm, 'correlation').mean()
    }
    
    return features

def compute_texture_maps(image: np.ndarray, block_size: int = 16) -> dict:
    """
    Compute texture feature maps for entire image using sliding blocks.
    
    Parameters:
        image: Normalized SAR image
        block_size: Size of analysis blocks in pixels
    
    Returns:
        Dictionary of 2D texture maps
    """
    h, w = image.shape
    h_blocks, w_blocks = h // block_size, w // block_size
    
    # Initialize maps
    contrast_map = np.zeros((h_blocks, w_blocks))
    homogeneity_map = np.zeros((h_blocks, w_blocks))
    energy_map = np.zeros((h_blocks, w_blocks))
    correlation_map = np.zeros((h_blocks, w_blocks))
    
    img_uint8 = (image * 255).astype(np.uint8)
    
    for i in range(h_blocks):
        for j in range(w_blocks):
            block = img_uint8[i*block_size:(i+1)*block_size,
                             j*block_size:(j+1)*block_size]
            
            if block.max() > block.min():
                features = compute_glcm_features(block)
                contrast_map[i, j] = features['contrast']
                homogeneity_map[i, j] = features['homogeneity']
                energy_map[i, j] = features['energy']
                correlation_map[i, j] = features['correlation']
    
    return {
        'contrast': contrast_map,
        'homogeneity': homogeneity_map,
        'energy': energy_map,
        'correlation': correlation_map
    }

EDGES - Edge Detection

Sobel gradients and Canny edge detection for structural analysis:

from skimage.filters import sobel_h, sobel_v
from skimage.feature import canny

def process_edges(image: np.ndarray, sigma: float = 2.0) -> dict:
    """
    Compute edge detection using Sobel and Canny operators.
    
    Parameters:
        image: Normalized SAR image
        sigma: Gaussian smoothing for Canny
    
    Returns:
        Dictionary with gradients and edge maps
    """
    # Sobel gradients
    grad_x = sobel_h(image)  # Horizontal gradient
    grad_y = sobel_v(image)  # Vertical gradient
    
    # Gradient magnitude
    gradient_magnitude = np.sqrt(grad_x**2 + grad_y**2)
    
    # Canny edge detection
    edges = canny(image, sigma=sigma, low_threshold=0.05, high_threshold=0.15)
    
    return {
        'grad_x': grad_x,
        'grad_y': grad_y,
        'magnitude': gradient_magnitude,
        'edges': edges
    }

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes[0,0].imshow(image, cmap='gray')
axes[0,0].set_title('Original')
axes[0,1].imshow(result['grad_x'], cmap='RdBu')
axes[0,1].set_title('Sobel X (horizontal edges)')
axes[1,0].imshow(result['grad_y'], cmap='RdBu')
axes[1,0].set_title('Sobel Y (vertical edges)')
axes[1,1].imshow(result['edges'], cmap='gray')
axes[1,1].set_title('Canny Edges')

CLUSTERING - K-Means Surface Units

Unsupervised texture clustering with PCA visualization:

from sklearn.cluster import MiniBatchKMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

def extract_patch_features(patch: np.ndarray) -> np.ndarray:
    """
    Extract 6D feature vector from a patch.
    
    Features: [mean, std, contrast, homogeneity, energy, correlation]
    """
    # Intensity statistics
    mean_i = patch.mean()
    std_i = patch.std()
    
    # GLCM features
    patch_uint8 = (patch * 255).astype(np.uint8)
    patch_scaled = ((patch_uint8 - patch_uint8.min()) / 
                    (patch_uint8.max() - patch_uint8.min() + 1e-8) * 15).astype(np.uint8)
    
    glcm = graycomatrix(patch_scaled, distances=[1], angles=[0],
                        levels=16, symmetric=True, normed=True)
    
    contrast = graycoprops(glcm, 'contrast')[0, 0]
    homogeneity = graycoprops(glcm, 'homogeneity')[0, 0]
    energy = graycoprops(glcm, 'energy')[0, 0]
    correlation = graycoprops(glcm, 'correlation')[0, 0]
    
    return np.array([mean_i, std_i, contrast, homogeneity, energy, correlation])

def cluster_surface_units(image: np.ndarray, 
                          patch_size: int = 32,
                          step: int = 16,
                          n_clusters: int = 5) -> dict:
    """
    Perform K-Means clustering on texture patches.
    
    Parameters:
        image: Normalized SAR image
        patch_size: Size of patches (pixels)
        step: Stride between patches
        n_clusters: Number of clusters
    
    Returns:
        Dictionary with cluster map, labels, and PCA coordinates
    """
    h, w = image.shape
    patches = []
    positions = []
    
    # Sample patches
    for y in range(0, h - patch_size, step):
        for x in range(0, w - patch_size, step):
            patch = image[y:y+patch_size, x:x+patch_size]
            patches.append(patch)
            positions.append((y, x))
    
    # Extract features
    features = np.array([extract_patch_features(p) for p in patches])
    
    # Standardize
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features)
    
    # K-Means clustering
    kmeans = MiniBatchKMeans(n_clusters=n_clusters, random_state=42)
    labels = kmeans.fit_predict(features_scaled)
    
    # PCA for visualization
    pca = PCA(n_components=2)
    features_2d = pca.fit_transform(features_scaled)
    
    # Create cluster map
    cluster_map = np.zeros((h, w))
    for idx, (y, x) in enumerate(positions):
        cluster_map[y:y+patch_size, x:x+patch_size] = labels[idx]
    
    return {
        'cluster_map': cluster_map,
        'labels': labels,
        'features_2d': features_2d,
        'positions': positions,
        'kmeans': kmeans,
        'scaler': scaler,
        'pca': pca
    }

BOUNDARIES - Surface Unit Boundaries

Detect boundaries between texture regions using feature-space gradients:

from scipy.ndimage import gaussian_filter
from skimage.transform import resize

def compute_boundary_map(image: np.ndarray,
                        window_size: int = 32,
                        step: int = 16) -> np.ndarray:
    """
    Compute boundary probability based on texture discontinuity.
    
    Parameters:
        image: Normalized SAR image
        window_size: Analysis window size
        step: Window stride
    
    Returns:
        Boundary probability map (0-1)
    """
    h, w = image.shape
    n_y = (h - window_size) // step + 1
    n_x = (w - window_size) // step + 1
    
    # Extract features for each window
    feature_grid = np.zeros((n_y, n_x, 6))
    
    for i in range(n_y):
        for j in range(n_x):
            y, x = i * step, j * step
            window = image[y:y+window_size, x:x+window_size]
            feature_grid[i, j] = extract_patch_features(window)
    
    # Standardize features
    features_flat = feature_grid.reshape(-1, 6)
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features_flat)
    feature_grid_scaled = features_scaled.reshape(n_y, n_x, 6)
    
    # Compute boundary scores (distance to neighbors)
    boundary_scores = np.zeros((n_y, n_x))
    neighbors = [(-1, 0), (1, 0), (0, -1), (0, 1)]
    
    for i in range(n_y):
        for j in range(n_x):
            feat = feature_grid_scaled[i, j]
            distances = []
            for dy, dx in neighbors:
                ni, nj = i + dy, j + dx
                if 0 <= ni < n_y and 0 <= nj < n_x:
                    neighbor_feat = feature_grid_scaled[ni, nj]
                    dist = np.linalg.norm(feat - neighbor_feat)
                    distances.append(dist)
            if distances:
                boundary_scores[i, j] = np.mean(distances)
    
    # Normalize to 0-1
    if boundary_scores.max() > boundary_scores.min():
        boundary_scores = (boundary_scores - boundary_scores.min()) / \
                         (boundary_scores.max() - boundary_scores.min())
    
    # Upsample and smooth
    boundary_map = resize(boundary_scores, image.shape, order=1)
    boundary_map = gaussian_filter(boundary_map, sigma=1.5)
    
    return boundary_map

ORIENTATION - Structure Tensor Analysis

Compute dominant structural orientation using the structure tensor:

from skimage.filters import scharr_h, scharr_v
from scipy.ndimage import gaussian_filter
import matplotlib.colors as mcolors

def compute_orientation(image: np.ndarray, 
                        sigma: float = 2.0) -> dict:
    """
    Compute dominant orientation using structure tensor.
    
    Parameters:
        image: Normalized SAR image
        sigma: Smoothing for structure tensor
    
    Returns:
        Dictionary with orientation, coherence, and visualization
    """
    # Despeckle first
    image_smooth = median(image, disk(2))
    
    # Compute gradients (Scharr for better rotation invariance)
    Ix = scharr_h(image_smooth)
    Iy = scharr_v(image_smooth)
    
    # Structure tensor components with smoothing
    Jxx = gaussian_filter(Ix * Ix, sigma=sigma)
    Jyy = gaussian_filter(Iy * Iy, sigma=sigma)
    Jxy = gaussian_filter(Ix * Iy, sigma=sigma)
    
    # Orientation angle
    theta = 0.5 * np.arctan2(2 * Jxy, Jxx - Jyy)
    
    # Eigenvalues for coherence
    sqrt_term = np.sqrt((Jxx - Jyy)**2 + 4 * Jxy**2)
    lambda1 = 0.5 * (Jxx + Jyy + sqrt_term)
    lambda2 = 0.5 * (Jxx + Jyy - sqrt_term)
    
    # Coherence (anisotropy measure)
    epsilon = 1e-10
    coherence = (lambda1 - lambda2) / (lambda1 + lambda2 + epsilon)
    
    # Mask low-energy regions
    energy = lambda1 + lambda2
    threshold = np.percentile(energy, 10)
    mask = energy > threshold
    coherence_masked = coherence * mask
    
    # HSV visualization: Hue=orientation, Value=coherence
    hue = (theta + np.pi/2) / np.pi  # Map to [0, 1]
    sat = np.ones_like(hue)
    val = coherence_masked
    hsv = np.stack([hue, sat, val], axis=-1)
    rgb = mcolors.hsv_to_rgb(hsv)
    
    return {
        'theta': theta,
        'coherence': coherence_masked,
        'rgb': rgb
    }

MORPHOSPACE - UMAP 3D Embedding

Create a 3D morphospace visualization using UMAP dimensionality reduction:

import umap
from skimage.feature import canny

def extract_8d_features(patch: np.ndarray) -> np.ndarray:
    """
    Extract 8D feature vector for morphospace analysis.
    
    Features: [mean, std, contrast, homogeneity, energy, entropy,
               edge_density, orientation_coherence]
    """
    # Basic statistics
    mean_i = patch.mean()
    std_i = patch.std()
    
    # GLCM features
    patch_uint8 = (patch * 255).astype(np.uint8)
    patch_scaled = ((patch_uint8 - patch_uint8.min()) / 
                    (patch_uint8.max() - patch_uint8.min() + 1e-8) * 15).astype(np.uint8)
    
    glcm = graycomatrix(patch_scaled, distances=[1], angles=[0],
                        levels=16, symmetric=True, normed=True)
    
    contrast = graycoprops(glcm, 'contrast')[0, 0]
    homogeneity = graycoprops(glcm, 'homogeneity')[0, 0]
    energy = graycoprops(glcm, 'energy')[0, 0]
    
    # Entropy
    glcm_flat = glcm[:, :, 0, 0].flatten()
    glcm_flat = glcm_flat[glcm_flat > 0]
    entropy = -np.sum(glcm_flat * np.log2(glcm_flat + 1e-10))
    
    # Edge density
    edges = canny(patch, sigma=1.5)
    edge_density = edges.mean()
    
    # Orientation coherence
    Ix = sobel_h(patch)
    Iy = sobel_v(patch)
    Jxx = gaussian_filter(Ix * Ix, sigma=1.5)
    Jyy = gaussian_filter(Iy * Iy, sigma=1.5)
    Jxy = gaussian_filter(Ix * Iy, sigma=1.5)
    sqrt_term = np.sqrt((Jxx - Jyy)**2 + 4 * Jxy**2)
    lambda1 = 0.5 * (Jxx + Jyy + sqrt_term)
    lambda2 = 0.5 * (Jxx + Jyy - sqrt_term)
    coherence = ((lambda1 - lambda2) / (lambda1 + lambda2 + 1e-10)).mean()
    
    return np.array([mean_i, std_i, contrast, homogeneity, 
                     energy, entropy, edge_density, coherence])

def create_morphospace(image: np.ndarray,
                       patch_size: int = 64,
                       step: int = 32,
                       n_clusters: int = 5) -> dict:
    """
    Create 3D morphospace using UMAP embedding.
    
    Parameters:
        image: Normalized SAR image
        patch_size: Patch size for feature extraction
        step: Stride between patches
        n_clusters: Number of K-Means clusters
    
    Returns:
        Dictionary with 3D coordinates, features, and clusters
    """
    h, w = image.shape
    patches = []
    positions = []
    
    for y in range(0, h - patch_size, step):
        for x in range(0, w - patch_size, step):
            patch = image[y:y+patch_size, x:x+patch_size]
            patches.append(patch)
            positions.append((y, x))
    
    # Extract 8D features
    features = np.array([extract_8d_features(p) for p in patches])
    
    # Standardize
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features)
    
    # UMAP to 3D
    reducer = umap.UMAP(n_components=3, 
                        n_neighbors=15,
                        min_dist=0.1, 
                        metric='euclidean',
                        random_state=42)
    embedding_3d = reducer.fit_transform(features_scaled)
    
    # K-Means clustering
    kmeans = MiniBatchKMeans(n_clusters=n_clusters, random_state=42)
    labels = kmeans.fit_predict(features_scaled)
    
    return {
        'embedding': embedding_3d,
        'features': features,
        'labels': labels,
        'positions': positions,
        'patches': patches
    }

# Interactive 3D visualization with Plotly
import plotly.graph_objects as go

def plot_morphospace_3d(result: dict):
    fig = go.Figure(data=[go.Scatter3d(
        x=result['embedding'][:, 0],
        y=result['embedding'][:, 1],
        z=result['embedding'][:, 2],
        mode='markers',
        marker=dict(
            size=4,
            color=result['labels'],
            colorscale='Viridis',
            opacity=0.8
        ),
        text=[f'Cluster {l+1}' for l in result['labels']]
    )])
    
    fig.update_layout(
        title='Venus SAR Texture Morphospace',
        scene=dict(
            xaxis_title='UMAP 1',
            yaxis_title='UMAP 2',
            zaxis_title='UMAP 3'
        )
    )
    
    return fig

Requirements

Python dependencies for running all processing modes:

# requirements.txt
numpy>=1.21.0
scipy>=1.7.0
scikit-image>=0.19.0
scikit-learn>=1.0.0
matplotlib>=3.5.0
umap-learn>=0.5.0
plotly>=5.0.0
Pillow>=9.0.0
tifffile>=2022.0.0
requests>=2.27.0

Installation: pip install -r requirements.txt

Áttekintés és Architektúra

Ez a fejlesztői útmutató másolásra kész Python kódpéldákat tartalmaz a Venus SAR Exploratory Lab-ban használt képfeldolgozási és gépi tanulási módszerek implementálásához.

Pipeline Architektúra

SAR Kép (GeoTIFF)
    │
    ├── Előfeldolgozás
    │   ├── Normalizálás (0-1 tartomány)
    │   └── Pöttyösség-csökkentés (medián szűrő)
    │
    ├── Jellemző Kinyerés
    │   ├── GLCM Haralick textúrák
    │   ├── Éldetektálás (Sobel, Canny)
    │   └── Struktúra tenzor (orientáció)
    │
    └── ML Elemzés
        ├── K-közép klaszterezés
        ├── PCA vizualizáció
        └── UMAP beágyazás

Megjegyzés: A kódpéldák Python nyelven íródtak és univerzálisak - ugyanazok működnek minden környezetben.

Adatforrás

A Venus Magellan C3-MDIR Globális Mozaik nyilvánosan elérhető a NASA/USGS szervereiről:

# Adat URL
url = "https://planetarymaps.usgs.gov/mosaic/Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif"

# Letöltés és betöltés requests + PIL/tifffile használatával
import requests
from PIL import Image
import io

response = requests.get(url)
image = Image.open(io.BytesIO(response.content))
data = np.array(image)

# Vagy tifffile használata GeoTIFF-hez
import tifffile
data = tifffile.imread('Venus_Magellan_C3-MDIR_Global_Mosaic_2025m.tif')

# Normalizálás 0-1 tartományra
data = data.astype(np.float32)
data = (data - data.min()) / (data.max() - data.min())

Megjegyzés: A teljes mozaik ~170 MB (9388 x 18775 pixel). Érdemes kisebb ROI régiót kivágni a gyorsabb feldolgozáshoz.

NYERS mód - Kép normalizálás

Alapvető előfeldolgozás a SAR amplitúdó értékek normalizálásához:

def process_raw(image: np.ndarray) -> np.ndarray:
    """
    SAR kép normalizálása 0-1 tartományra.
    
    Paraméterek:
        image: 2D numpy tömb SAR visszaverődési értékekkel
    
    Visszatérési érték:
        Normalizált kép (float32, 0-1 tartomány)
    """
    img = image.astype(np.float32)
    img_min, img_max = img.min(), img.max()
    
    if img_max - img_min > 0:
        normalized = (img - img_min) / (img_max - img_min)
    else:
        normalized = np.zeros_like(img)
    
    return normalized

# Vizualizáció
fig, ax = plt.subplots(figsize=(10, 10))
im = ax.imshow(normalized_image, cmap='gray', vmin=0, vmax=1)
plt.colorbar(im, label='Normalizált radar visszaverődés')
ax.set_title('Venus Magellan SAR - Nyers')
plt.show()

ZAJSZŰRÉS - Pöttyösség csökkentés

A SAR képek koherens pöttyösséget (speckle) tartalmaznak. A mediánszűrő csökkenti a zajt az élek megőrzése mellett:

from skimage.filters import median
from skimage.morphology import disk

def process_despeckle(image: np.ndarray, radius: int = 3) -> np.ndarray:
    """
    Mediánszűrő alkalmazása pöttyösség csökkentésre.
    
    Paraméterek:
        image: Normalizált SAR kép (0-1 tartomány)
        radius: Körkörös struktúra elem sugara (alapértelmezett: 3)
    
    Visszatérési érték:
        Zajszűrt kép
    """
    despeckled = median(image, disk(radius))
    return despeckled

# Eredeti és zajszűrt összehasonlítása
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
axes[0].imshow(original, cmap='gray')
axes[0].set_title('Eredeti SAR kép')
axes[1].imshow(despeckled, cmap='gray')
axes[1].set_title(f'Zajszűrt (sugár={radius})')
plt.tight_layout()

Tipp: Nagyobb sugár = erősebb simítás, de elmoshatja a finom részleteket. Tipikus érték: 2-5 pixel.

TEXTÚRA - GLCM Haralick jellemzők

A Gray-Level Co-occurrence Matrix (GLCM) textúra leírókat számít ki:

from skimage.feature import graycomatrix, graycoprops

def compute_glcm_features(block: np.ndarray, levels: int = 16) -> dict:
    """
    GLCM Haralick textúra jellemzők kiszámítása egy képblokkhoz.
    
    Paraméterek:
        block: 2D képblokk (uint8, 0-255 vagy szintekre skálázott)
        levels: Szürkeségi szintek száma a GLCM-hez
    
    Visszatérési érték:
        Szótár kontraszt, homogenitás, energia, korreláció értékekkel
    """
    # Skálázás a kívánt szintekre
    block_scaled = ((block - block.min()) / 
                    (block.max() - block.min() + 1e-8) * (levels-1)).astype(np.uint8)
    
    # GLCM számítás
    glcm = graycomatrix(block_scaled, 
                        distances=[1], 
                        angles=[0, np.pi/4, np.pi/2, 3*np.pi/4],
                        levels=levels, 
                        symmetric=True, 
                        normed=True)
    
    # Haralick jellemzők kinyerése (átlag a szögeken)
    features = {
        'contrast': graycoprops(glcm, 'contrast').mean(),
        'homogeneity': graycoprops(glcm, 'homogeneity').mean(),
        'energy': graycoprops(glcm, 'energy').mean(),
        'correlation': graycoprops(glcm, 'correlation').mean()
    }
    
    return features

def compute_texture_maps(image: np.ndarray, block_size: int = 16) -> dict:
    """
    Textúra jellemző térképek számítása a teljes képre csúszó blokkokkal.
    
    Paraméterek:
        image: Normalizált SAR kép
        block_size: Elemző blokkok mérete pixelben
    
    Visszatérési érték:
        2D textúra térképek szótára
    """
    h, w = image.shape
    h_blocks, w_blocks = h // block_size, w // block_size
    
    # Térképek inicializálása
    contrast_map = np.zeros((h_blocks, w_blocks))
    homogeneity_map = np.zeros((h_blocks, w_blocks))
    energy_map = np.zeros((h_blocks, w_blocks))
    correlation_map = np.zeros((h_blocks, w_blocks))
    
    img_uint8 = (image * 255).astype(np.uint8)
    
    for i in range(h_blocks):
        for j in range(w_blocks):
            block = img_uint8[i*block_size:(i+1)*block_size,
                             j*block_size:(j+1)*block_size]
            
            if block.max() > block.min():
                features = compute_glcm_features(block)
                contrast_map[i, j] = features['contrast']
                homogeneity_map[i, j] = features['homogeneity']
                energy_map[i, j] = features['energy']
                correlation_map[i, j] = features['correlation']
    
    return {
        'contrast': contrast_map,
        'homogeneity': homogeneity_map,
        'energy': energy_map,
        'correlation': correlation_map
    }

ÉLEK - Éldetektálás

Sobel gradiensek és Canny éldetektálás struktúra elemzéshez:

from skimage.filters import sobel_h, sobel_v
from skimage.feature import canny

def process_edges(image: np.ndarray, sigma: float = 2.0) -> dict:
    """
    Éldetektálás Sobel és Canny operátorokkal.
    
    Paraméterek:
        image: Normalizált SAR kép
        sigma: Gauss simítás a Canny-hez
    
    Visszatérési érték:
        Gradiensek és éltérképek szótára
    """
    # Sobel gradiensek
    grad_x = sobel_h(image)  # Vízszintes gradiens
    grad_y = sobel_v(image)  # Függőleges gradiens
    
    # Gradiens nagyság
    gradient_magnitude = np.sqrt(grad_x**2 + grad_y**2)
    
    # Canny éldetektálás
    edges = canny(image, sigma=sigma, low_threshold=0.05, high_threshold=0.15)
    
    return {
        'grad_x': grad_x,
        'grad_y': grad_y,
        'magnitude': gradient_magnitude,
        'edges': edges
    }

# Vizualizáció
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes[0,0].imshow(image, cmap='gray')
axes[0,0].set_title('Eredeti')
axes[0,1].imshow(result['grad_x'], cmap='RdBu')
axes[0,1].set_title('Sobel X (vízszintes élek)')
axes[1,0].imshow(result['grad_y'], cmap='RdBu')
axes[1,0].set_title('Sobel Y (függőleges élek)')
axes[1,1].imshow(result['edges'], cmap='gray')
axes[1,1].set_title('Canny élek')

KLASZTEREZÉS - K-közép felszíni egységek

Felügyelet nélküli textúra klaszterezés PCA vizualizációval:

from sklearn.cluster import MiniBatchKMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

def extract_patch_features(patch: np.ndarray) -> np.ndarray:
    """
    6D jellemzővektor kinyerése egy régióból.
    
    Jellemzők: [átlag, szórás, kontraszt, homogenitás, energia, korreláció]
    """
    # Intenzitás statisztikák
    mean_i = patch.mean()
    std_i = patch.std()
    
    # GLCM jellemzők
    patch_uint8 = (patch * 255).astype(np.uint8)
    patch_scaled = ((patch_uint8 - patch_uint8.min()) / 
                    (patch_uint8.max() - patch_uint8.min() + 1e-8) * 15).astype(np.uint8)
    
    glcm = graycomatrix(patch_scaled, distances=[1], angles=[0],
                        levels=16, symmetric=True, normed=True)
    
    contrast = graycoprops(glcm, 'contrast')[0, 0]
    homogeneity = graycoprops(glcm, 'homogeneity')[0, 0]
    energy = graycoprops(glcm, 'energy')[0, 0]
    correlation = graycoprops(glcm, 'correlation')[0, 0]
    
    return np.array([mean_i, std_i, contrast, homogeneity, energy, correlation])

def cluster_surface_units(image: np.ndarray, 
                          patch_size: int = 32,
                          step: int = 16,
                          n_clusters: int = 5) -> dict:
    """
    K-közép klaszterezés végrehajtása textúra régiókon.
    
    Paraméterek:
        image: Normalizált SAR kép
        patch_size: Régiók mérete (pixel)
        step: Lépésköz a régiók között
        n_clusters: Klaszterek száma
    
    Visszatérési érték:
        Klasztertérkép, címkék és PCA koordináták szótára
    """
    h, w = image.shape
    patches = []
    positions = []
    
    # Régiók mintavételezése
    for y in range(0, h - patch_size, step):
        for x in range(0, w - patch_size, step):
            patch = image[y:y+patch_size, x:x+patch_size]
            patches.append(patch)
            positions.append((y, x))
    
    # Jellemzők kinyerése
    features = np.array([extract_patch_features(p) for p in patches])
    
    # Standardizálás
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features)
    
    # K-közép klaszterezés
    kmeans = MiniBatchKMeans(n_clusters=n_clusters, random_state=42)
    labels = kmeans.fit_predict(features_scaled)
    
    # PCA vizualizációhoz
    pca = PCA(n_components=2)
    features_2d = pca.fit_transform(features_scaled)
    
    # Klasztertérkép létrehozása
    cluster_map = np.zeros((h, w))
    for idx, (y, x) in enumerate(positions):
        cluster_map[y:y+patch_size, x:x+patch_size] = labels[idx]
    
    return {
        'cluster_map': cluster_map,
        'labels': labels,
        'features_2d': features_2d,
        'positions': positions,
        'kmeans': kmeans,
        'scaler': scaler,
        'pca': pca
    }

HATÁROK - Felszíni egység határok

Határok detektálása jellemzőtér-gradiensek alapján:

from scipy.ndimage import gaussian_filter
from skimage.transform import resize

def compute_boundary_map(image: np.ndarray,
                        window_size: int = 32,
                        step: int = 16) -> np.ndarray:
    """
    Határ-valószínűség számítása textúra diszkontinuitás alapján.
    
    Paraméterek:
        image: Normalizált SAR kép
        window_size: Elemzési ablak mérete
        step: Ablak lépésköz
    
    Visszatérési érték:
        Határ-valószínűségi térkép (0-1)
    """
    h, w = image.shape
    n_y = (h - window_size) // step + 1
    n_x = (w - window_size) // step + 1
    
    # Jellemzők kinyerése minden ablakra
    feature_grid = np.zeros((n_y, n_x, 6))
    
    for i in range(n_y):
        for j in range(n_x):
            y, x = i * step, j * step
            window = image[y:y+window_size, x:x+window_size]
            feature_grid[i, j] = extract_patch_features(window)
    
    # Jellemzők standardizálása
    features_flat = feature_grid.reshape(-1, 6)
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features_flat)
    feature_grid_scaled = features_scaled.reshape(n_y, n_x, 6)
    
    # Határ pontszámok számítása (távolság a szomszédoktól)
    boundary_scores = np.zeros((n_y, n_x))
    neighbors = [(-1, 0), (1, 0), (0, -1), (0, 1)]
    
    for i in range(n_y):
        for j in range(n_x):
            feat = feature_grid_scaled[i, j]
            distances = []
            for dy, dx in neighbors:
                ni, nj = i + dy, j + dx
                if 0 <= ni < n_y and 0 <= nj < n_x:
                    neighbor_feat = feature_grid_scaled[ni, nj]
                    dist = np.linalg.norm(feat - neighbor_feat)
                    distances.append(dist)
            if distances:
                boundary_scores[i, j] = np.mean(distances)
    
    # Normalizálás 0-1 tartományra
    if boundary_scores.max() > boundary_scores.min():
        boundary_scores = (boundary_scores - boundary_scores.min()) / \
                         (boundary_scores.max() - boundary_scores.min())
    
    # Felskálázás és simítás
    boundary_map = resize(boundary_scores, image.shape, order=1)
    boundary_map = gaussian_filter(boundary_map, sigma=1.5)
    
    return boundary_map

ORIENTÁCIÓ - Struktúra tenzor elemzés

Domináns szerkezeti orientáció számítása struktúra tenzorral:

from skimage.filters import scharr_h, scharr_v
from scipy.ndimage import gaussian_filter
import matplotlib.colors as mcolors

def compute_orientation(image: np.ndarray, 
                        sigma: float = 2.0) -> dict:
    """
    Domináns orientáció számítása struktúra tenzorral.
    
    Paraméterek:
        image: Normalizált SAR kép
        sigma: Simítás a struktúra tenzorhoz
    
    Visszatérési érték:
        Orientáció, koherencia és vizualizáció szótára
    """
    # Zajszűrés először
    image_smooth = median(image, disk(2))
    
    # Gradiensek számítása (Scharr jobb forgási invarianciához)
    Ix = scharr_h(image_smooth)
    Iy = scharr_v(image_smooth)
    
    # Struktúra tenzor komponensek simítással
    Jxx = gaussian_filter(Ix * Ix, sigma=sigma)
    Jyy = gaussian_filter(Iy * Iy, sigma=sigma)
    Jxy = gaussian_filter(Ix * Iy, sigma=sigma)
    
    # Orientációs szög
    theta = 0.5 * np.arctan2(2 * Jxy, Jxx - Jyy)
    
    # Sajátértékek a koherenciához
    sqrt_term = np.sqrt((Jxx - Jyy)**2 + 4 * Jxy**2)
    lambda1 = 0.5 * (Jxx + Jyy + sqrt_term)
    lambda2 = 0.5 * (Jxx + Jyy - sqrt_term)
    
    # Koherencia (anizotrópia mérték)
    epsilon = 1e-10
    coherence = (lambda1 - lambda2) / (lambda1 + lambda2 + epsilon)
    
    # Alacsony energia régiók maszkolása
    energy = lambda1 + lambda2
    threshold = np.percentile(energy, 10)
    mask = energy > threshold
    coherence_masked = coherence * mask
    
    # HSV vizualizáció: Árnyalat=orientáció, Érték=koherencia
    hue = (theta + np.pi/2) / np.pi  # Leképezés [0, 1]-re
    sat = np.ones_like(hue)
    val = coherence_masked
    hsv = np.stack([hue, sat, val], axis=-1)
    rgb = mcolors.hsv_to_rgb(hsv)
    
    return {
        'theta': theta,
        'coherence': coherence_masked,
        'rgb': rgb
    }

MORFOSZFÉRA - UMAP 3D beágyazás

3D morfoszféra vizualizáció UMAP dimenziócsökkentéssel:

import umap
from skimage.feature import canny

def extract_8d_features(patch: np.ndarray) -> np.ndarray:
    """
    8D jellemzővektor kinyerése morfoszféra elemzéshez.
    
    Jellemzők: [átlag, szórás, kontraszt, homogenitás, energia, entrópia,
               élsűrűség, orientációs_koherencia]
    """
    # Alapstatisztikák
    mean_i = patch.mean()
    std_i = patch.std()
    
    # GLCM jellemzők
    patch_uint8 = (patch * 255).astype(np.uint8)
    patch_scaled = ((patch_uint8 - patch_uint8.min()) / 
                    (patch_uint8.max() - patch_uint8.min() + 1e-8) * 15).astype(np.uint8)
    
    glcm = graycomatrix(patch_scaled, distances=[1], angles=[0],
                        levels=16, symmetric=True, normed=True)
    
    contrast = graycoprops(glcm, 'contrast')[0, 0]
    homogeneity = graycoprops(glcm, 'homogeneity')[0, 0]
    energy = graycoprops(glcm, 'energy')[0, 0]
    
    # Entrópia
    glcm_flat = glcm[:, :, 0, 0].flatten()
    glcm_flat = glcm_flat[glcm_flat > 0]
    entropy = -np.sum(glcm_flat * np.log2(glcm_flat + 1e-10))
    
    # Élsűrűség
    edges = canny(patch, sigma=1.5)
    edge_density = edges.mean()
    
    # Orientációs koherencia
    Ix = sobel_h(patch)
    Iy = sobel_v(patch)
    Jxx = gaussian_filter(Ix * Ix, sigma=1.5)
    Jyy = gaussian_filter(Iy * Iy, sigma=1.5)
    Jxy = gaussian_filter(Ix * Iy, sigma=1.5)
    sqrt_term = np.sqrt((Jxx - Jyy)**2 + 4 * Jxy**2)
    lambda1 = 0.5 * (Jxx + Jyy + sqrt_term)
    lambda2 = 0.5 * (Jxx + Jyy - sqrt_term)
    coherence = ((lambda1 - lambda2) / (lambda1 + lambda2 + 1e-10)).mean()
    
    return np.array([mean_i, std_i, contrast, homogeneity, 
                     energy, entropy, edge_density, coherence])

def create_morphospace(image: np.ndarray,
                       patch_size: int = 64,
                       step: int = 32,
                       n_clusters: int = 5) -> dict:
    """
    3D morfoszféra létrehozása UMAP beágyazással.
    
    Paraméterek:
        image: Normalizált SAR kép
        patch_size: Régióméret a jellemzőkinyeréshez
        step: Lépésköz a régiók között
        n_clusters: K-közép klaszterek száma
    
    Visszatérési érték:
        3D koordináták, jellemzők és klaszterek szótára
    """
    h, w = image.shape
    patches = []
    positions = []
    
    for y in range(0, h - patch_size, step):
        for x in range(0, w - patch_size, step):
            patch = image[y:y+patch_size, x:x+patch_size]
            patches.append(patch)
            positions.append((y, x))
    
    # 8D jellemzők kinyerése
    features = np.array([extract_8d_features(p) for p in patches])
    
    # Standardizálás
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features)
    
    # UMAP 3D-be
    reducer = umap.UMAP(n_components=3, 
                        n_neighbors=15,
                        min_dist=0.1, 
                        metric='euclidean',
                        random_state=42)
    embedding_3d = reducer.fit_transform(features_scaled)
    
    # K-közép klaszterezés
    kmeans = MiniBatchKMeans(n_clusters=n_clusters, random_state=42)
    labels = kmeans.fit_predict(features_scaled)
    
    return {
        'embedding': embedding_3d,
        'features': features,
        'labels': labels,
        'positions': positions,
        'patches': patches
    }

# Interaktív 3D vizualizáció Plotly-val
import plotly.graph_objects as go

def plot_morphospace_3d(result: dict):
    fig = go.Figure(data=[go.Scatter3d(
        x=result['embedding'][:, 0],
        y=result['embedding'][:, 1],
        z=result['embedding'][:, 2],
        mode='markers',
        marker=dict(
            size=4,
            color=result['labels'],
            colorscale='Viridis',
            opacity=0.8
        ),
        text=[f'Klaszter {l+1}' for l in result['labels']]
    )])
    
    fig.update_layout(
        title='Venus SAR Textúra Morfoszféra',
        scene=dict(
            xaxis_title='UMAP 1',
            yaxis_title='UMAP 2',
            zaxis_title='UMAP 3'
        )
    )
    
    return fig

Követelmények

Python csomagok az összes feldolgozási mód futtatásához:

# requirements.txt
numpy>=1.21.0
scipy>=1.7.0
scikit-image>=0.19.0
scikit-learn>=1.0.0
matplotlib>=3.5.0
umap-learn>=0.5.0
plotly>=5.0.0
Pillow>=9.0.0
requests>=2.27.0

Telepítés: pip install -r requirements.txt

Venus SAR Exploratory Lab

Error

Feature Space Visualization

Coherence Map

Rose Diagram

Error

Feature Space Visualization

Coherence Map

Rose Diagram

Select Region of Interest

Image Viewer

Venus SAR Analysis Pipeline

Venus 3D Surface Model

About - Venus SAR Exploratory Lab

Project Video Presentation

The Broader Goal of the Project

What is the Venus SAR Exploratory Lab?

Important Note on Scientific Validity

Data Source

How Should SAR Imagery Be Interpreted?

Processing Modes

RAW - Raw SAR Image

DESPECKLE - Noise Reduction

TEXTURE - Texture Analysis (GLCM / Haralick)

EDGES - Edge Detection (Sobel + Canny)

CLUSTERING - Unsupervised Texture Grouping (K-Means + PCA)

BOUNDARIES - Surface Unit Boundary Probability

ORIENTATION - Dominant Structural Orientation

Geological "Candidate" Detection - Experimental Module

Morphospace Explorer

How It Works

Interactive Features

Interactive Features

ROI Selection

Local Patch Analysis

Processing Log

Technical and Methodological Limitations

What the System Is Suitable For - And What It Is Not

Suitable for:

Not suitable for:

Contact and Contributors

Projekt videó bemutató

A projekt átfogó célja

Mi a Vénusz SAR Kutató Laboratórium?

Fontos megjegyzés a tudományos érvényességről

Adatforrás

Hogyan értelmezzük a SAR képeket?

Feldolgozási módok

NYERS - Nyers SAR kép

ZAJSZŰRÉS - Zajcsökkentés

TEXTÚRA - Textúraelemzés (GLCM / Haralick)

ÉLEK - Éldetektálás (Sobel + Canny)

KLASZTEREZÉS - Felügyelet nélküli textúra-csoportosítás (K-közép + PCA)

HATÁROK - Felszíni egység határ-valószínűség

ORIENTÁCIÓ - Domináns szerkezeti orientáció

Geológiai "jelölt" detektálás - Kísérleti modul

Morfoszféra böngésző

Hogyan működik

Interaktív funkciók

Interaktív funkciók

ROI kiválasztás

Lokális régió elemzés

Feldolgozási napló

Technikai és módszertani korlátozások

Mire alkalmas a rendszer - és mire nem

Alkalmas:

Nem alkalmas:

Kapcsolat és közreműködők

Developer Guide

Overview & Architecture

Pipeline Architecture

Required Libraries

Data Source

RAW Mode - Image Normalization

DESPECKLE - Speckle Reduction

TEXTURE - GLCM Haralick Features

EDGES - Edge Detection

CLUSTERING - K-Means Surface Units

BOUNDARIES - Surface Unit Boundaries

ORIENTATION - Structure Tensor Analysis

MORPHOSPACE - UMAP 3D Embedding

Requirements

Áttekintés és Architektúra

Pipeline Architektúra

Adatforrás