📘 DBSCAN (Density-Based Spatial Clustering of Applications with Noise)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn.datasets import make_blobs, make_moons
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score, adjusted_rand_score
import seaborn as sns

# Generate sample data
np.random.seed(42)

# Create datasets with different characteristics
# Dataset 1: Circular blobs
X_blobs, y_true_blobs = make_blobs(n_samples=300, centers=4, 
                                   n_features=2, cluster_std=0.5, 
                                   random_state=42)

# Dataset 2: Non-linear shapes (moons)
X_moons, y_true_moons = make_moons(n_samples=200, noise=0.1, 
                                   random_state=42)

# Dataset 3: Varying densities
X_varied = np.random.rand(250, 2) * 10
# Add dense regions
dense_region1 = np.random.multivariate_normal([2, 2], [[0.1, 0], [0, 0.1]], 50)
dense_region2 = np.random.multivariate_normal([7, 7], [[0.2, 0], [0, 0.2]], 30)
X_varied = np.vstack([X_varied, dense_region1, dense_region2])

datasets = [
    (X_blobs, "Circular Blobs"),
    (X_moons, "Non-linear Shapes"),
    (X_varied, "Varying Densities")
]

# Apply DBSCAN to each dataset
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for idx, (X, title) in enumerate(datasets):
    # Standardize data
    X_scaled = StandardScaler().fit_transform(X)

    # Apply DBSCAN
    dbscan = DBSCAN(eps=0.3, min_samples=5)
    cluster_labels = dbscan.fit_predict(X_scaled)

    # Number of clusters (excluding noise)
    n_clusters = len(set(cluster_labels)) - (1 if -1 in cluster_labels else 0)
    n_noise = list(cluster_labels).count(-1)

    # Plot original data
    axes[0, idx].scatter(X[:, 0], X[:, 1], c='blue', alpha=0.6)
    axes[0, idx].set_title(f'Original: {title}')
    axes[0, idx].set_xlabel('Feature 1')
    axes[0, idx].set_ylabel('Feature 2')

    # Plot clustered data
    unique_labels = set(cluster_labels)
    colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

    for k, col in zip(unique_labels, colors):
        if k == -1:
            # Noise points in black
            col = 'black'
            marker = 'x'
            label = 'Noise'
        else:
            marker = 'o'
            label = f'Cluster {k}'

        class_member_mask = (cluster_labels == k)
        xy = X[class_member_mask]
        axes[1, idx].scatter(xy[:, 0], xy[:, 1], c=[col], 
                           marker=marker, alpha=0.6, s=50, label=label)

    axes[1, idx].set_title(f'DBSCAN: {n_clusters} clusters, {n_noise} noise points')
    axes[1, idx].set_xlabel('Feature 1')
    axes[1, idx].set_ylabel('Feature 2')
    axes[1, idx].legend()

plt.tight_layout()
plt.show()

# Performance metrics example
X_sample, y_true = make_blobs(n_samples=200, centers=3, 
                              n_features=2, random_state=42)
X_sample = StandardScaler().fit_transform(X_sample)

dbscan_sample = DBSCAN(eps=0.5, min_samples=5)
y_pred = dbscan_sample.fit_predict(X_sample)

# Remove noise points for silhouette score calculation
mask = y_pred != -1
if np.sum(mask) > 1:
    silhouette = silhouette_score(X_sample[mask], y_pred[mask])
    print(f"Silhouette Score: {silhouette:.3f}")

# If we have true labels
ari = adjusted_rand_score(y_true, y_pred)
print(f"Adjusted Rand Index: {ari:.3f}")
print(f"Number of clusters found: {len(set(y_pred)) - (1 if -1 in y_pred else 0)}")
print(f"Number of noise points: {list(y_pred).count(-1)}")

from sklearn.neighbors import NearestNeighbors
from kneed import KneeLocator

