
function m = momc(f, R, parent, i, j)
    m = zeros(size(f));
    
    [M, N] = size(f);
    [_, s] = sort(f(:), "descend");
    
    [mx, my] = mxy(f, R, parent);
    
    for p = s'
        q = p;
        
        while true
            x = mod(q - 1, M) + 1;
            y = ceil(q / M);

            m(q) = m(q) + (x - mx(q)) ^ i * (y - my(q)) ^ j * f(q);
            
            if parent(q) == q
                break;
            end
            
            q = parent(q);
        end
    end
end
