Image Inpainting Using Fractional Calculus: A MATLAB Simulation Example

This is an example of using Fractional Calculus for image inpainting in MATLAB simulation. We will use a classical model based on Fractional Partial Differential Equations (FPDE).

1. Introduction to the Principles

Traditional image inpainting models based on integer-order Partial Differential Equations (PDEs), such as the TV model, have succeeded in smoothing noise and preserving edges, but can sometimes produce the “staircase effect.” Fractional Calculus extends the order of calculus from integers to fractions, inheriting the properties of integer-order calculus while possessing non-locality and weak singularity.

In image processing, fractional-order differentiation:

• Enhances Texture: Fractional-order differentiation has a stronger enhancement effect on the details and texture information of high-frequency signals in images.
• Preserves Edges: At the same time, it has a good ability to preserve mid to high-frequency information such as edges.
• Smooths Flat Areas: The processing effect on smooth areas of the image (low-frequency information) is similar to that of integer-order differentiation.

Therefore, the inpainting model based on Fractional Calculus can better preserve the texture details of the image while removing noise, alleviating the staircase effect.

A commonly used fractional-order image inpainting model is the Fractional Total Variation (FTV) model. The core idea is to minimize a fractional-order total variation energy functional.

2. Simulation Steps and MATLAB Code

We will implement a discretization method based on the Grunwald-Letnikov (G-L) definition of fractional-order differentiation. The G-L definition is very suitable for digital implementation.

Step 1: Define the Fractional Order Differentiation Mask (Kernel)First, we need a fractional-order differentiation operator. Here we create an approximate mask based on the G-L definition.

function mask = fractional_mask(alpha, n)
    % Generate fractional-order differentiation mask (1D)
    % alpha: fractional order (0 < alpha < 2)
    % n: mask length
    mask = zeros(1, n);
    mask(1) = 1; % g_0 = 1
    for k = 1:n-1
        % Calculate coefficients using the Gamma function, or use recursive approximation
        mask(k+1) = (-1)^k * gamma(alpha+1) / (gamma(k+1) * gamma(alpha - k + 1));
        % Recursive formula (more efficient, avoids Gamma function)
        % mask(k+1) = (1 - (alpha+1)/k) * mask(k);
    end
    % Normalization (optional, depends on the model)
    % mask = mask / sum(abs(mask));
end

Step 2: Core of the Image Inpainting Algorithm (FPDE Diffusion)We will use an explicit iterative diffusion process. The fractional gradient norm is part of the diffusion coefficient.