def find_optimal_eps(X, min_samples=5, plot=True):
    """
    Find optimal eps parameter using k-distance graph
    """
    # Calculate distances to k-th nearest neighbor
    neighbors = NearestNeighbors(n_neighbors=min_samples)
    neighbors_fit = neighbors.fit(X)
    distances, indices = neighbors_fit.kneighbors(X)

    # Sort distances
    distances = np.sort(distances[:, min_samples-1], axis=0)

    if plot:
        plt.figure(figsize=(10, 6))
        plt.plot(range(len(distances)), distances)
        plt.xlabel('Data Points sorted by distance')
        plt.ylabel(f'{min_samples}-NN Distance')
        plt.title('K-Distance Graph for Optimal Eps Selection')
        plt.grid(True)

        # Find knee point
        kneedle = KneeLocator(range(len(distances)), distances, 
                             curve="convex", direction="increasing")
        if kneedle.knee:
            optimal_eps = distances[kneedle.knee]
            plt.axhline(y=optimal_eps, color='red', linestyle='--', 
                       label=f'Optimal eps ≈ {optimal_eps:.3f}')
            plt.legend()
            plt.show()
            return optimal_eps
        else:
            plt.show()
            return None
    else:
        kneedle = KneeLocator(range(len(distances)), distances, 
                             curve="convex", direction="increasing")
        return distances[kneedle.knee] if kneedle.knee else None

# Parameter sensitivity analysis
def analyze_parameter_sensitivity(X, eps_range, min_samples_range):
    """
    Analyze how different parameter combinations affect clustering
    """
    results = []

    for eps in eps_range:
        for min_samples in min_samples_range:
            dbscan = DBSCAN(eps=eps, min_samples=min_samples)
            labels = dbscan.fit_predict(X)

            n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
            n_noise = list(labels).count(-1)

            results.append({
                'eps': eps,
                'min_samples': min_samples,
                'n_clusters': n_clusters,
                'n_noise': n_noise,
                'noise_ratio': n_noise / len(X)
            })

    return results

# Example usage
X_analysis, _ = make_blobs(n_samples=300, centers=4, random_state=42)
X_analysis = StandardScaler().fit_transform(X_analysis)

# Find optimal eps
optimal_eps = find_optimal_eps(X_analysis, min_samples=5)

# Parameter sensitivity analysis
eps_range = np.arange(0.1, 1.0, 0.1)
min_samples_range = range(3, 15, 2)

results = analyze_parameter_sensitivity(X_analysis, eps_range, min_samples_range)

# Visualize parameter sensitivity
import pandas as pd

df_results = pd.DataFrame(results)
pivot_clusters = df_results.pivot(index='min_samples', columns='eps', values='n_clusters')
pivot_noise = df_results.pivot(index='min_samples', columns='eps', values='noise_ratio')

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

sns.heatmap(pivot_clusters, annot=True, fmt='d', cmap='viridis', ax=ax1)
ax1.set_title('Number of Clusters')
ax1.set_xlabel('Eps')
ax1.set_ylabel('Min Samples')

sns.heatmap(pivot_noise, annot=True, fmt='.2f', cmap='Reds', ax=ax2)
ax2.set_title('Noise Ratio')
ax2.set_xlabel('Eps')
ax2.set_ylabel('Min Samples')

plt.tight_layout()
plt.show()

# Simulate customer data for segmentation
np.random.seed(42)

# Create synthetic customer data
n_customers = 1000

# Customer features
age = np.random.normal(35, 12, n_customers)
age = np.clip(age, 18, 80)

income = np.random.lognormal(10.5, 0.5, n_customers)
income = np.clip(income, 20000, 200000)

spending_score = np.random.beta(2, 5, n_customers) * 100

# Add some correlation
spending_score += (income / 2000) + np.random.normal(0, 5, n_customers)
spending_score = np.clip(spending_score, 0, 100)

# Create customer dataset
customer_data = np.column_stack([age, income/1000, spending_score])
feature_names = ['Age', 'Income (k$)', 'Spending Score']

# Standardize the data
scaler = StandardScaler()
customer_data_scaled = scaler.fit_transform(customer_data)

