function finelt

beta1 = 1E-2;
beta2 = 4;
lambda = 2;

% create triangles
img = imread ( 'filter2d.pgm' );
[ N, M ] = size ( img );
h = 1 / ( N - 1 );

np = N * M;

TI = zeros ( np, 1 );
TI ( : ) = img;

TX = zeros ( np, 1 );
TY = zeros ( np, 1 );

k = 1;
for i = 1 : M
    for j = 1 : N
        TX ( k ) = ( i - 1 ) * h;
        TY ( k ) = ( j - 1 ) * h;
        k = k + 1;
    end
end

nt = 2 * ( N - 1 ) * ( M - 1 );

TT = zeros ( nt, 3 );

k = 1;
for i = 1 : M - 1
    for j = 1 : N - 1
        t = ( i - 1 ) * N + j;
        TT ( k, : ) = [ ( t + N ) ( t ) ( t + 1 ) ] ;
        TT ( k + 1, : ) = [ ( t + 1 ) ( t + N + 1 ) ( t + N ) ] ;
        k = k + 2;
    end
end

% preallocate matrix data
nijv = 9 * nt;
LI = zeros ( nijv, 1 );
LJ = zeros ( nijv, 1 );
LV = zeros ( nijv, 1 );
MI = zeros ( nijv, 1 );
MJ = zeros ( nijv, 1 );
MV = zeros ( nijv, 1 );

% precompute some values
A = zeros ( nt, 1 ); % area of triangle
C = zeros ( nt, 3 ); % cosine of angle at node
S = zeros ( nt, 3 ); % sine of angle at node
H = zeros ( nt, 3 ); % slope of basis function
G = zeros ( nt, 3, 2 ); % gradient of basis function
for t = 1 : nt
    T = TT ( t, : );
    for s = 1 : 3
        % vector pointing from node s to next node
        dx = TX ( T ( s ) ) - TX ( T ( n ( s ) ) );
        dy = TY ( T ( s ) ) - TY ( T ( n ( s ) ) );
        E ( s, : ) = [ dx dy ]; 
        L ( s ) = norm ( E ( s, : ) );
        D ( s, : ) = E ( s, : ) * 1 / L ( s );
    end
    for s = 1 : 3
        C ( t, s ) = D ( p ( s ), : ) * D ( s, : )';
        S ( t, s ) = det ( D ( [ ( p ( s ) ) s ], : ) );
    end
    for s = 1 : 3
        H ( t, s ) = 1 / ( L ( s ) * S ( t, n ( s ) ) );
        N = E ( n ( s ), : );
        N = N / sqrt ( N ( 1 ) ^ 2 + N ( 2 ) ^ 2 );
        N = [ ( - N ( 2 ) ) ( N ( 1 ) ) ];
        N = H ( t, s ) * N;
        G ( t, s, : ) = N;
    end
    A ( t ) = 0.5 * abs ( det ( E ( 1 : 2, : ) ) );
end

% assemble mass matrix
Mloc = [ 2 1 1; 1 2 1; 1 1 2 ];
k = 1;
for t = 1 : nt
    T = TT ( t, : );
    fac = A ( t ) / 12.;
    for i = 1 : 3
        for j = 1 : 3
            val = fac * Mloc ( i, j );
            MI ( k ) = T ( i ); MJ ( k ) = T ( j ); MV ( k ) = val; k = k + 1;
        end
    end
end
MI ( k : nijv ) = [];
MJ ( k : nijv ) = [];
MV ( k : nijv ) = [];
M = sparse ( MI, MJ, MV );

% assemble stiffness matrix
O = [ 0 3 2; 3 0 1; 2 1 0 ]; % O ( i, j ) is index of other vertex for i ~= j
k = 1;
for t = 1 : nt
    T = TT ( t, : );
    fac = A ( t );
    for i = 1 : 3
        val = fac * H ( t, i ) ^ 2;
        LI ( k ) = T ( i ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        for j = i + 1 : 3 
            val = fac * H ( t, i ) * H ( t , j ) * C ( t, O ( i, j ) );
            LI ( k ) = T ( i ); LJ ( k ) = T ( j ); LV ( k ) = val; k = k + 1;
            LI ( k ) = T ( j ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        end
    end
end
LI ( k : nijv ) = [];
LJ ( k : nijv ) = [];
LV ( k : nijv ) = [];
L = sparse ( LI, LJ, LV );

% presmoothing
MTI = M * TI;
LPM = beta1 * L + M;
TIs = LPM \ MTI;

% perona-malik
k = 1;
for t = 1 : nt
    T = TT ( t, : );
    g = [ 0 0 ];
    for i = 1 : 3
        g ( 1 ) = g ( 1 ) + ( G ( t, i, 1 ) * TIs ( T ( i ) ) );
        g ( 2 ) = g ( 2 ) + ( G ( t, i, 2 ) * TIs ( T ( i ) ) );
    end
    gg = sqrt ( g ( 1 ) * g ( 1 ) + g ( 2 ) * g ( 2 ) );
    fac = A ( t ) / ( 1 + ( gg / lambda ) ^ 2 );
    for i = 1 : 3
        val = fac * H ( t, i ) ^ 2;
        LI ( k ) = T ( i ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        for j = i + 1 : 3 
            val = fac * H ( t, i ) * H ( t , j ) * C ( t, O ( i, j ) );
            LI ( k ) = T ( i ); LJ ( k ) = T ( j ); LV ( k ) = val; k = k + 1;
            LI ( k ) = T ( j ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        end
    end
end
LI ( k : nijv ) = [];
LJ ( k : nijv ) = [];
LV ( k : nijv ) = [];
L = sparse ( LI, LJ, LV );

LPM = beta2 * L + M;
Zpm = LPM \ MTI;

% visualize
colormap ( gray );
h1 = subplot ( 221 );
trisurf ( TT, TX, TY, TI );
view ( 2 );
h2 = subplot ( 222 );
imagesc ( img );
axis off;
h3 = subplot ( 223 );
img ( : ) = TIs;
imagesc ( img );
axis off;
h4 = subplot ( 224 );
img ( : ) = Zpm;
imagesc ( img );
axis off;
axis ( [ h1 h2 h3 h4 ], 'square' );

% index of next and previous vertex
function n = n ( s )
n = mod ( s, 3 ) + 1;
function p = p ( s )
p = mod ( s + 1, 3 ) + 1;