Fractional Calculus Edge Detection MATLAB Simulation Examples

Fractional Calculus Edge Detection MATLAB Simulation Examples

Fractional calculus edge detection can better preserve the texture details of images while effectively suppressing noise. Below are several implementations of different fractional edge detection methods.

1. Basic Fractional Differential Edge Detection

% Basic Fractional Differential Edge Detection
function basic_fractional_edge_detection()
    % Clear environment
    clear; clc; close all;
    
    % Read image
    I = imread('lena.png');
    if size(I, 3) == 3
        I = rgb2gray(I);
    end
    I = im2double(I);
    
    % Add noise to test robustness
    I_noisy = imnoise(I, 'gaussian', 0, 0.01);
    
    % Different fractional orders detection
    v_values = [0.3, 0.5, 0.7, 0.9];
    
    figure('Position', [100, 100, 1200, 800]);
    
    % Original image
    subplot(3, 5, 1);
    imshow(I);
    title('Original Image');
    
    % Noisy image
    subplot(3, 5, 2);
    imshow(I_noisy);
    title('Noisy Image');
    
    % Sobel operator comparison
    subplot(3, 5, 3);
    edges_sobel = edge(I_noisy, 'sobel');
    imshow(edges_sobel);
    title('Sobel Edge Detection');
    
    % Canny operator comparison
    subplot(3, 5, 4);
    edges_canny = edge(I_noisy, 'canny');
    imshow(edges_canny);
    title('Canny Edge Detection');
    
    % Results of different fractional orders detection
    for i = 1:4
        v = v_values(i);
        
        % Construct fractional differential mask
        mask_size = 5;
        mask = construct_frac_mask(v, mask_size);
        
        % Apply fractional differential
        I_frac = imfilter(I_noisy, mask, 'replicate');
        
        % Binarize edges
        threshold = 0.15;
        edges_frac = imbinarize(abs(I_frac), threshold);
        
        % Display results
        subplot(3, 5, 5 + i);
        imshow(I_frac);
        title(['Fractional Differential v=', num2str(v)]);
        
        subplot(3, 5, 10 + i);
        imshow(edges_frac);
        title(['Binarized Edges v=', num2str(v)]);
    end
    
    % Save results
    saveas(gcf, 'fractional_edge_detection_results.png');
end

% Construct fractional differential mask
function mask = construct_frac_mask(v, mask_size)
    mask = zeros(mask_size);
    center = ceil(mask_size / 2);
    
    % Based on Grünwald-Letnikov definition
    for i = 1:mask_size
        for j = 1:mask_size
            dx = i - center;
            dy = j - center;
            
            if dx == 0 && dy == 0
                % Center point coefficient
                mask(i, j) = gamma(2 - v) / (gamma(1 - v)^2);
            else
                % Other point coefficients
                r = sqrt(dx^2 + dy^2);
                mask(i, j) = (-1)^(abs(dx) + abs(dy)) * gamma(1 - v) / ...
                            (gamma(abs(dx) - v + 1) * gamma(abs(dy) - v + 1));
            end
        end
    end
    
    % Normalize
    mask = mask / sum(abs(mask(:)));
end

2. Multi-Directional Fractional Edge Detection

% Multi-Directional Fractional Edge Detection
function multi_direction_fractional_edge()
    % Clear environment
    clear; clc; close all;
    
    % Read image
    I = im2double(rgb2gray(imread('cameraman.tif')));
    
    % Parameter settings
    v = 0.6;        % Fractional order
    mask_size = 7;  % Mask size
    directions = 8; % Number of detection directions
    
    % Construct multi-directional fractional masks
    masks = create_multi_direction_masks(v, mask_size, directions);
    
    % Multi-directional edge detection
    edge_maps = zeros([size(I), directions]);
    
    figure('Position', [100, 100, 1200, 600]);
    
    for d = 1:directions
        % Apply fractional differential
        edge_response = imfilter(I, masks{d}, 'replicate');
        edge_maps(:, :, d) = abs(edge_response);
        
        % Display single direction result
        subplot(2, 4, d);
        imshow(edge_maps(:, :, d), []);
        title(['Direction ', num2str(d*45), '°']);
        colormap(jet);
        colorbar;
    end
    
    % Fuse multi-directional results
    fused_edges = max(edge_maps, [], 3);
    
    % Binarize
    threshold = 0.12;
    binary_edges = fused_edges > threshold;
    
    % Display fusion results
    figure;
    subplot(1, 3, 1);
    imshow(I);
    title('Original Image');
    
    subplot(1, 3, 2);
    imshow(fused_edges, []);
    title('Multi-Directional Fractional Edge Response');
    colormap(jet);
    colorbar;
    
    subplot(1, 3, 3);
    imshow(binary_edges);
    title('Binarized Edge Map');
    
    % Compare with traditional methods
    compare_with_traditional(I, fused_edges);
