```procedure DrawDisk(Bitmap: TBitmap; CenterX, CenterY, Radius,
Feather: single);
// Draw a disk on Bitmap. Bitmap must be a 256 color (pf8bit)
// palette bitmap, and parts outside the disk will get palette
// index 0, parts inside will get palette index 255, and in the
// antialiased area (feather), the pixels will get values
// inbetween.
// ***Parameters***
// Bitmap:
//   The bitmap to draw on
// CenterX, CenterY:
//   The center of the disk (float precision). Note that [0, 0]
//   would be the center of the first pixel. To draw in the
//   exact middle of a 100x100 bitmap, use CenterX = 49.5 and
//   CenterY = 49.5
//   The radius of the drawn disk in pixels (float precision)
// Feather:
//   The feather area. Use 1 pixel for a 1-pixel antialiased
//   area. Pixel centers outside 'Radius + Feather / 2' become
//   0, pixel centers inside 'Radius - Feather/2' become 255.
//   Using a value of 0 will yield a bilevel image.
// Copyright (c) 2003 Nils Haeck M.Sc. www.simdesign.nl
var
x, y: integer;
LX, RX, LY, RY: integer;
Fact: integer;
RPF2, RMF2: single;
P: PByteArray;
SqY, SqDist: single;
sqX: array of single;
begin
// Determine some helpful values (singles)

// Determine bounds:
LX := Max(floor(CenterX - RPF2), 0);
RX := Min(ceil (CenterX + RPF2), Bitmap.Width - 1);
LY := Max(floor(CenterY - RPF2), 0);
RY := Min(ceil (CenterY + RPF2), Bitmap.Height - 1);

// Optimization run: find squares of X first
SetLength(SqX, RX - LX + 1);
for x := LX to RX do
SqX[x - LX] := sqr(x - CenterX);

// Loop through Y values
for y := LY to RY do begin
P := Bitmap.Scanline[y];
SqY := Sqr(y - CenterY);
// Loop through X values
for x := LX to RX do begin

// determine squared distance from center for this pixel
SqDist := SqY + SqX[x - LX];

// inside inner circle? Most often..
if sqdist < RMF2 then begin
// inside the inner circle.. just give the scanline the
// new color
P[x] := 255
end else begin
// inside outer circle?
if sqdist < RPF2 then begin
// We are inbetween the inner and outer bound, now
// mix the color
Fact := round(((Radius - sqrt(sqdist)) * 2 / Feather)
* 127.5 + 127.5);
// just in case limit to [0, 255]
P[x] := Max(0, Min(Fact, 255));
end else begin
P[x] := 0;
end;
end;
end;
end;
end;```