Bug: division returning zero

global L = 256; % grayscale.
global K = 0.01;
global C = (L * K) * (L * K); % constant used to avoid divide by zero.

function mu = mean(img)
    % number of lines and columns
    N = size(img)(1);
    M = size(img)(2);

    s = sum(img(:));
    mu = s / (N*M);
endfunction

function dp = standard_deviation(img)
    N = size(img)(1);
    M = size(img)(2);

    mi = mean(img);
    soma = 0.0;

    for i=1:N
        for j=1:M
            soma = soma + (img(i, j) - mi) * (img(i, j) - mi);
        end
    end

    variancia = soma / (N*M - 1);
    dp = sqrt(variancia);

endfunction


function corr = correlation(img1, img2)
    N1 = size(img1)(1);
    M1 = size(img1)(2);

    N2 = size(img2)(1);
    M2 = size(img2)(2);

    N = min(N1, N2);
    M = min(M1, M2);

    mi1 = mean(img1);
    mi2 = mean(img2);
    soma = 0.0;

    for i=1:N
        for j=1:M
            soma = soma + (img1(i, j) - mi1) * (img2(i, j) - mi2);
        end
    end

    corr = soma / (N*M - 1);

endfunction


function ec = struct_comp (img1, img2)
    global C;
    dp1 = standard_deviation(img1);
    dp2 = standard_deviation(img2);
    corr = correlation(img1, img2);
    disp('C:');
    disp(C);
    disp('dp1');
    disp(dp1);
    disp('dp2');
    disp(dp2);
    disp('corr');
    disp(corr);

    disp('(corr + C)');
    disp((corr + C));
    disp('(dp1 * dp2 + C)');
    disp(dp1 * dp2 + C);
    disp('(corr + C) / (dp1 * dp2 + C)');
    disp((corr + C) / (dp1 * dp2 + C));
    ec = (corr + C) / (dp1 * dp2 + C);
endfunction

name1 = 'tux.jpg';
name2 = 'monalisa.jpg';

img1 = int32(imread(name1));
img2 = int32(imread(name2));

e = struct_comp(img1, img2);
disp('Struct comp: ');
disp(e);