end

% Create multi-directional fractional masks
function masks = create_multi_direction_masks(v, mask_size, num_directions)
    masks = cell(num_directions, 1);
    center = ceil(mask_size / 2);
    
    for d = 1:num_directions
        angle = (d - 1) * pi / num_directions;
        mask = zeros(mask_size);
        
        for i = 1:mask_size
            for j = 1:mask_size
                dx = i - center;
                dy = j - center;
                
                % Rotate coordinates to current direction
                x_rot = dx * cos(angle) + dy * sin(angle);
                y_rot = -dx * sin(angle) + dy * cos(angle);
                
                if abs(x_rot) < 0.5 && abs(y_rot) < 0.5
                    mask(i, j) = gamma(2 - v) / (gamma(1 - v)^2);
                else
                    mask(i, j) = (-1)^(abs(round(x_rot)) + abs(round(y_rot))) * ...
                                gamma(1 - v) / (gamma(abs(round(x_rot)) - v + 1) * ...
                                gamma(abs(round(y_rot)) - v + 1));
                end
            end
        end
        
        masks{d} = mask / sum(abs(mask(:)));
    end
end

% Compare with traditional methods
function compare_with_traditional(I, frac_edges)
    % Traditional edge detection methods
    sobel_edges = edge(I, 'sobel');
    prewitt_edges = edge(I, 'prewitt');
    log_edges = edge(I, 'log');
    canny_edges = edge(I, 'canny');
    
    % Fractional edge binarization
    threshold = 0.15;
    frac_binary = frac_edges > threshold;
    
    figure('Position', [100, 100, 1200, 400]);
    
    subplot(2, 5, 1);
    imshow(I);
    title('Original Image');
    
    subplot(2, 5, 2);
    imshow(sobel_edges);
    title('Sobel');
    
    subplot(2, 5, 3);
    imshow(prewitt_edges);
    title('Prewitt');
    
    subplot(2, 5, 4);
    imshow(log_edges);
    title('LoG');
    
    subplot(2, 5, 5);
    imshow(canny_edges);
    title('Canny');
    
    subplot(2, 5, 6);
    imshow(frac_edges, []);
    title('Fractional Response');
    colormap(jet);
    
    subplot(2, 5, 7);
    imshow(frac_binary);
    title('Fractional Binarization');
    
    % Calculate evaluation metrics
    [~, gt_edges] = edge(I, 'canny'); % Use Canny as reference
    
    % Calculate F1 scores
    f1_scores = zeros(5, 1);
    methods = {sobel_edges, prewitt_edges, log_edges, canny_edges, frac_binary};
    method_names = {'Sobel', 'Prewitt', 'LoG', 'Canny', 'Fractional'};
    
    for i = 1:5
        f1_scores(i) = calculate_f1_score(gt_edges, methods{i});
    end
    
    % Display F1 scores
    subplot(2, 5, 8);
    bar(f1_scores);
    set(gca, 'XTickLabel', method_names);
    title('F1 Score Comparison');
    ylabel('F1 Score');
    rotateXLabels(gca, 45);
    
    disp('=== Edge Detection Performance Comparison ===');
    for i = 1:5
        fprintf('%s: F1 = %.4f\n', method_names{i}, f1_scores(i));
    end
