+ Coordinate_Product type for extended range computations ans Ellipse cleanup

This commit is contained in:
Vovanium 2023-05-22 12:30:45 +03:00
parent 1891961f63
commit 135317eb5a
2 changed files with 37 additions and 31 deletions

View File

@ -2,6 +2,9 @@ package Video.Integer_Geometry with Pure is
subtype Coordinate is Integer;
type Coordinate_Product is range -(2 * Coordinate'First**2) .. 2 * Coordinate'Last**2;
-- An extended range type for geometry computations
subtype Distance is Coordinate'Base range 0 .. (
if Coordinate'Last > Coordinate'Base'Last + Coordinate'First then
Coordinate'Base'Last else Coordinate'Last - Coordinate'First);

View File

@ -104,10 +104,12 @@ package body Video.Rasters.Generic_Blits is
Target : in out Raster;
Bounds : in Integer_Geometry.Box)
is
X, Y : Coordinate; -- Offset from center in half-pixels
D : Coordinate; -- Deviation from ideal ellipse
A, B : Coordinate; -- Axis lengths
S : Coordinate;
subtype Prod is Coordinate_Product;
X, Y : Coordinate; -- Offset from center in half-pixels
D : Prod; -- Deviation from ideal ellipse
A, B : Distance; -- Axis lengths
A2, B2 : Prod; -- Squares of axis length
S : Prod;
procedure Draw is
XS : constant Coordinate := Bounds.X.First + Bounds.X.Last;
YS : constant Coordinate := Bounds.Y.First + Bounds.Y.Last;
@ -122,8 +124,8 @@ package body Video.Rasters.Generic_Blits is
return; -- No pixels to draw
end if;
A := Length (Bounds.X); -- +1 to have 1/2 pixel offset
B := Length (Bounds.Y);
A := Length (Bounds.X); A2 := Prod (A)**2;
B := Length (Bounds.Y); B2 := Prod (B)**2;
-- Horizontal part
if A mod 2 = 0 then
@ -137,45 +139,46 @@ package body Video.Rasters.Generic_Blits is
Y := B;
-- D := ((DX * B)**2 + (DY * A)**2 - A**2 * B) / 4;
D := ((X * B)**2) / 4; -- + 2 is for rounding
while (X + 1) * B**2 < (Y - 1) * A**2 loop
D := D + (X + 1) * B**2;
D := ((Prod (X))**2 * B2 + 2) / 4; -- + 2 is for rounding
while Prod (X) * B2 < Prod (Y) * A2 loop
D := D + (Prod (X) + 1) * B2;
X := X + 2;
pragma Assert (D = (X * B)**2 + (Y * A)**2 - A**2 * B**2);
--pragma Assert (D = (X * B)**2 + (Y * A)**2 - (A * B)**2);
-- Step Y -> Y - 2
-- D(i) = ((X * B)**2 + ( Y * A)**2 - A**2 * B**2) / 4
-- D(i+1) = ((X * B)**2 + ((Y - 2) * A)**2 - A**2 * B**2) / 4
-- D(i+1) - D(i) = (((Y - 2) * A)**2) / 4 - ((Y * A)**2) / 4
-- = A**2 * (1 - Y)
S := (1 - Y) * A**2;
S := (1 - Prod (Y)) * A2;
-- Step down when its deviation is of opposite sign and less in abs. value
if (D + S / 2) >= 0 then
if D + S >= -D then
Y := Y - 2;
D := D + S;
pragma Assert (D = (X * B)**2 + (Y * A)**2 - A**2 * B**2);
end if;
Draw;
end loop;
while Y > 0 loop
D := D + (1 - Y) * A**2;
Y := Y - 2;
pragma Assert (D = (X * B)**2 + (Y * A)**2 - A**2 * B**2);
-- Step X -> X + 2
-- D(i) = (( X * B)**2 + (Y * A)**2 - A**2 * B**2) / 4
-- D(i+1) = (((X + 2) * B)**2 + (Y * A)**2 - A**2 * B**2) / 4
-- D(i+1) - D(i) = (((X + 2) * B)**2) / 4 - ((X * B)**2) / 4
-- = B**2 * (X + 1)
S := (X + 1) * B**2;
if (D + S / 2) < 0 then
X := X + 2;
D := D + S;
pragma Assert (D = (X * B)**2 + (Y * A)**2 - A**2 * B**2);
--pragma Assert (D = (X * B)**2 + (Y * A)**2 - (A * B)**2);
end if;
Draw;
end loop;
-- Vertical part
while Y > 2 loop
D := D + (1 - Prod (Y)) * A2;
Y := Y - 2;
--pragma Assert (D = (X * B)**2 + (Y * A)**2 - (A * B)**2);
-- Step X -> X + 2
-- D(i) = (( X * B)**2 + (Y * A)**2 - A**2 * B**2) / 4
-- D(i+1) = (((X + 2) * B)**2 + (Y * A)**2 - A**2 * B**2) / 4
-- D(i+1) - D(i) = (((X + 2) * B)**2) / 4 - ((X * B)**2) / 4
-- = B**2 * (X + 1)
S := (Prod (X) + 1) * B2;
if D + S / 2 < -D then
X := X + 2;
D := D + S;
--pragma Assert (D = (X * B)**2 + (Y * A)**2 - (A * B)**2);
end if;
Draw;
end loop;
if B mod 2 = 0 then
Target (Bounds.Y.First, Center (Bounds.X)) := Color;
Target (Bounds.Y.Last, Center (Bounds.X)) := Color;