
function m = mom(f, R, parent, i, j)
    m = zeros(size(f));
    
    [M, N] = size(f);
    [_, s] = sort(f(:), "descend");
    
    for x = 1:M
        for y = 1:N
            m(x, y) = x ^ i * y ^ j * f(x, y);
        end
    end
    
    for p = s'
        if parent(p) == p
            continue
        end
        
        m(parent(p)) = m(parent(p)) + m(p);
    end
end