end

% Calculate F1 score
function f1 = calculate_f1_score(gt, detected)
    tp = sum(gt(:) && detected(:));      % True positives
    fp = sum(~gt(:) && detected(:));     % False positives
    fn = sum(gt(:) && ~detected(:));     % False negatives
    
    precision = tp / (tp + fp + eps);
    recall = tp / (tp + fn + eps);
    f1 = 2 * precision * recall / (precision + recall + eps);
end

3. Adaptive Fractional Edge Detection

% Adaptive Fractional Edge Detection
function adaptive_fractional_edge()
    % Clear environment
    clear; clc; close all;
    
    % Read image
    I = im2double(rgb2gray(imread('peppers.png')));
    
    % Add noise
    I_noisy = imnoise(I, 'gaussian', 0, 0.02);
    
    % Adaptive fractional edge detection
    [adaptive_edges, v_map] = adaptive_frac_edge_detection(I_noisy);
    
    % Display results
    figure('Position', [100, 100, 1200, 400]);
    
    subplot(1, 4, 1);
    imshow(I);
    title('Original Image');
    
    subplot(1, 4, 2);
    imshow(I_noisy);
    title('Noisy Image');
    
    subplot(1, 4, 3);
    imshow(v_map, []);
    title('Adaptive Fractional Order Map');
    colormap(jet);
    colorbar;
    
    subplot(1, 4, 4);
    imshow(adaptive_edges);
    title('Adaptive Fractional Edges');
    
    % Fixed fractional order comparison
    compare_fixed_vs_adaptive(I_noisy);
end

% Adaptive fractional edge detection function
function [edges, v_map] = adaptive_frac_edge_detection(I)
    [rows, cols] = size(I);
    edges = zeros(rows, cols);
    v_map = zeros(rows, cols);
    
    % Local window parameters
    window_size = 11;
    half_win = floor(window_size / 2);
    
    % Fractional order range
    v_min = 0.3;
    v_max = 0.9;
    
    % Traverse image
    for i = half_win + 1:rows - half_win
        for j = half_win + 1:cols - half_win
            % Extract local window
            window = I(i-half_win:i+half_win, j-half_win:j+half_win);
            
            % Calculate local statistical features
            local_std = std(window(:));
            local_grad = max([abs(I(i+1,j)-I(i-1,j)), abs(I(i,j+1)-I(i,j-1))]);
            
            % Adaptive selection of fractional order
            v = v_min + (v_max - v_min) * (1 - exp(-local_std / 0.1));
            v_map(i, j) = v;
            
            % Construct fractional mask
            mask = construct_adaptive_mask(v, local_grad);
            
            % Apply fractional differential
            edge_response = sum(sum(window .* mask));
            edges(i, j) = abs(edge_response);
        end
    end
    
    % Binarize
    threshold = 0.1;
    edges = edges > threshold;
end

% Adaptive mask construction
function mask = construct_adaptive_mask(v, gradient)
    mask_size = 5;
    mask = zeros(mask_size);
    center = ceil(mask_size / 2);
    
    % Adjust mask shape based on gradient
    if gradient > 0.1
        % Strong edge region, use sharper mask
        beta = 0.8;
    else
        % Flat region, use smoother mask
        beta = 0.3;
    end
    
    for i = 1:mask_size
        for j = 1:mask_size
            dx = i - center;
            dy = j - center;
            
            if dx == 0 && dy == 0
                mask(i, j) = gamma(2 - v) / (gamma(1 - v)^2);
            else
                r = sqrt(dx^2 + dy^2);
                mask(i, j) = (-1)^(abs(dx) + abs(dy)) * gamma(1 - v) / ...
                            (gamma(abs(dx) - v + beta) * gamma(abs(dy) - v + beta));
            end
        end
    end
    
    mask = mask / sum(abs(mask(:)));
end