# Apply DBSCAN
dbscan_customers = DBSCAN(eps=0.5, min_samples=20)
customer_clusters = dbscan_customers.fit_predict(customer_data_scaled)

# Analyze results
n_clusters = len(set(customer_clusters)) - (1 if -1 in customer_clusters else 0)
n_noise = list(customer_clusters).count(-1)

print(f"Customer Segmentation Results:")
print(f"Number of customer segments: {n_clusters}")
print(f"Number of outlier customers: {n_noise}")
print(f"Percentage of outliers: {n_noise/len(customer_data)*100:.1f}%")

# Visualize customer segments
fig = plt.figure(figsize=(15, 5))

# 2D projections
feature_pairs = [(0, 1), (0, 2), (1, 2)]
pair_names = [('Age', 'Income'), ('Age', 'Spending'), ('Income', 'Spending')]

for i, ((f1, f2), (name1, name2)) in enumerate(zip(feature_pairs, pair_names)):
    ax = plt.subplot(1, 3, i+1)

    unique_labels = set(customer_clusters)
    colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

    for k, col in zip(unique_labels, colors):
        if k == -1:
            col = 'black'
            marker = 'x'
            label = 'Outliers'
            alpha = 0.3
        else:
            marker = 'o'
            label = f'Segment {k}'
            alpha = 0.7

        class_member_mask = (customer_clusters == k)
        xy = customer_data[class_member_mask]
        plt.scatter(xy[:, f1], xy[:, f2], c=[col], marker=marker, 
                   alpha=alpha, s=30, label=label)

    plt.xlabel(name1)
    plt.ylabel(name2)
    plt.title(f'{name1} vs {name2}')
    plt.legend()

plt.tight_layout()
plt.show()

# Segment analysis
print("\nCustomer Segment Analysis:")
for cluster_id in sorted(set(customer_clusters)):
    if cluster_id == -1:
        continue

    mask = customer_clusters == cluster_id
    segment_data = customer_data[mask]

    print(f"\nSegment {cluster_id} (n={np.sum(mask)}):")
    print(f"  Average Age: {np.mean(segment_data[:, 0]):.1f} years")
    print(f"  Average Income: ${np.mean(segment_data[:, 1]*1000):,.0f}")
    print(f"  Average Spending Score: {np.mean(segment_data[:, 2]):.1f}")

import numpy as np
from scipy.spatial.distance import pdist, squareform

