function phm = pixelize(E, xVec, yVec, oversampling);

Nx = length(xVec) * oversampling;
deltaX = (xVec(2) - xVec(1)) / oversampling;
xLen = deltaX * Nx;
Ny = length(yVec) * oversampling;
deltaY = (yVec(2) - yVec(1)) / oversampling;
yLen = deltaY * Ny;

xVec = (0.5 : Nx - 0.5) * xLen / Nx - xLen / 2;
yVec = (0.5 : Ny - 0.5) * yLen / Ny - yLen / 2;

[xGrid, yGrid] = meshgrid(xVec, yVec);
phm = zeros(size(xGrid));

for i = 1 : size(E, 1)
    x0 = E(i, 1);           % x offset
    y0 = E(i, 2);           % y offset
     a = E(i, 3);            % radius
  dens = E(i, 4);        % density
     x = xGrid - x0;
     y = yGrid - y0;
   idx = find(x.^2 + y.^2 < a^2);
   phm(idx) = phm(idx) + dens;
end

phm = conv2(phm, ones(oversampling) / oversampling^2, 'same');
phm = phm(round(oversampling/2):oversampling:end, ...
      round(oversampling/2):oversampling:end);


