perform_lic.m

image denoising toolbox in matlab

function M = perform_lic(v, w, options)

% perform_lic - perform line integral convolution
%
%   M = perform_lic(v, w, options);
%
%   v is a vector field (should be approximately of unit norm).
%   w is the length of the convolution (in pixels)
%   M is an image of filtered noise that illustrate well the flow of v.
%
%   options.spot_size set the size of the features.
%
%   Set options.flow_correction=1 in order to fix problem
%   around singularity points of the flow.
%
% The method is described in
%   Imaging vector fields using line integral convolution
%   Brian Cabral, Leith Casey Leedom
%   Siggraph 1993
%   
%   Copyright (c) Gabriel Peyre 2007

n = size(v,1);

options.null = 0;
if isfield(options, 'M0')
    M0 = options.M0;
else
    M0 = randn(n);
    if isfield(options, 'spot_size')
        sigma = options.spot_size;
    else
        sigma = 2;
    end
    M0 = perform_blurring(M0, sigma, options);
end

if size(M0,3)>1
    % lic on each channel
    for i=1:size(M0,3)
        options.M0 = M0(:,:,i);
        M(:,:,i) = perform_lic(v, w, options);
    end
    return;
end

if isfield(options, 'niter_lic') && options.niter_lic>1
    M = options.M0;
    niter_lic = options.niter_lic;
    options.niter_lic = 1;
    for i=1:niter_lic
        options.M0 = M;
        M = perform_lic(v, w, options);
    end
    return;
end

if isfield(options, 'histogram')
    histogram = options.histogram;
else
    histogram = 'gaussian';
end

if isfield(options, 'flow_correction')
    flow_correction = options.flow_correction;
else
    flow_correction = 1;
end

if isstr(histogram)
    switch histogram
        case 'linear'
            hist = linspace(0,1, 100^2);
        case 'gaussian'
            hist = randn(100); hist = hist(:);
        otherwise
            error('Unkown kind of histograms');
    end
else
    hist = histogram;
end

if isfield(options, 'dt')
    dt = options.dt;
else
    dt = 0.5;
end

% perform integration of the vector field: Forward
T_list = 0:dt:w/2;
H = perform_vf_integration(v, dt, T_list, options );

% perform integration of the vector field: Backward
T_list(1) = [];
H1 = perform_vf_integration(-v, dt, T_list, options );
H = cat(3, H1(:,:,end:-1:1), H );
p = size(H,3);

% try to remove sampling problems
if flow_correction
    A = H(:,:,(end+1)/2);
    dX = real(H)-repmat(real(A), [1 1 p]);
    dY = imag(H)-repmat(imag(A), [1 1 p]);
    dX(dX>n/2) = dX(dX>n/2)-(n-1); dX(dX<-n/2) = dX(dX<-n/2)+(n-1);
    dY(dY>n/2) = dY(dY>n/2)-(n-1); dY(dY<-n/2) = dY(dY<-n/2)+(n-1);
    d = sqrt(dX.^2 + dY.^2);
    d = repmat( mean(mean(d)), [n n 1] ) - d;
    % threshold on the distance map deviation
    eta = 0.5;
end

% compute averaging
M = zeros(n);
[Y,X] = meshgrid(1:n,1:n);
W = zeros(n);
for i=1:size(H,3)
    A = interp2(1:n,1:n, M0, imag(H(:,:,i)), real(H(:,:,i)) );
    if flow_correction
        w = d(:,:,i)<eta;
    else
        w = ones(n);
    end
    M = M + A.*w;
    W = W+w;
end
W(W==0) = 1;
M = M./W;

if not(isempty(hist))
    M = perform_histogram_equalization(M, hist);
end