Line Integral Convolution

lic.m

function licimg = lic(field)

[height, width, ~] = size(field);

noise = rand(height, width);

% lic length
L = 20.0;

licimg = zeros(height, width);
parfor i=1:height
    i
    for j=1:width
        licimg(i, j) = computeLIC(field, noise, i, j, L);
    end
end
end

function v = sampleField(field, x, y, w, h)
    lx = max(floor(x), 1);
    rx = min(lx + 1, w);
    ty = max(floor(y), 1);
    dy = min(ty + 1, h);
    
    fx = x - lx;
    fy = y - ty;
    
    v = field(ty, lx, :) * (1-fx) * (1-fy) ...
      + field(ty, rx, :) * fx * (1-fy) ...
      + field(dy, lx, :) * (1-fx) * fy ...
      + field(dy, rx, :) * fx * fy;
   
    [l, m, n] = size(v);
    v = reshape(v, l*m*n, 1);
end

function [x, y] = RK(field, x, y, w, h, s, factor)
  v = sampleField(field, x, y, w, h);
  k1 = v * s;

  v = sampleField(field, x + k1(1)*0.5, y + k1(2)*0.5, w, h);
  k2 = v * s;

  v = sampleField(field, x + k2(1)*0.5, y + k2(2)*0.5, w, h);
  k3 = v * s;

  v = sampleField(field, x + k3(1), y + k3(2), w, h);
  k4 = v * s;

  p = [x; y];
  p = p + (k1/6 + k2/3 + k3/3 + k4/6) * factor;
  x = p(1);
  y = p(2);
end

function flag = inbound(x, y, w, h)
    flag = (x>0) && (x<=w) && (y>0) && (y<=h);
end

function I = computeLIC(field, noise, i, j, L)
    [h, w, ~] = size(field);
    I = 0.0;
    weights = 0.0;
    
    % forward stream line
    x = j + 0.5;
    y = i + 0.5;
   
    s = 0.0;
    steps = 0;
    while s < L
        steps = steps + 1;
        if steps > 128 || ~inbound(x, y, w, h)
            break;
        end
        
        s = s + 1.0;
        wi = 1.0;
        weights = weights + wi;
        I = I + sampleField(noise, x, y, w, h)*wi;

        [x, y] = RK(field, x, y, w, h, 0.5, 1);
    end
    
    % backward stream line
    x = j + 0.5;
    y = i + 0.5;
   
    s = 0.0;
    steps = 0;
    while s < L
        steps = steps + 1;
        if steps > 128 || ~inbound(x, y, w, h)
            break;
        end
        
        s = s + 1.0;
        wi = 1.0;
        weights = weights + wi;
        I = I + sampleField(noise, x, y, w, h)*wi;

        [x, y] = RK(field, x, y, w, h, 0.5, -1);
    end
    
    I = I / weights;
end

%% test lic
% w = 256;
% h = 256;
% field = zeros(h, w, 2);
% for i=1:h
% y = i - h/2;
% for j=1:w
% x = j - w/2;
% L =sqrt(x_x+y_y);
% if L == 0
% field(i, j, :) = [0, 0];
% else
% field(i, j, 1) = -y/L;
% field(i, j, 2) = x/L;
% end
% end
% end

% [X,Y] = meshgrid(-2:4.0/127:2);
% Z = X._exp(-X.^2 - Y.^2);
% [Dx, Dy] = gradient(Z,.2,.2);
% Gmag = sqrt(Dx._Dx + Dy.*Dy);
% [h, w] = size(Dx);
% field = zeros(h, w, 2);
% field(:, :, 1) = Dx./Gmag;
% field(:, :, 2) = Dy./Gmag;

I = imread('bake.jpg');
I = rgb2gray(I);
[Gmag, ~] = imgradient(I);
[Gx, Gy] = imgradientxy(I);
[w, h] = size(I);
field = zeros(h, w, 2);
field(:, :, 1) = -Gy./Gmag;
field(:, :, 2) = Gx./Gmag;
quiver(field(:,:,1), field(:,:,2));

licimg = lic(field);
imagesc(licimg);colormap(gray);