procedure ConvexHull(var TPA: TPointArray);
var
  p, Lower: TPointArray;
  LH, H, I, UH, c, j, idx, x, y: Integer;
  area: TBox;
  m: array of TBoolArray;
  b: Boolean;
begin
  h := High(TPA);
  if (h > 0) then
  begin
    area := GetTPABounds(TPA);
    SetLength(m, ((area.X2 - area.X1) + 1));
    for i := 0 to (area.X2 - area.X1) do
      SetLength(m[i], ((area.Y2 - area.Y1) + 1));
    C := 0;
    for i := 0 to h do
    begin
      x := (TPA[i].X - area.X1);
      y := (TPA[i].Y - area.Y1);
      if m[x][y] then
        Continue;
      m[x][y] := True;
      Inc(C);
    end;
    SetLength(p, C);
    idx := 0;
    b := False;
    for y := area.Y1 to area.Y2 do
    begin
      for x := area.X1 to area.X2 do
        if m[(x - area.X1)][(y - area.Y1)] then
        begin
          p[idx] := Point(x, y);
          Inc(idx);
          b := (idx >= C);
          if b then
            Break;
        end;
      if b then
        Break;
    end;
    h := High(p);
    if (h > 0) then
    begin
      UH := 2;
      SetLength(TPA, (h + 1));
      TPA[0] := p[0];
      TPA[1] := p[1];
      for i := 2 to h do
      begin
        TPA[UH] := p[i];
        Inc(UH);
        while ((UH > 2) and not (((TPA[(UH - 2)].x * TPA[(UH - 1)].y + TPA[(UH - 3)].x * TPA[(UH - 2)].y + TPA[(UH - 1)].x * TPA[(UH - 3)].y) - (TPA[(UH - 2)].x * TPA[(UH - 3)].y + TPA[(UH - 1)].x * TPA[(UH - 2)].y + TPA[(UH - 3)].x * TPA[(UH - 1)].y)) < 0)) do
        begin
          Dec(UH);
          TPA[(UH - 1)] := TPA[UH];
        end;
      end;
      LH := 2;
      SetLength(Lower, (h + 1));
      Lower[0] := p[h];
      Lower[1] := p[(h - 1)];
      for i := 2 to h do
      begin
        Lower[LH] := p[(h - i)];
        Inc(LH);
        while ((LH > 2) and not (((Lower[(LH - 2)].x * Lower[(LH - 1)].y + Lower[(LH - 3)].x * Lower[(LH - 2)].y + Lower[(LH - 1)].x * Lower[(LH - 3)].y) - (Lower[(LH - 2)].x * Lower[(LH - 3)].y + Lower[(LH - 1)].x * Lower[(LH - 2)].y + Lower[(LH - 3)].x * Lower[(LH - 1)].y)) < 0)) do
        begin
          Dec(LH);
          Lower[(LH - 1)] := Lower[LH];
        end;
      end;
      Dec(LH);
      SetLength(TPA, (UH + LH));
      for i := UH to ((UH + LH) - 1) do
        TPA[i] := Lower[(i - UH)];
    end;
  end;
end;

{==============================================================================]
  @action: Outputs polygon by given outlines (shape)
  @note: None
  @contributors: slacky, Janilabo
[==============================================================================}
function TPAFromPolygon(shape: TPointArray): TPointArray;
var
  b: TBox;
  x, y, h, i, l, r, z: Integer;
  o: TPointArray;
  f: Boolean;
  t: array of TBoolArray;
  e: Extended;
  q, p, d: TPoint;
begin
  h := High(shape);
  if (h > -1) then
  begin
    SetLength(o, 0);
    b := IntToBox(shape[0].X, shape[0].Y, shape[0].X, shape[0].Y);
    for i := 0 to h do
    begin
      q := shape[i];
      if (i < h) then
        p := shape[(i + 1)]
      else
        p := shape[0];
      r := Length(o);
      if ((q.X <> p.X) or (q.Y <> p.Y)) then
      begin
        l := Max(Round(Abs(q.X - p.X)), Round(Abs(q.Y - p.Y)));
        SetLength(o, ((r + l) + 1));
        for z := 0 to l do
          o[(r + z)] := Point((q.X + Round((p.X - q.X) * (z / Extended(l)))), (q.Y + Round((p.Y - q.Y) * (z / Extended(l)))));
      end else
      begin
        SetLength(o, (r + 1));
        o[r] := q;
      end;
      if (shape[i].X < b.X1) then
        b.X1 := shape[i].X
      else
        if (shape[i].X > b.X2) then
          b.X2 := shape[i].X;
      if (shape[i].Y < b.Y1) then
        b.Y1 := shape[i].Y
      else
        if (shape[i].Y > b.Y2) then
          b.Y2 := shape[i].Y;
    end;
    SetLength(t, ((b.X2 - b.X1) + 1));
    for i := 0 to (b.X2 - b.X1) do
      SetLength(t[i], ((b.Y2 - b.Y1) + 1));
    l := Length(o);
    for i := 0 to (l - 1) do
      if not t[(o[i].X - b.X1)][(o[i].Y - b.Y1)] then
        t[(o[i].X - b.X1)][(o[i].Y - b.Y1)] := True;
    for y := 0 to (b.Y2 - b.Y1) do
      for x := 0 to (b.X2 - b.X1) do
      begin
        f := t[x][y];
        if not f then
        begin
          d := Point((x + b.X1), (y + b.Y1));
          q := shape[0];
          for i := 0 to (h + 1) do
          begin
            p := shape[(i mod (h + 1))];
            if (d.Y > Min(q.Y, p.Y)) then
              if (d.Y <= Max(q.Y, p.Y)) then
                if (d.X <= Max(q.X, p.X)) then
                begin
                  if (q.y <> p.y) then
                    e := ((d.Y - q.Y) * (p.X - q.X) / Extended((p.Y - q.Y)) + q.X);
                  if ((q.X = p.X) or (d.X < e)) then
                    f := not f;
                end;
            q := p;
          end;
        end;
        if f then
        begin
          l := Length(Result);
          SetLength(Result, (l + 1));
          Result[l] := Point((x + b.X1), (y + b.Y1));
        end;
      end;
  end else
    SetLength(Result, 0);
end;

//******* Plot a few shapes using this algorithm *******//
var
  TPA,CTPA: TPointArray;
  bmp, t: Integer;
begin
  bmp := CreateBitmap(700, 700);

  TPA := TPAFromCircle(350, 350, 150);
  CTPA := CopyTPA(TPA);
  SetLength(TPA, 0);
  t := GetSystemTime;
  ConvexHull(CTPA);
  WriteLn('ConvexHull took ' + IntToStr(GetSystemTime - t) + ' ms.');
  t := GetSystemTime;
  TPA := TPAFromPolygon(CTPA);
  WriteLn('TPAFromPolygon took ' + IntToStr(GetSystemTime - t) + ' ms.');
  DrawTPABitmap(bmp, TPA, 357435);
  DrawTPABitmap(bmp, CTPA, 255);

  DisplayDebugImgWindow(700,700);
  DrawBitmapDebugImg(bmp);
  FreeBitmap(bmp);
end.