Is there an efficient implementation of the fractional-order laplacian

Question

The fractional order laplacian is defined by a continuity argument to generalise the more familiar 2D or 3D laplacian, abla^2, to non-integer values of differentiation, i.e. away from the 2nd-order derivative.

This is of use to the image analysis community because the p-th order derivative is nonlocal if p is not an integer, and fractional-order derivatives have found some use in image analysis problems with the presence of noise: i.e. they are marginally more robust to single bright voxel problems.

I have written a quick-and-dirty implementation of this for my problems in Matlab, and I know that it is undoubtedly possible to do far better than I have done. Are there any efficient algorithms for the computation of the fractional order laplacian, in either Python or Matlab? My terrible code is below – I'm deeply aware there is probably an applied mathematician out there somewhere cursing under his or her breath and painfully more aware of the numerical aspects of this than I am:

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) .^ J;
    R  = toeplitz( data(:)' );
    T  = meshgrid( (gamma(order+1)/(1^order)) * signs ./ (G1.*G2) );
    out  = reshape(sum( R .* M .* T, 2 ), size(data)); 
end

This gives a similar answer to 4*del2 if order=2, without the artefact at the edges and it is (perhaps understandably) far slower. I'd like to know if there is a numerically stable, efficient implementation (ideally to >2D) -- I haven't been able to find any obvious ones (unlike, say, the "ordinary" laplacian).

Some example images of this are below.

$Example fractional order laplacian of order 0.75$

A brief example of why this is perhaps useful appears well with noise added:

Is there an efficient implementation of the fractional-order laplacian

0 Answers