class DBSCAN_FromScratch:
    """
    From-scratch implementation of DBSCAN clustering algorithm
    """

    def __init__(self, eps=0.5, min_samples=5, metric='euclidean'):
        """
        Initialize DBSCAN parameters

        Parameters:
        eps: float, maximum distance between two samples for one to be 
             considered as in the neighborhood of the other
        min_samples: int, number of samples in a neighborhood for a point
                    to be considered as a core point
        metric: str, distance metric to use
        """
        self.eps = eps
        self.min_samples = min_samples
        self.metric = metric
        self.labels_ = None
        self.core_sample_indices_ = None

    def _get_neighbors(self, X, point_idx):
        """
        Find all neighbors within eps distance of a point

        Parameters:
        X: array-like, shape (n_samples, n_features)
        point_idx: int, index of the point to find neighbors for

        Returns:
        neighbors: list of indices of neighboring points
        """
        neighbors = []
        point = X[point_idx]

        for i in range(len(X)):
            if i != point_idx:
                distance = np.linalg.norm(X[i] - point)
                if distance <= self.eps:
                    neighbors.append(i)

        # Include the point itself
        neighbors.append(point_idx)
        return neighbors

    def _expand_cluster(self, X, point_idx, neighbors, cluster_id, labels, visited):
        """
        Expand cluster by adding density-reachable points

        Parameters:
        X: array-like, input data
        point_idx: int, index of core point
        neighbors: list, indices of neighbors
        cluster_id: int, current cluster identifier
        labels: array, cluster labels for all points
        visited: set, set of visited points

        Returns:
        bool: True if cluster was expanded successfully
        """
        labels[point_idx] = cluster_id

        i = 0
        while i < len(neighbors):
            neighbor_idx = neighbors[i]

            if neighbor_idx not in visited:
                visited.add(neighbor_idx)
                neighbor_neighbors = self._get_neighbors(X, neighbor_idx)

                # If neighbor is also a core point, add its neighbors
                if len(neighbor_neighbors) >= self.min_samples:
                    # Add new neighbors to the list
                    for new_neighbor in neighbor_neighbors:
                        if new_neighbor not in neighbors:
                            neighbors.append(new_neighbor)

            # If neighbor is not assigned to any cluster, assign to current cluster
            if labels[neighbor_idx] == -2:  # -2 means unassigned
                labels[neighbor_idx] = cluster_id

            i += 1

        return True

    def fit_predict(self, X):
        """
        Perform DBSCAN clustering

        Parameters:
        X: array-like, shape (n_samples, n_features)

        Returns:
        labels: array, cluster labels for each point (-1 for noise)
        """
        X = np.array(X)
        n_points = len(X)

        # Initialize labels: -2 = unassigned, -1 = noise, ≥0 = cluster id
        labels = np.full(n_points, -2, dtype=int)
        visited = set()
        cluster_id = 0
        core_samples = []

        for point_idx in range(n_points):
            if point_idx in visited:
                continue

            visited.add(point_idx)

            # Find neighbors
            neighbors = self._get_neighbors(X, point_idx)

            # Check if point is a core point
            if len(neighbors) < self.min_samples:
                # Mark as noise (may change later if it becomes border point)
                labels[point_idx] = -1
            else:
                # Point is a core point
                core_samples.append(point_idx)

                # Expand cluster from this core point
                self._expand_cluster(X, point_idx, neighbors, cluster_id, 
                                   labels, visited)
                cluster_id += 1

        self.labels_ = labels
        self.core_sample_indices_ = np.array(core_samples)

        return labels

    def fit(self, X):
        """
        Fit DBSCAN clustering

        Parameters:
        X: array-like, shape (n_samples, n_features)

        Returns:
        self: object
        """
        self.fit_predict(X)
        return self

    def get_cluster_info(self):
        """
        Get information about clustering results

        Returns:
        dict: clustering information
        """
        if self.labels_ is None:
            raise ValueError("Model has not been fitted yet.")

        unique_labels = set(self.labels_)
        n_clusters = len(unique_labels) - (1 if -1 in unique_labels else 0)
        n_noise = list(self.labels_).count(-1)

        cluster_sizes = {}
        for label in unique_labels:
            if label != -1:
                cluster_sizes[f'cluster_{label}'] = list(self.labels_).count(label)

        return {
            'n_clusters': n_clusters,
            'n_noise_points': n_noise,
            'n_core_points': len(self.core_sample_indices_),
            'cluster_sizes': cluster_sizes
        }