% Compare fixed and adaptive
function compare_fixed_vs_adaptive(I)
    % Fixed fractional order detection
    v_fixed = 0.6;
    mask_fixed = construct_frac_mask(v_fixed, 5);
    edges_fixed = imfilter(I, mask_fixed, 'replicate');
    edges_fixed_binary = abs(edges_fixed) > 0.1;
    
    % Adaptive detection
    [edges_adaptive, v_map] = adaptive_frac_edge_detection(I);
    
    figure('Position', [100, 100, 1200, 400]);
    
    subplot(1, 4, 1);
    imshow(I);
    title('Input Image');
    
    subplot(1, 4, 2);
    imshow(abs(edges_fixed), []);
    title(['Fixed v=', num2str(v_fixed), ' Response']);
    
    subplot(1, 4, 3);
    imshow(edges_fixed_binary);
    title('Fixed Fractional Edges');
    
    subplot(1, 4, 4);
    imshow(edges_adaptive);
    title('Adaptive Fractional Edges');
    
    % Performance comparison
    [~, gt_edges] = edge(I, 'canny');
    f1_fixed = calculate_f1_score(gt_edges, edges_fixed_binary);
    f1_adaptive = calculate_f1_score(gt_edges, edges_adaptive);
    
    fprintf('Fixed Fractional (v=%.1f) F1: %.4f\n', v_fixed, f1_fixed);
    fprintf('Adaptive Fractional F1: %.4f\n', f1_adaptive);
end
end

4. Fractional Canny Edge Detection

% Fractional Canny Edge Detection
function fractional_canny_edge()
    % Clear environment
    clear; clc; close all;
    
    % Read image
    I = im2double(rgb2gray(imread('coins.png')));
    
    % Add noise
    I_noisy = imnoise(I, 'gaussian', 0, 0.03);
    
    % Fractional Canny edge detection
    frac_canny_edges = fractional_canny(I_noisy, 0.7);
    
    % Traditional Canny
    trad_canny_edges = edge(I_noisy, 'canny');
    
    % Display results
    figure('Position', [100, 100, 1000, 300]);
    subplot(1, 3, 1);
    imshow(I_noisy);
    title('Noisy Image');
    
    subplot(1, 3, 2);
    imshow(trad_canny_edges);
    title('Traditional Canny');
    
    subplot(1, 3, 3);
    imshow(frac_canny_edges);
    title('Fractional Canny');
    
    % Performance evaluation
    [~, gt_edges] = edge(I, 'canny');
    f1_trad = calculate_f1_score(gt_edges, trad_canny_edges);
    f1_frac = calculate_f1_score(gt_edges, frac_canny_edges);
    
    fprintf('Traditional Canny F1: %.4f\n', f1_trad);
    fprintf('Fractional Canny F1: %.4f\n', f1_frac);
end

% Fractional Canny implementation
function edges = fractional_canny(I, v)
    % Step 1: Fractional Gaussian smoothing
    sigma = 1.5;
    gaussian_size = 5;
    I_smooth = imgaussfilt(I, sigma);
    
    % Step 2: Fractional gradient calculation
    [Gx, Gy] = fractional_gradient_2d(I_smooth, v);
    
    % Step 3: Gradient magnitude and direction
    G_mag = sqrt(Gx.^2 + Gy.^2);
    G_dir = atan2(Gy, Gx);
    
    % Step 4: Non-maximum suppression
    G_suppressed = non_maximum_suppression(G_mag, G_dir);
    
    % Step 5: Double threshold detection
    low_thresh = 0.05;
    high_thresh = 0.15;
    edges = double_threshold(G_suppressed, low_thresh, high_thresh);
end

% 2D fractional gradient
function [Gx, Gy] = fractional_gradient_2d(I, v)
    [rows, cols] = size(I);
    Gx = zeros(rows, cols);
    Gy = zeros(rows, cols);
    
    % x-direction fractional differential
    for i = 1:rows
        for j = 2:cols
            Gx(i, j) = I(i, j) - v * I(i, j-1);
            for k = 2:min(5, j-1)
                Gx(i, j) = Gx(i, j) + (-1)^k * gamma(v+1) / ...
                          (gamma(k+1) * gamma(v-k+1)) * I(i, j-k);
            end
        end
    end
    
    % y-direction fractional differential
    for j = 1:cols
        for i = 2:rows
            Gy(i, j) = I(i, j) - v * I(i-1, j);
            for k = 2:min(5, i-1)
                Gy(i, j) = Gy(i, j) + (-1)^k * gamma(v+1) / ...
                          (gamma(k+1) * gamma(v-k+1)) * I(i-k, j);
            end
        end
    end
