Not sure where this is going, but here's one way.
Instead of doing normal line-circle intersection finding, I'm going to find the intersections between the lines and the polygonal representation of the circle that's actually being drawn. For sake of generalization, the circle is an ellipse at some center.
% an ellipse
r = [5 3]; % radii [x y]
c = [1 1]; % center [x y]
npts = 15; % nominal vertex count
% an arbitrary grid
gridx = [-4.9 -4.8 -4.5 -3.1 -1.6 -1.1 0.4 0.9 2.8 4.2];
gridy = [-4.3 -2.6 -2.3 -1.3 0.86 1.3 1.8 1.9 2.5 4];
% generate symmetric ellipse
npts = ceil((npts-1)/4)*4+1;
th = linspace(0,2*pi,npts);
xell = r(1)*cos(th) + c(1);
yell = r(2)*sin(th) + c(2);
% find intersection points between grid and polygonal circle approx
bp = ceil([1 npts/4 npts/2 npts]);
gxint1 = interp1(xell(bp(1):bp(3)),yell(bp(1):bp(3)),gridx);
gxint2 = interp1(xell(bp(3):bp(4)),yell(bp(3):bp(4)),gridx);
xell = circshift(xell(2:end),bp(2)); xell = [xell xell(1)];
yell = circshift(yell(2:end),bp(2)); yell = [yell yell(1)];
gyint1 = interp1(yell(bp(1):bp(3)),xell(bp(1):bp(3)),gridy);
gyint2 = interp1(yell(bp(3):bp(4)),xell(bp(3):bp(4)),gridy);
% plot ellipse
plot(xell,yell,'r'); hold on
axis equal
% plot grid
% exterior portion plotted for clarity
for k = 1:numel(gridx)
xline(gridx(k),'b:'); % draw the whole grid
plot([1 1]*gridx(k),[gxint1(k) gxint2(k)],'b') % the inscribed portion
end
for k = 1:numel(gridy)
yline(gridy(k),'b:'); % draw the whole grid
plot([gyint1(k) gyint2(k)],[1 1]*gridy(k),'b') % the inscribed portion
end