# Example usage and comparison with sklearn
if __name__ == "__main__":
    # Generate test data
    np.random.seed(42)
    X_test, _ = make_blobs(n_samples=150, centers=3, 
                          n_features=2, cluster_std=0.8, 
                          random_state=42)

    # Standardize data
    X_test = StandardScaler().fit_transform(X_test)

    # Our implementation
    dbscan_custom = DBSCAN_FromScratch(eps=0.3, min_samples=5)
    labels_custom = dbscan_custom.fit_predict(X_test)

    # Sklearn implementation
    dbscan_sklearn = DBSCAN(eps=0.3, min_samples=5)
    labels_sklearn = dbscan_sklearn.fit_predict(X_test)

    # Compare results
    print("Comparison of implementations:")
    print(f"Custom DBSCAN - Clusters: {len(set(labels_custom)) - (1 if -1 in labels_custom else 0)}, "
          f"Noise: {list(labels_custom).count(-1)}")
    print(f"Sklearn DBSCAN - Clusters: {len(set(labels_sklearn)) - (1 if -1 in labels_sklearn else 0)}, "
          f"Noise: {list(labels_sklearn).count(-1)}")

    # Check if results are identical (may differ due to tie-breaking)
    agreement = np.mean(labels_custom == labels_sklearn)
    print(f"Agreement between implementations: {agreement:.1%}")

    # Visualize both results
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

    # Plot custom implementation results
    unique_labels = set(labels_custom)
    colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

    for k, col in zip(unique_labels, colors):
        if k == -1:
            col = 'black'
            marker = 'x'
        else:
            marker = 'o'

        class_member_mask = (labels_custom == k)
        xy = X_test[class_member_mask]
        ax1.scatter(xy[:, 0], xy[:, 1], c=[col], marker=marker, alpha=0.7, s=50)

    ax1.set_title('Custom DBSCAN Implementation')
    ax1.set_xlabel('Feature 1')
    ax1.set_ylabel('Feature 2')

    # Plot sklearn results
    unique_labels = set(labels_sklearn)
    colors = plt.cm.Spectral(np.linspace(0, 1, len(unique_labels)))

    for k, col in zip(unique_labels, colors):
        if k == -1:
            col = 'black'
            marker = 'x'
        else:
            marker = 'o'

        class_member_mask = (labels_sklearn == k)
        xy = X_test[class_member_mask]
        ax2.scatter(xy[:, 0], xy[:, 1], c=[col], marker=marker, alpha=0.7, s=50)

    ax2.set_title('Sklearn DBSCAN')
    ax2.set_xlabel('Feature 1')
    ax2.set_ylabel('Feature 2')

    plt.tight_layout()
    plt.show()

    # Show detailed cluster info
    info = dbscan_custom.get_cluster_info()
    print("\nDetailed clustering information:")
    for key, value in info.items():
        print(f"{key}: {value}")

# Simulate network traffic data for anomaly detection
np.random.seed(42)

# Generate normal network traffic patterns
n_normal = 800
normal_packet_size = np.random.normal(1500, 300, n_normal)  # Bytes
normal_frequency = np.random.exponential(0.1, n_normal)     # Packets/sec
normal_duration = np.random.gamma(2, 2, n_normal)          # Connection duration

# Generate anomalous patterns
n_anomalies = 50

# DDoS attack - high frequency, small packets
ddos_packet_size = np.random.normal(64, 10, 20)
ddos_frequency = np.random.normal(100, 20, 20)
ddos_duration = np.random.normal(1, 0.2, 20)

# Port scanning - many short connections
scan_packet_size = np.random.normal(40, 5, 15)
scan_frequency = np.random.normal(50, 10, 15)
scan_duration = np.random.normal(0.1, 0.05, 15)

# Data exfiltration - large packets, sustained
exfil_packet_size = np.random.normal(5000, 500, 15)
exfil_frequency = np.random.normal(0.5, 0.1, 15)
exfil_duration = np.random.normal(300, 50, 15)

# Combine all data
packet_sizes = np.concatenate([normal_packet_size, ddos_packet_size, 
                              scan_packet_size, exfil_packet_size])
frequencies = np.concatenate([normal_frequency, ddos_frequency, 
                             scan_frequency, exfil_frequency])
durations = np.concatenate([normal_duration, ddos_duration, 
                           scan_duration, exfil_duration])

# Create feature matrix
network_data = np.column_stack([packet_sizes, frequencies, durations])

# Standardize features
scaler = StandardScaler()
network_data_scaled = scaler.fit_transform(network_data)

# Apply DBSCAN for anomaly detection
dbscan_network = DBSCAN(eps=0.6, min_samples=10)
network_labels = dbscan_network.fit_predict(network_data_scaled)

# Analyze results
n_clusters = len(set(network_labels)) - (1 if -1 in network_labels else 0)
n_anomalies_detected = list(network_labels).count(-1)