end

% Non-maximum suppression
function suppressed = non_maximum_suppression(mag, dir)
    suppressed = zeros(size(mag));
    [rows, cols] = size(mag);
    
    % Quantize direction to 4 main directions
    dir_quant = mod(round(dir * 4 / pi), 4);
    
    for i = 2:rows-1
        for j = 2:cols-1
            switch dir_quant(i, j)
                case 0  % 0° direction
                    neighbors = [mag(i, j-1), mag(i, j+1)];
                case 1  % 45° direction
                    neighbors = [mag(i-1, j+1), mag(i+1, j-1)];
                case 2  % 90° direction
                    neighbors = [mag(i-1, j), mag(i+1, j)];
                case 3  % 135° direction
                    neighbors = [mag(i-1, j-1), mag(i+1, j+1)];
            end
            
            if mag(i, j) >= max(neighbors)
                suppressed(i, j) = mag(i, j);
            end
        end
    end
end

% Double threshold detection
function edges = double_threshold(mag, low, high)
    strong_edges = mag > high;
    weak_edges = (mag >= low) && (mag <= high);
    
    % Connect weak edges
    edges = strong_edges;
    [rows, cols] = size(mag);
    
    for i = 2:rows-1
        for j = 2:cols-1
            if weak_edges(i, j) && any(any(strong_edges(i-1:i+1, j-1:j+1)))
                edges(i, j) = true;
            end
        end
    end
end

Usage Instructions

  1. Image Preparation: Ensure the test images (lena.png, cameraman.tif, peppers.png, coins.png) are in the MATLAB path

  2. Parameter Adjustment:

  • Fractional Order: Typically 0.3-0.9, smaller values yield smoother results, larger values enhance edge response
  • Mask Size: Typically 3×3, 5×5, 7×7; larger sizes increase computation but improve results
  • Threshold: Adjust based on image contrast, usually between 0.1-0.2

  • Performance Evaluation:

    • Quantitatively evaluate edge detection performance using F1 score
    • Compare with traditional methods (Sobel, Prewitt, Canny, etc.)

  • Applicable Scenarios:

    • Texture-rich images
    • Edge detection in noisy environments
    • Edge detection tasks requiring detail preservation

    These methods demonstrate the advantages of fractional calculus in edge detection, particularly in noise robustness and detail preservation.

    Related posts

    1. Channel Estimation Based on MATLAB Millimeter-Wave Hybrid Massive MIMO System GAMP Channel Estimation (Including MATLAB Source Code)
    2. Fast Channel Estimation in Transform Domain for Millimeter Wave Systems Based on MATLAB [Includes MATLAB Source Code]
    3. An Open-Source Multi-Function Tester Built with STM32: Logic Analyzer, Function Generator, Oscilloscope, SPI/I2C/USART Bridge, 8-Bit Parallel Port, PWM
    4. China Telecom’s 5G Black Technology: Multi-RAT, Modular MM, and PON Transmission
    5. The Future of c-MET: Dual Antibodies, ADCs, or New Technologies?
    6. Apple Watch Set for Major Transformation in 2026 with Doubled Sensors and Enhanced Battery Life
    7. Design of a Multifunctional Electronic Scale Based on the 51 Microcontroller
    8. MATLAB Functions Supporting Variable Number of Outputs
    9. Introduction to Computing and Networking Convergence
    10. Edge Artificial Intelligence: Transforming Industrial Processes
    11. Exploration of Intelligent Highway Edge Computing Architecture and Network System
    12. The Boundaries of Operations Automation: A Comparative Analysis of Ansible and Terraform

    Leave a Comment