function restored_image = fractional_inpainting(damaged_image, mask, alpha, lambda, dt, iterations)
    % damaged_image: damaged image (grayscale, double type), damaged area pixel values should be NaN or 0, and identified with a binary mask.
    % mask: binary matrix, 1 indicates damaged area, 0 indicates intact area
    % alpha: fractional order
    % lambda: fidelity term weight
    % dt: iteration time step
    % iterations: number of iterations
    
    % Initialization
    u = damaged_image;
    [rows, cols] = size(u);
    
    % Create fractional-order differentiation mask
    n = 5; % mask length, can be adjusted as needed
    GL_mask = fractional_mask(alpha, n);
    
    % Iteration process
    for iter = 1:iterations
        u_old = u;
        
        % Calculate fractional-order gradient (using convolution)
        % Note: here we need to calculate the fractional partial derivatives in x and y directions
        % For simplicity, we use an approximate isotropic operator
        % A more rigorous approach is to calculate D_x^alpha and D_y^alpha separately
        
        % Use imfilter for convolution, handle boundaries ('symmetric')
        grad_alpha_x = imfilter(u, GL_mask,  'symmetric');
        grad_alpha_y = imfilter(u, GL_mask', 'symmetric');
        
        % Calculate fractional-order gradient norm
        grad_alpha_mag = sqrt(grad_alpha_x.^2 + grad_alpha_y.^2 + 1e-10); % avoid division by zero
        
        % Construct diffusion coefficient c(g) = 1 / (1 + (g^2 / k^2)) 
        k = 0.1; % contrast parameter
        c = 1 ./ (1 + (grad_alpha_mag / k).^2);
        
        % Calculate diffusion term div( c * grad_alpha(u) )
        % This is an approximation, strict divergence calculation requires the adjoint form of the fractional operator
        term_x = imfilter(c .* grad_alpha_x, GL_mask,  'symmetric');
        term_y = imfilter(c .* grad_alpha_y, GL_mask', 'symmetric');
        diffusion_term = term_x + term_y;
        
        % Update equation: u^{n+1} = u^n + dt * [ diffusion_term + lambda * mask .* (u0 - u^n) ]
        % The fidelity term only acts on the damaged area (mask==1), forcing the repair result to approach the known part of the original damaged image
        fidelity_term = mask .* (damaged_image - u); % Note: damaged_image has original values in intact areas, damaged areas are 0 or NaN
        
        u = u + dt * (diffusion_term + lambda * fidelity_term);
        
        % Display iteration progress
        if mod(iter, 50) == 0
            fprintf('Iteration: %d / %d\n', iter, iterations);
            diff = norm(u - u_old, 'fro') / norm(u, 'fro');
            fprintf('Change: %.6f\n', diff);
        end
    end
    
    restored_image = u;
end

Step 3: Main Program (Load Image, Create Damaged Area, Perform Inpainting)

% Clear workspace
clear; close all; clc;

% 1. Load image and convert to grayscale double
originalImage = im2double(imread('lena.png'));
if size(originalImage, 3) == 3
    originalImage = rgb2gray(originalImage);
end

% 2. Simulate damaged area (e.g., randomly lost blocks or scratches)
damageMask = false(size(originalImage));
% Example 1: Randomly lost pixels (20%)
% damageMask(rand(size(originalImage)) < 0.2) = true;

% Example 2: Scratch
[rows, cols] = size(originalImage);
damageMask(rows/2-10:rows/2+10, :) = true; % A horizontal scratch
% damageMask(:, cols/2-10:cols/2+10) = true; % A vertical scratch

damagedImage = originalImage;
damagedImage(damageMask) = 0; % Set damaged area pixels to 0 (or NaN)

% Display original and damaged images
figure;
subplot(1, 3, 1); imshow(originalImage); title('Original Image');
subplot(1, 3, 2); imshow(damageMask);    title('Damage Mask');
subplot(1, 3, 3); imshow(damagedImage);  title('Damaged Image');

% 3. Set fractional-order inpainting parameters
alpha = 1.2;     % Fractional order (try 0.5, 1.0, 1.5, 1.8, etc.)
lambda = 0.1;    % Fidelity term weight
dt = 0.01;       % Time step (must be small to ensure stability)
iterations = 300; % Number of iterations

% 4. Perform image inpainting
restoredResult = fractional_inpainting(damagedImage, damageMask, alpha, lambda, dt, iterations);

% 5. Display restoration results
figure;
imshow(restoredResult);
title(['Restored Image (α = ', num2str(alpha), ')']);

% 6. Calculate evaluation metrics (PSNR, SSIM)
psnrVal = psnr(restoredResult, originalImage);
ssimVal = ssim(restoredResult, originalImage);
fprintf('PSNR: %.4f dB\n', psnrVal);
fprintf('SSIM: %.4f\n', ssimVal);

3. Considerations and Advanced Topics

1. Performance and Optimization: The above code uses explicit iteration, which is less efficient. For practical applications, consider using Fast Fourier Transform (FFT) or Alternating Direction Implicit (ADI) methods for efficiency.
2. Stability: The time step must be sufficiently small to ensure the stability of the explicit format.
3. Parameter Selection: The fractional order is the most critical parameter. It degenerates to the traditional integer-order TV model. Typically, values between 1 and 2 yield better texture preservation. Parameters such as lambda, dt, and iterations need to be adjusted experimentally.
4. More Accurate Models: The divergence calculation above is a simplification. More rigorous mathematical derivation requires the adjoint form of the fractional operator, which is more complex to implement.
5. Color Images: For color images, processing is usually done separately for each channel (R, G, B).
6. Mask Generation: The recursive formula in the <span>fractional_mask</span> function is more efficient than using the <span>gamma</span> function and avoids large factorial calculations.

4. Expected Results

By running the above code (it is recommended to use standard test images such as <span>lena.png</span> or <span>cameraman.tif</span>), you will observe:

• Scratches or lost blocks in the damaged image are gradually repaired.
• Compared to integer-order models, fractional-order models (such as α) restore areas with richer texture, resulting in a more natural appearance and alleviating the “blocky” or “smearing” effect.
• You can quantitatively find the most suitable fractional order for the current image by comparing PSNR and SSIM metrics under different α values.

This example provides a beginner-level framework. To delve deeper, further academic papers on numerical implementations of fractional calculus and applications in image processing should be consulted.