print(f"Network Anomaly Detection Results:")
print(f"Total connections analyzed: {len(network_data)}")
print(f"Normal behavior clusters found: {n_clusters}")
print(f"Anomalies detected: {n_anomalies_detected}")
print(f"Anomaly detection rate: {n_anomalies_detected/len(network_data)*100:.1f}%")

# Visualize results
fig = plt.figure(figsize=(15, 10))

# 3D visualization
ax1 = fig.add_subplot(221, projection='3d')

unique_labels = set(network_labels)
colors = plt.cm.Set1(np.linspace(0, 1, len(unique_labels)))

for k, col in zip(unique_labels, colors):
    if k == -1:
        col = 'red'
        marker = '^'
        label = f'Anomalies (n={list(network_labels).count(k)})'
        alpha = 0.8
        size = 60
    else:
        marker = 'o'
        label = f'Normal Cluster {k} (n={list(network_labels).count(k)})'
        alpha = 0.6
        size = 30

    class_member_mask = (network_labels == k)
    data_subset = network_data[class_member_mask]

    ax1.scatter(data_subset[:, 0], data_subset[:, 1], data_subset[:, 2],
               c=[col], marker=marker, alpha=alpha, s=size, label=label)

ax1.set_xlabel('Packet Size (bytes)')
ax1.set_ylabel('Frequency (packets/sec)')
ax1.set_zlabel('Duration (seconds)')
ax1.set_title('3D Network Traffic Analysis')
ax1.legend()

# 2D projections
projections = [(0, 1, 'Packet Size', 'Frequency'),
               (0, 2, 'Packet Size', 'Duration'),
               (1, 2, 'Frequency', 'Duration')]

for i, (f1, f2, name1, name2) in enumerate(projections):
    ax = fig.add_subplot(2, 2, i+2)

    for k, col in zip(unique_labels, colors):
        if k == -1:
            col = 'red'
            marker = '^'
            alpha = 0.8
            size = 60
        else:
            marker = 'o'
            alpha = 0.6
            size = 30

        class_member_mask = (network_labels == k)
        data_subset = network_data[class_member_mask]

        if len(data_subset) > 0:
            ax.scatter(data_subset[:, f1], data_subset[:, f2],
                      c=[col], marker=marker, alpha=alpha, s=size)

    ax.set_xlabel(name1)
    ax.set_ylabel(name2)
    ax.set_title(f'{name1} vs {name2}')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Detailed anomaly analysis
anomaly_indices = np.where(network_labels == -1)[0]
normal_indices = np.where(network_labels != -1)[0]

if len(anomaly_indices) > 0:
    print(f"\nAnomaly Characteristics:")
    anomaly_data = network_data[anomaly_indices]
    normal_data = network_data[normal_indices]

    features = ['Packet Size', 'Frequency', 'Duration']

    for i, feature in enumerate(features):
        anomaly_mean = np.mean(anomaly_data[:, i])
        normal_mean = np.mean(normal_data[:, i])

        print(f"{feature}:")
        print(f"  Anomalies - Mean: {anomaly_mean:.2f}, Std: {np.std(anomaly_data[:, i]):.2f}")
        print(f"  Normal - Mean: {normal_mean:.2f}, Std: {np.std(normal_data[:, i]):.2f}")
        print(f"  Difference: {(anomaly_mean - normal_mean)/normal_mean*100:+.1f}%")

from sklearn.datasets import load_sample_image
from skimage import segmentation, color

