分数阶拉普拉斯算子的高效实现方案咨询

分数阶拉普拉斯算子（fractional order laplacian）通过连续性论证定义，用于将常见的2D或3D拉普拉斯算子∇²推广至非整数阶导数（即二阶导数之外的场景）。

该算子在图像分析领域极具实用性：当阶数p为非整数时，p阶导数具备非局部性，在含噪图像分析中表现出对单亮体素问题更强的鲁棒性，已得到实际应用。

我已在Matlab中编写了一个简易实现，但清楚存在更优方案。现寻求Python或Matlab中数值稳定、高效的分数阶拉普拉斯算子计算算法（最好支持2维以上）——目前未找到类似普通拉普拉斯算子那样的现成实现。

我的Matlab实现代码如下：

function out=fractional_del(image, order)
%Compute the fractional lapacian of order (order)
%on the input 2D, equally-sampled scalar field (image)

   partialX = zeros(size(image)); 
   partialY = zeros(size(image)); 

   %This is undoubtedly not the best way of doing this... 
   for ix =1:size(image,1)
       partialX(ix,:)=one_dee_deriv(order,image(ix,:));
   end 
   for ix =1:size(image,2)
       partialY(:,ix)=one_dee_deriv(order,image(:,ix));
   end 

   out=partialX + partialY; 

end

function out=one_dee_deriv(order, data) 
    n_points  = numel(data);
    M  = tril( ones(n_points) );
    kernel  = 0:(n_points-1);
    G1 = gamma( kernel+1 ); G2 = gamma( order-kernel+1 );
    signs  = (-1) .^ kernel; % 修正原代码笔误：J改为kernel
    R  = toeplitz( data(:)' );
    T  = meshgrid( (gamma(order+1)/(1^order)) * signs ./ (G1.*G2) );
    out  = reshape(sum( R .* M .* T, 2 ), size(data)); 
end 

现有代码的表现：当order=2时，计算结果与4*del2类似且无边缘伪影，但运行速度明显更慢。

无噪声图像示例 #

• 原始图像：
• 0.75阶分数阶拉普拉斯结果：
• 分数阶代码计算的2阶拉普拉斯结果：
• Matlab原生del2()计算结果：

含噪声图像示例 #

• 含噪图像：
• 含噪图像的0.75阶分数阶拉普拉斯结果：
• 含噪图像的del2()计算结果：

内容的提问来源于stack exchange，提问作者Landak