def image_segmentation_dbscan(image_path=None, n_segments=100):
    """
    Perform image segmentation using DBSCAN on SLIC superpixels
    """
    # Load sample image (or use sklearn's sample)
    if image_path is None:
        image = load_sample_image("flower.jpg")
    else:
        from PIL import Image
        image = np.array(Image.open(image_path))

    # Resize for faster processing
    if image.shape[0] > 300:
        from skimage.transform import resize
        image = resize(image, (300, 400), anti_aliasing=True)
        image = (image * 255).astype(np.uint8)

    print(f"Image shape: {image.shape}")

    # Convert to LAB color space for better segmentation
    image_lab = color.rgb2lab(image)

    # Generate superpixels using SLIC
    segments = segmentation.slic(image, n_segments=n_segments, compactness=10, 
                                sigma=1, start_label=1)

    # Extract features for each superpixel
    n_superpixels = np.max(segments)
    features = []

    for segment_id in range(1, n_superpixels + 1):
        mask = segments == segment_id
        if np.sum(mask) == 0:
            continue

        # Color features (mean LAB values)
        l_mean = np.mean(image_lab[mask, 0])
        a_mean = np.mean(image_lab[mask, 1])
        b_mean = np.mean(image_lab[mask, 2])

        # Texture features (standard deviation)
        l_std = np.std(image_lab[mask, 0])
        a_std = np.std(image_lab[mask, 1])
        b_std = np.std(image_lab[mask, 2])

        # Spatial features (centroid)
        y_coords, x_coords = np.where(mask)
        centroid_y = np.mean(y_coords)
        centroid_x = np.mean(x_coords)

        # Size feature
        size = np.sum(mask)

        features.append([l_mean, a_mean, b_mean, l_std, a_std, b_std,
                        centroid_y/image.shape[0], centroid_x/image.shape[1], 
                        np.log(size)])

    features = np.array(features)

    # Standardize features
    scaler = StandardScaler()
    features_scaled = scaler.fit_transform(features)

    # Apply DBSCAN clustering
    dbscan_img = DBSCAN(eps=0.5, min_samples=3)
    cluster_labels = dbscan_img.fit_predict(features_scaled)

    # Create segmented image
    segmented_image = np.zeros_like(image)
    unique_labels = set(cluster_labels)
    n_clusters = len(unique_labels) - (1 if -1 in unique_labels else 0)

    # Assign colors to clusters
    colors_palette = plt.cm.Set3(np.linspace(0, 1, n_clusters + 1))

    color_map = {}
    color_idx = 0
    for label in unique_labels:
        if label == -1:
            color_map[label] = [0, 0, 0]  # Black for noise
        else:
            color_map[label] = (colors_palette[color_idx][:3] * 255).astype(int)
            color_idx += 1

    # Apply colors to segments
    for i, (segment_id, cluster_label) in enumerate(zip(range(1, n_superpixels + 1), 
                                                        cluster_labels)):
        if i >= len(cluster_labels):
            break
        mask = segments == segment_id
        segmented_image[mask] = color_map[cluster_label]

    # Visualization
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))

    # Original image
    axes[0, 0].imshow(image)
    axes[0, 0].set_title('Original Image')
    axes[0, 0].axis('off')

    # SLIC superpixels
    axes[0, 1].imshow(segmentation.mark_boundaries(image, segments))
    axes[0, 1].set_title(f'SLIC Superpixels (n={n_superpixels})')
    axes[0, 1].axis('off')

    # DBSCAN segmentation
    axes[1, 0].imshow(segmented_image)
    axes[1, 0].set_title(f'DBSCAN Segmentation (n={n_clusters} regions)')
    axes[1, 0].axis('off')

    # Combined overlay
    overlay = image.copy()
    boundaries = segmentation.find_boundaries(segments, mode='thick')
    overlay[boundaries] = [255, 0, 0]  # Red boundaries
    axes[1, 1].imshow(overlay)
    axes[1, 1].set_title('Superpixel Boundaries')
    axes[1, 1].axis('off')

    plt.tight_layout()
    plt.show()

    print(f"Segmentation Results:")
    print(f"Number of superpixels: {n_superpixels}")
    print(f"Number of regions found: {n_clusters}")
    print(f"Number of noise superpixels: {list(cluster_labels).count(-1)}")

    return segmented_image, cluster_labels, features

# Run image segmentation
try:
    segmented_img, labels, features = image_segmentation_dbscan()
except ImportError:
    print("Skipping image segmentation example - requires additional dependencies")
    print("Install with: pip install scikit-image pillow")