mirrors/doublecmd
2
0
Fork 0
mirror of https://github.com/doublecmd/doublecmd.git synced 2026-06-21 09:58:13 +00:00
Code Issues Projects Releases Packages Wiki Activity
doublecmd/components/Image32/source/Img32.Transform.pas
Alexander Koblov 656b459b2a UPD: Image32 library
2025-03-23 21:11:06 +03:00

1375 lines
42 KiB
ObjectPascal
Raw Blame History

unit Img32.Transform;
(*******************************************************************************
* Author : Angus Johnson *
* Version : 4.7 *
* Date : 6 January 2025 *
* Website : http://www.angusj.com *
* Copyright : Angus Johnson 2019-2025 *
* Purpose : Affine and projective transformation routines for TImage32 *
* License : http://www.boost.org/LICENSE_1_0.txt *
*******************************************************************************)
interface
{$I Img32.inc}
uses
SysUtils, Classes, Math, Types, Img32, Img32.Vector;
type
TMatrixD = array [0..2, 0..2] of double;
//Matrix functions
function IsIdentityMatrix(const matrix: TMatrixD): Boolean;
{$IFDEF INLINE} inline; {$ENDIF}
function IsValidMatrix(const matrix: TMatrixD): Boolean;
{$IFDEF INLINE} inline; {$ENDIF}
function Matrix(const m00, m01, m02, m10, m11, m12, m20, m21, m22: double): TMatrixD;
function MatrixDeterminant(const matrix: TMatrixD): double;
function MatrixAdjugate(const matrix: TMatrixD): TMatrixD;
function MatrixInvert(var matrix: TMatrixD): Boolean;
// Note: Matrix multiplication IS NOT commutative hence ...
procedure MatrixMultiply(var matrix1: TMatrixD; const matrix2: TMatrixD);
procedure MatrixMultiply2(const matrix1: TMatrixD; var matrix2: TMatrixD);
procedure MatrixApply(const matrix: TMatrixD;
var x, y: double); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure MatrixApply(const matrix: TMatrixD;
var pt: TPointD); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure MatrixApply(const matrix: TMatrixD; var rec: TRect); overload;
procedure MatrixApply(const matrix: TMatrixD; var rec: TRectD); overload;
procedure MatrixApply(const matrix: TMatrixD; var path: TPathD); overload;
procedure MatrixApply(const matrix: TMatrixD; var paths: TPathsD); overload;
procedure MatrixApply(const matrix: TMatrixD;
img: TImage32; scaleAdjust: Boolean = false); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure MatrixApply(const matrix: TMatrixD;
img, targetImg: TImage32; scaleAdjust: Boolean = false); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure MatrixSkew(var matrix: TMatrixD; angleX, angleY: double);
procedure MatrixScale(var matrix: TMatrixD; scale: double); overload;
procedure MatrixScale(var matrix: TMatrixD; scaleX, scaleY: double); overload;
procedure MatrixRotate(var matrix: TMatrixD; angRad: double); overload;
procedure MatrixRotate(var matrix: TMatrixD; const center: TPointD; angRad: double); overload;
procedure MatrixTranslate(var matrix: TMatrixD; dx, dy: double);
// The following MatrixExtract routines assume here is no skew
procedure MatrixExtractScale(const mat: TMatrixD; out scale: double); overload;
procedure MatrixExtractScale(const mat: TMatrixD; out X, Y: double); overload;
procedure MatrixExtractTranslation(const mat: TMatrixD; out dx, dy: double);
procedure MatrixExtractRotation(const mat: TMatrixD; out angle: double);
// MatrixExtractAll - except skew :)
function MatrixExtractAll(const mat: TMatrixD; out angle: double;
out scale, trans: TPointD): Boolean;
// AffineTransformImage: will automagically translate the image
// Note: when the "scaleAdjust" parameter is enabled, it prevents antialiasing
// from extending way outside of images when they are being enlarged
// significantly (> 2 times) and rotated concurrently
function AffineTransformImage(img: TImage32; matrix: TMatrixD;
scaleAdjust: Boolean = false): TPoint; overload; {$IFDEF INLINE} inline; {$ENDIF}
function AffineTransformImage(img, targetImg: TImage32; matrix: TMatrixD;
scaleAdjust: Boolean = false): TPoint; overload;
// ProjectiveTransform:
// srcPts, dstPts => each path must contain 4 points
// margins => the margins around dstPts (in the dest. projective).
// Margins are only meaningful when srcPts are inside the image.
function ProjectiveTransform(img: TImage32;
const srcPts, dstPts: TPathD; const margins: TRect): Boolean; overload; {$IFDEF INLINE} inline; {$ENDIF}
function ProjectiveTransform(img, targetImg: TImage32;
const srcPts, dstPts: TPathD; const margins: TRect): Boolean; overload;
function SplineVertTransform(img: TImage32; const topSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean; overload; {$IFDEF INLINE} inline; {$ENDIF}
function SplineVertTransform(img, targetImg: TImage32; const topSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean; overload;
function SplineHorzTransform(img: TImage32; const leftSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean; overload; {$IFDEF INLINE} inline; {$ENDIF}
function SplineHorzTransform(img, targetImg: TImage32; const leftSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean; overload;
type
PWeightedColor = ^TWeightedColor;
TWeightedColor = {$IFDEF RECORD_METHODS} record {$ELSE} object {$ENDIF}
private
fAddCount : Integer;
fAlphaTot : Int64;
fColorTotR: Int64;
fColorTotG: Int64;
fColorTotB: Int64;
function GetColor: TColor32;
public
procedure Reset; overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure Reset(c: TColor32; w: Integer = 1); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure Add(c: TColor32; w: Integer); overload; {$IFDEF INLINE_COMPATIBLE} inline; {$ENDIF}
procedure Add(c: TColor32); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure Add(const other: TWeightedColor); overload;
{$IFDEF INLINE} inline; {$ENDIF}
procedure Subtract(c: TColor32; w: Integer); overload;
procedure Subtract(c: TColor32); overload; {$IFDEF INLINE} inline; {$ENDIF}
procedure Subtract(const other: TWeightedColor); overload;
{$IFDEF INLINE} inline; {$ENDIF}
function AddSubtract(addC, subC: TColor32): Boolean; {$IFDEF INLINE_COMPATIBLE} inline; {$ENDIF}
function AddNoneSubtract(c: TColor32): Boolean; {$IFDEF INLINE_COMPATIBLE} inline; {$ENDIF}
procedure AddWeight(w: Integer); {$IFDEF INLINE} inline; {$ENDIF}
property AddCount: Integer read fAddCount;
property Color: TColor32 read GetColor;
property Weight: integer read fAddCount;
end;
TArrayOfWeightedColor = array of TWeightedColor;
const
IdentityMatrix: TMatrixD = ((1, 0, 0),(0, 1, 0),(0, 0, 1));
implementation
uses Img32.Resamplers;
resourcestring
rsInvalidScale = 'Invalid matrix scaling factor (0)';
const
DivOneByXTableSize = 1024;
{$IFDEF CPUX86}
// Use faster Trunc for x86 code in this unit.
Trunc: function(Value: Double): Integer = __Trunc;
{$ENDIF CPUX86}
var
// DivOneByXTable[x] = 1/x
DivOneByXTable: array[0 .. DivOneByXTableSize -1] of Double;
//------------------------------------------------------------------------------
// Matrix functions
//------------------------------------------------------------------------------
function IsIdentityMatrix(const matrix: TMatrixD): Boolean;
begin
result := (matrix[0,0] = 1) and (matrix[0,1] = 0) and (matrix[0,2] = 0) and
(matrix[1,0] = 0) and (matrix[1,1] = 1) and (matrix[1,2] = 0) and
(matrix[2,0] = 0) and (matrix[2,1] = 0) and (matrix[2,2] = 1);
end;
//------------------------------------------------------------------------------
function IsValidMatrix(const matrix: TMatrixD): Boolean;
begin
result := matrix[2][2] = 1.0;
end;
//------------------------------------------------------------------------------
function Matrix(const m00, m01, m02, m10, m11, m12, m20, m21, m22: double): TMatrixD;
begin
Result[0,0] := m00; Result[0,1] := m01; Result[0,2] := m02;
Result[1,0] := m10; Result[1,1] := m11; Result[1,2] := m12;
Result[2,0] := m20; Result[2,1] := m21; Result[2,2] := m22;
end;
//------------------------------------------------------------------------------
function Det4(a1, a2, b1, b2: double): double; {$IFDEF INLINE} inline; {$ENDIF}
begin
Result := a1 * b2 - a2 * b1;
end;
//------------------------------------------------------------------------------
function Det9(a1, a2, a3, b1, b2, b3, c1, c2, c3: double): double;
{$IFDEF INLINE} inline; {$ENDIF}
begin
Result := a1 * Det4(b2, b3, c2, c3) -
b1 * Det4(a2, a3, c2, c3) +
c1 * Det4(a2, a3, b2, b3);
end;
//------------------------------------------------------------------------------
function MatrixDeterminant(const matrix: TMatrixD): double;
{$IFDEF INLINE} inline; {$ENDIF}
begin
Result := Det9(matrix[0,0], matrix[1,0], matrix[2,0],
matrix[0,1], matrix[1,1], matrix[2,1],
matrix[0,2], matrix[1,2], matrix[2,2]);
end;
//------------------------------------------------------------------------------
function MatrixAdjugate(const matrix: TMatrixD): TMatrixD;
begin
//https://en.wikipedia.org/wiki/Adjugate_matrix
Result[0,0] := Det4(matrix[1,1], matrix[1,2], matrix[2,1], matrix[2,2]);
Result[0,1] := -Det4(matrix[0,1], matrix[0,2], matrix[2,1], matrix[2,2]);
Result[0,2] := Det4(matrix[0,1], matrix[0,2], matrix[1,1], matrix[1,2]);
Result[1,0] := -Det4(matrix[1,0], matrix[1,2], matrix[2,0], matrix[2,2]);
Result[1,1] := Det4(matrix[0,0], matrix[0,2], matrix[2,0], matrix[2,2]);
Result[1,2] := -Det4(matrix[0,0], matrix[0,2], matrix[1,0], matrix[1,2]);
Result[2,0] := Det4(matrix[1,0], matrix[1,1], matrix[2,0], matrix[2,1]);
Result[2,1] := -Det4(matrix[0,0], matrix[0,1], matrix[2,0], matrix[2,1]);
Result[2,2] := Det4(matrix[0,0], matrix[0,1], matrix[1,0], matrix[1,1]);
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var x, y: double);
var
tmpX: double;
begin
tmpX := x;
x := tmpX * matrix[0, 0] + y * matrix[1, 0] + matrix[2, 0];
y := tmpX * matrix[0, 1] + y * matrix[1, 1] + matrix[2, 1];
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var pt: TPointD);
var
tmpX: double;
begin
tmpX := pt.x;
pt.X := tmpX * matrix[0, 0] + pt.Y * matrix[1, 0] + matrix[2, 0];
pt.Y := tmpX * matrix[0, 1] + pt.Y * matrix[1, 1] + matrix[2, 1];
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var rec: TRect);
var
path: TPathD;
begin
if not IsValidMatrix(matrix) then Exit;
path := Rectangle(rec);
MatrixApply(matrix, path);
rec := GetBounds(path);
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var rec: TRectD);
var
path: TPathD;
begin
if not IsValidMatrix(matrix) then Exit;
path := Rectangle(rec);
MatrixApply(matrix, path);
rec := GetBoundsD(path);
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var path: TPathD);
var
i, len: integer;
tmpX: double;
pp: PPointD;
begin
len := Length(path);
if (len = 0) or IsIdentityMatrix(matrix) or
not IsValidMatrix(matrix) then Exit;
pp := @path[0];
for i := 0 to len -1 do
begin
tmpX := pp.X;
pp.X := tmpX * matrix[0, 0] + pp.Y * matrix[1, 0] + matrix[2, 0];
pp.Y := tmpX * matrix[0, 1] + pp.Y * matrix[1, 1] + matrix[2, 1];
inc(pp);
end;
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD; var paths: TPathsD);
var
i,j,len: integer;
tmpX: double;
pp: PPointD;
begin
if not Assigned(paths) or IsIdentityMatrix(matrix) or
not IsValidMatrix(matrix) then Exit;
for i := 0 to High(paths) do
begin
len := Length(paths[i]);
if len = 0 then Continue;
pp := @paths[i][0];
for j := 0 to High(paths[i]) do
begin
tmpX := pp.X;
pp.X := tmpX * matrix[0, 0] + pp.Y * matrix[1, 0] + matrix[2, 0];
pp.Y := tmpX * matrix[0, 1] + pp.Y * matrix[1, 1] + matrix[2, 1];
inc(pp);
end;
end;
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD;
img: TImage32; scaleAdjust: Boolean);
begin
AffineTransformImage(img, matrix, scaleAdjust);
end;
//------------------------------------------------------------------------------
procedure MatrixApply(const matrix: TMatrixD;
img, targetImg: TImage32; scaleAdjust: Boolean);
begin
AffineTransformImage(img, targetImg, matrix, scaleAdjust);
end;
//------------------------------------------------------------------------------
procedure MatrixMultiply(var matrix1: TMatrixD; const matrix2: TMatrixD);
var
i, j: Integer;
m: TMatrixD;
begin
for i := 0 to 2 do
for j := 0 to 2 do
m[i, j] :=
(matrix1[i, 0] * matrix2[0, j]) +
(matrix1[i, 1] * matrix2[1, j]) +
(matrix1[i, 2] * matrix2[2, j]);
matrix1 := m;
end;
//------------------------------------------------------------------------------
procedure MatrixMultiply2(const matrix1: TMatrixD; var matrix2: TMatrixD);
var
i, j: Integer;
m: TMatrixD;
begin
for i := 0 to 2 do
for j := 0 to 2 do
m[i, j] :=
(matrix1[i, 0] * matrix2[0, j]) +
(matrix1[i, 1] * matrix2[1, j]) +
(matrix1[i, 2] * matrix2[2, j]);
matrix2 := m;
end;
//------------------------------------------------------------------------------
procedure MatrixScale(var matrix: TMatrixD; scaleX, scaleY: double);
var
m: TMatrixD;
begin
m := IdentityMatrix;
if (scaleX = 0) or (scaleY = 0) then
raise Exception(rsInvalidScale);
if ValueAlmostOne(scaleX) and ValueAlmostOne(scaleY) then Exit;
m[0, 0] := scaleX;
m[1, 1] := scaleY;
MatrixMultiply(matrix, m);
end;
//------------------------------------------------------------------------------
procedure MatrixScale(var matrix: TMatrixD; scale: double);
begin
if (scale = 0) or (scale = 1) then Exit;
MatrixScale(matrix, scale, scale);
end;
//------------------------------------------------------------------------------
procedure MatrixRotate(var matrix: TMatrixD;
const center: TPointD; angRad: double);
var
m: TMatrixD;
sinA, cosA: double;
begin
if not PointsEqual(center, NullPointD) then
begin
NormalizeAngle(angRad);
if angRad = 0 then Exit;
{$IFNDEF CLOCKWISE_ROTATION_WITH_NEGATIVE_ANGLES}
angRad := -angRad; //negated angle because of inverted Y-axis.
{$ENDIF}
m := IdentityMatrix;
MatrixTranslate(matrix, -center.X, -center.Y);
GetSinCos(angRad, sinA, cosA);
m := IdentityMatrix;
m[0, 0] := cosA; m[1, 0] := sinA;
m[0, 1] := -sinA; m[1, 1] := cosA;
MatrixMultiply(matrix, m);
MatrixTranslate(matrix, center.X, center.Y);
end
else
MatrixRotate(matrix, angRad)
end;
//------------------------------------------------------------------------------
procedure MatrixRotate(var matrix: TMatrixD; angRad: double);
var
m: TMatrixD;
sinA, cosA: double;
begin
NormalizeAngle(angRad);
if angRad = 0 then Exit;
{$IFNDEF CLOCKWISE_ROTATION_WITH_NEGATIVE_ANGLES}
angRad := -angRad; //negated angle because of inverted Y-axis.
{$ENDIF}
m := IdentityMatrix;
GetSinCos(angRad, sinA, cosA);
m := IdentityMatrix;
m[0, 0] := cosA; m[1, 0] := sinA;
m[0, 1] := -sinA; m[1, 1] := cosA;
MatrixMultiply(matrix, m);
end;
//------------------------------------------------------------------------------
procedure MatrixTranslate(var matrix: TMatrixD; dx, dy: double);
var
m: TMatrixD;
begin
if ValueAlmostZero(dx) and ValueAlmostZero(dy) then Exit;
m := IdentityMatrix;
m[2, 0] := dx;
m[2, 1] := dy;
MatrixMultiply(matrix, m);
end;
//------------------------------------------------------------------------------
procedure ScaleInternal(var matrix: TMatrixD; s: double);
var
i, j: Integer;
begin
for i := 0 to 2 do
for j := 0 to 2 do
matrix[i,j] := matrix[i,j] * s;
end;
//------------------------------------------------------------------------------
function MatrixInvert(var matrix: TMatrixD): Boolean;
var
d: double;
const
tolerance = 1.0E-5;
begin
d := MatrixDeterminant(matrix);
Result := abs(d) > tolerance;
if Result then
begin
matrix := MatrixAdjugate(matrix);
ScaleInternal(matrix, 1/d);
end;
end;
//------------------------------------------------------------------------------
procedure MatrixSkew(var matrix: TMatrixD; angleX, angleY: double);
var
m: TMatrixD;
begin
if ValueAlmostZero(angleX) and ValueAlmostZero(angleY) then Exit;
m := IdentityMatrix;
m[1, 0] := tan(angleX);
m[0, 1] := tan(angleY);
MatrixMultiply(matrix, m);
end;
//------------------------------------------------------------------------------
procedure MatrixExtractScale(const mat: TMatrixD; out X, Y: double);
begin
// https://stackoverflow.com/a/32125700/359538
X := Sqrt(Sqr(mat[0,0]) + Sqr(mat[0,1]));
//Y := Sqrt(Sqr(mat[1,0]) + Sqr(mat[1,1]));
Y := Abs((mat[0,0] * mat[1,1] - mat[1,0] * mat[0,1]) / X);
end;
//------------------------------------------------------------------------------
procedure MatrixExtractScale(const mat: TMatrixD; out scale: double);
var
x,y: double;
begin
MatrixExtractScale(mat, x, y);
scale := Average(x,y);
end;
//------------------------------------------------------------------------------
procedure MatrixExtractTranslation(const mat: TMatrixD; out dx, dy: double);
begin
dx := mat[2,0];
dy := mat[2,1];
end;
//------------------------------------------------------------------------------
procedure MatrixExtractRotation(const mat: TMatrixD; out angle: double);
begin
angle := ArcTan2(mat[0,1], mat[0,0]);
end;
//------------------------------------------------------------------------------
function MatrixExtractAll(const mat: TMatrixD;
out angle: double; out scale, trans: TPointD): Boolean;
var
m00, m01, m10, m11: double;
begin
m00 := mat[0][0]; m10 := mat[1][0];
m01 := mat[0][1]; m11 := mat[1][1];
trans.X := mat[2][0];
trans.Y := mat[2][1];
angle := 0;
scale := PointD(1,1);
Result := (m00 <> 0) or (m01 <> 0);
if not Result then Exit;
angle := ArcTan2(m01, m00);
// https://stackoverflow.com/a/32125700/359538
scale.X := Sqrt(Sqr(mat[0,0]) + Sqr(mat[0,1]));
scale.Y := (m00 * m11 - m10 * m01) / scale.X;
end;
//------------------------------------------------------------------------------
{$IFDEF USE_DOWNSAMPLER_AUTOMATICALLY}
function CanUseBoxDownsampler(const mat: TMatrixD; sx, sy: double): Boolean;
begin
// If the matrix looks like this after removing the scale,
// the box downsampler can be used. (only translation and scaling)
// cos(0) -sin(0) tx 1 0 tx
// sin(0) cos(0) ty => 0 1 ty
// 0 0 1 0 0 1
{
Result := (mat[0,0]/sx = 1) and (mat[0,1]/sx = 0) and
(mat[1,0]/sy = 0) and (mat[1,1]/sy = 1) and
(mat[2,0] = 0) and (mat[2,1] = 0) and
(mat[2,2] = 1);
}
// We can skip the divisions, because m/s is only zero if m is zero
// and m/s=1 is the same as m=s
Result := (SameValue(mat[0,0], sx)) and (mat[0,1] = 0) and
(mat[1,0] = 0) and (SameValue(mat[1,1], sy)) and
(mat[2,0] = 0) and (mat[2,1] = 0) and
(mat[2,2] = 1);
end;
{$ENDIF USE_DOWNSAMPLER_AUTOMATICALLY}
//------------------------------------------------------------------------------
// Affine Transformation
//------------------------------------------------------------------------------
function GetTransformBounds(img: TImage32; const matrix: TMatrixD): TRect;
var
pts: TPathD;
begin
pts := Rectangle(img.Bounds);
MatrixApply(matrix, pts);
Result := GetBounds(pts);
end;
//------------------------------------------------------------------------------
function AffineTransformImage(img: TImage32; matrix: TMatrixD;
scaleAdjust: Boolean): TPoint;
begin
Result := AffineTransformImage(img, img, matrix, scaleAdjust);
end;
//------------------------------------------------------------------------------
function AffineTransformImage(img, targetImg: TImage32; matrix: TMatrixD;
scaleAdjust: Boolean): TPoint;
var
i, j: integer;
newWidth, newHeight: integer;
sx, sy, x,y: double;
xLimLo, yLimLo, xLimHi, yLimHi: double;
pc: PColor32;
tmp: TArrayOfColor32;
dstRec: TRect;
resampler: TResamplerFunction;
{$IFDEF USE_DOWNSAMPLER_AUTOMATICALLY}
useBoxDownsampler: Boolean;
{$ENDIF}
begin
Result := NullPoint;
if IsIdentityMatrix(matrix) or
img.IsEmpty or (targetImg.Resampler = 0) then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
resampler := GetResampler(targetImg.Resampler);
if not Assigned(resampler) then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
//auto-resize the image so it'll fit transformed image
dstRec := img.Bounds;
MatrixApply(matrix, dstRec);
RectWidthHeight(dstRec, newWidth, newHeight);
MatrixExtractScale(matrix, sx, sy);
{$IFDEF USE_DOWNSAMPLER_AUTOMATICALLY}
if (sx < 1.0) and (sy < 1.0) then
begin
//only use box downsampling when downsizing
useBoxDownsampler := CanUseBoxDownsampler(matrix, sx, sy);
end else
useBoxDownsampler := false;
if useBoxDownsampler then
begin
BoxDownSampling(img, targetImg, sx, sy);
Exit;
end;
{$ENDIF}
if scaleAdjust then
begin
sx := Max(1, sx * 0.5);
sy := Max(1, sy * 0.5);
end;
//auto-translate the image too
Result := dstRec.TopLeft;
//starting with the result pixel coords, reverse lookup
//the fractional coordinates in the untransformed image
if not MatrixInvert(matrix) then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
NewColor32Array(tmp, newWidth * newHeight, True);
pc := @tmp[0];
xLimLo := -0.5/sx;
xLimHi := img.Width + 0.5/sx;
yLimLo := -0.5/sy;
yLimHi := img.Height + 0.5/sy;
for i := dstRec.Top to dstRec.Bottom -1 do
begin
for j := dstRec.Left to dstRec.Right -1 do
begin
x := j; y := i;
MatrixApply(matrix, x, y);
if (x <= xLimLo) or (x >= xLimHi) or (y <= yLimLo) or (y >= yLimHi) then
pc^ := clNone32
else
// nb: -0.5 below is needed to properly center the transformed image
// (and this is most obviously needed when there is large scaling)
pc^ := resampler(img, x - 0.5, y - 0.5);
inc(pc);
end;
end;
targetImg.AssignPixelArray(tmp, newWidth, newHeight);
end;
//------------------------------------------------------------------------------
// Projective Transformation
//------------------------------------------------------------------------------
procedure MatrixMulCoord(const matrix: TMatrixD; var x,y,z: double);
{$IFDEF INLINE} inline; {$ENDIF}
var
xx, yy: double;
begin
xx := x; yy := y;
x := matrix[0,0] *xx + matrix[0,1] *yy + matrix[0,2] *z;
y := matrix[1,0] *xx + matrix[1,1] *yy + matrix[1,2] *z;
z := matrix[2,0] *xx + matrix[2,1] *yy + matrix[2,2] *z;
end;
//------------------------------------------------------------------------------
function BasisToPoints(x1, y1, x2, y2, x3, y3, x4, y4: double): TMatrixD;
var
m2: TMatrixD;
z4: double;
begin
Result := Matrix(x1, x2, x3, y1, y2, y3, 1, 1, 1);
m2 := MatrixAdjugate(Result);
z4 := 1;
MatrixMulCoord(m2, x4, y4, z4);
m2 := Matrix(x4, 0, 0, 0, y4, 0, 0, 0, z4);
MatrixMultiply(Result, m2);
end;
//------------------------------------------------------------------------------
procedure GetSrcCoords256(const matrix: TMatrixD; var x, y: integer);
{$IFDEF INLINE} inline; {$ENDIF}
var
xx,yy,zz: double;
const
Q: integer = MaxInt div 256;
begin
//returns coords multiplied by 256 in anticipation of the following
//GetWeightedPixel function call which in turn expects the lower 8bits
//of the integer coord value to represent a fraction.
xx := x; yy := y; zz := 1;
MatrixMulCoord(matrix, xx, yy, zz);
if zz = 0 then
begin
if xx >= 0 then x := Q else x := -MaxInt;
if yy >= 0 then y := Q else y := -MaxInt;
end else
begin
xx := xx/zz;
if xx > Q then x := MaxInt
else if xx < -Q then x := -MaxInt
else x := Round(xx *256);
yy := yy/zz;
if yy > Q then y := MaxInt
else if yy < -Q then y := -MaxInt
else y := Round(yy *256);
end;
end;
//------------------------------------------------------------------------------
procedure GetSrcCoords(const matrix: TMatrixD; var x, y: double);
{$IFDEF INLINE} inline; {$ENDIF}
var
zz: double;
const
Q: integer = MaxInt div 256;
begin
//returns coords multiplied by 256 in anticipation of the following
//GetWeightedPixel function call which in turn expects the lower 8bits
//of the integer coord value to represent a fraction.
zz := 1;
MatrixMulCoord(matrix, x, y, zz);
if zz = 0 then
begin
if x >= 0 then x := Q else x := -MaxDouble;
if y >= 0 then y := Q else y := -MaxDouble;
end else
begin
x := x/zz;
if x > Q then x := MaxDouble
else if x < -Q then x := -MaxDouble;
y := y/zz;
if y > Q then y := MaxDouble
else if y < -Q then y := -MaxDouble
end;
end;
//------------------------------------------------------------------------------
function GetProjectionMatrix(const srcPts, dstPts: TPathD): TMatrixD;
var
dstMat: TMatrixD;
begin
if (length(srcPts) <> 4) or (length(dstPts) <> 4) then
begin
Result := IdentityMatrix;
Exit;
end;
Result := BasisToPoints(srcPts[0].X, srcPts[0].Y,
srcPts[1].X, srcPts[1].Y, srcPts[2].X, srcPts[2].Y, srcPts[3].X, srcPts[3].Y);
dstMat := BasisToPoints(dstPts[0].X, dstPts[0].Y,
dstPts[1].X, dstPts[1].Y, dstPts[2].X, dstPts[2].Y, dstPts[3].X, dstPts[3].Y);
MatrixMultiply(Result, MatrixAdjugate(dstMat));
end;
//------------------------------------------------------------------------------
function ProjectiveTransform(img: TImage32;
const srcPts, dstPts: TPathD; const margins: TRect): Boolean;
begin
Result := ProjectiveTransform(img, img, srcPts, dstPts, margins);
end;
//------------------------------------------------------------------------------
function ProjectiveTransform(img, targetImg: TImage32;
const srcPts, dstPts: TPathD; const margins: TRect): Boolean;
var
w,h,i,j: integer;
x,y: double;
xLimLo, yLimLo, xLimHi, yLimHi: double;
rec: TRect;
dstPts2: TPathD;
mat: TMatrixD;
tmp: TArrayOfColor32;
pc: PColor32;
resampler: TResamplerFunction;
begin
//https://math.stackexchange.com/a/339033/384709
if targetImg.Resampler = 0 then
resampler := nil else
resampler := GetResampler(targetImg.Resampler);
Result := Assigned(resampler) and not img.IsEmpty and
(Length(dstPts) = 4) and IsPathConvex(dstPts);
if not Result then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
rec := GetBounds(dstPts);
dec(rec.Left, margins.Left);
dec(rec.Top, margins.Top);
inc(rec.Right, margins.Right);
inc(rec.Bottom, margins.Bottom);
dstPts2 := TranslatePath(dstPts, -rec.Left, -rec.Top);
xLimLo := -0.5;
xLimHi := img.Width + 0.5;
yLimLo := -0.5;
yLimHi := img.Height + 0.5;
mat := GetProjectionMatrix(srcPts, dstPts2);
RectWidthHeight(rec, w, h);
NewColor32Array(tmp, w * h, True);
pc := @tmp[0];
for i := 0 to h -1 do
for j := 0 to w -1 do
begin
x := j; y := i;
GetSrcCoords(mat, x, y);
if (x <= xLimLo) or (x >= xLimHi) or (y <= yLimLo) or (y >= yLimHi) then
pc^ := clNone32
else
pc^ := resampler(img, x -0.5, y -0.5);
inc(pc);
end;
targetImg.BlockNotify;
targetImg.AssignPixelArray(tmp, w, h);
targetImg.UnblockNotify;
end;
//------------------------------------------------------------------------------
// Spline transformations
//------------------------------------------------------------------------------
function ReColor(color, newColor: TColor32): TColor32;
{$IFDEF INLINE} inline; {$ENDIF}
begin
Result := (color and $FF000000) or newColor;
end;
//------------------------------------------------------------------------------
function InterpolateSegX(const pt1, pt2: TPointD): TPathD;
var
i, x1, x2: integer;
xo,dydx: double;
begin
Result := nil;
if pt2.X > pt1.X then
begin
x1 := Ceil(pt1.X);
x2 := Ceil(pt2.X);
if x1 = x2 then Exit;
dydx := (pt2.Y - pt1.Y)/(pt2.X - pt1.X);
xo := x1 -pt1.X;
NewPointDArray(Result, x2-x1, True);
for i:= 0 to x2 - x1 -1 do
begin
Result[i].X := x1 +i;
Result[i].Y := pt1.Y + dydx * (xo +i);
end;
end else
begin
x1 := Floor(pt1.X);
x2 := Floor(pt2.X);
if x1 = x2 then Exit;
dydx := (pt2.Y - pt1.Y)/(pt2.X - pt1.X);
xo := x1 -pt1.X;
NewPointDArray(Result, x1-x2, True);
for i:= 0 to x1 - x2 -1 do
begin
Result[i].X := x1 -i;
Result[i].Y := pt1.Y + dydx * (xo -i);
end;
end;
end;
//------------------------------------------------------------------------------
function InterpolateSegY(const pt1, pt2: TPointD): TPathD;
var
i, y1,y2: integer;
yo,dxdy: double;
begin
Result := nil;
if pt2.Y > pt1.Y then
begin
y1 := Ceil(pt1.Y);
y2 := Ceil(pt2.Y);
if y1 = y2 then Exit;
dxdy := (pt2.X - pt1.X)/(pt2.Y - pt1.Y);
yo := y1 -pt1.Y;
NewPointDArray(Result, y2-y1, True);
for i:= 0 to y2 - y1 -1 do
begin
Result[i].Y := y1 +i;
Result[i].X := pt1.X + dxdy * (yo +i);
end;
end else
begin
y1 := Floor(pt1.Y);
y2 := Floor(pt2.Y);
if y1 = y2 then Exit;
dxdy := (pt2.X - pt1.X)/(pt2.Y - pt1.Y);
yo := y1 -pt1.Y;
NewPointDArray(Result, y1-y2, True);
for i:= 0 to y1 - y2 -1 do
begin
Result[i].Y := y1 -i;
Result[i].X := pt1.X + dxdy * (yo -i);
end;
end;
end;
//------------------------------------------------------------------------------
function InterpolatePathForX(const path: TPathD): TPathD;
var
i,len: integer;
tmp: TPathD;
begin
Result := nil;
len := length(path);
if len < 2 then Exit;
for i := 1 to len -1 do
begin
tmp := InterpolateSegX(path[i-1], path[i]);
ConcatPaths(Result, tmp);
end;
end;
//------------------------------------------------------------------------------
function InterpolatePathForY(const path: TPathD): TPathD;
var
i, len: integer;
tmp: TPathD;
begin
Result := nil;
len := length(path);
if len < 2 then Exit;
for i := 1 to len -1 do
begin
tmp := InterpolateSegY(path[i-1], path[i]);
ConcatPaths(Result, tmp);
end;
end;
//------------------------------------------------------------------------------
function SplineVertTransform(img: TImage32; const topSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean;
begin
Result := SplineVertTransform(img, img, topSpline, splineType, backColor, offset);
end;
//------------------------------------------------------------------------------
function SplineVertTransform(img, targetImg: TImage32; const topSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean;
var
i,j, w,h, len: integer;
x,y, yy, q: double;
yLimLo, yLimHi: double;
distances: TArrayOfDouble;
pc: PColor32;
rec: TRect;
tmp: TArrayOfColor32;
topPath: TPathD;
prevX: double;
resampler: TResamplerFunction;
backColoring, allowBackColoring: Boolean;
begin
offset := NullPoint;
if targetImg.Resampler = 0 then
resampler := nil else
resampler := GetResampler(targetImg.Resampler);
//convert the top spline control points into a flattened path
if splineType = stQuadratic then
topPath := FlattenQSpline(topSpline) else
topPath := FlattenCSpline(topSpline);
rec := GetBounds(topPath);
//return false if the spline is invalid or there's no vertical transformation
Result := Assigned(resampler) and not IsEmptyRect(rec);
if not Result then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
offset := rec.TopLeft;
topPath := InterpolatePathForX(topPath);
len := Length(topPath);
inc(rec.Bottom, img.Height);
RectWidthHeight(rec, w, h);
NewColor32Array(tmp, (w+1) * h, False);
prevX := topPath[0].X;
allowBackColoring := GetAlpha(backColor) > 2;
backColor := backColor and $00FFFFFF;
distances := GetCumulativeDistances(topPath);
q := img.Width / distances[High(distances)];
yLimLo := -0.5;
yLimHi := img.Height + 0.5;
for i := 0 to len -1 do
begin
pc := @tmp[Round(topPath[i].X)-rec.Left];
backColoring := allowBackColoring and (prevX >= topPath[i].X);
prevX := topPath[i].X;
yy := topPath[i].Y;
for j := rec.top to rec.bottom -1 do
begin
x := Distances[i]*q;
y := j - yy;
if (y < yLimLo) or (y > yLimHi) then
// do nothing !
else if backColoring then
pc^ := BlendToAlpha(pc^, ReColor(resampler(img, x -0.5, y -0.5), backColor))
else
pc^ := BlendToAlpha(pc^, resampler(img, x -0.5, y -0.5));
inc(pc, w);
end;
end;
// tmp was creates with "(w+1)*h". We take advantage of the
// memory manager's inplace shrink.
SetLength(tmp, w * h);
targetImg.AssignPixelArray(tmp, w, h);
end;
//------------------------------------------------------------------------------
function SplineHorzTransform(img: TImage32; const leftSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean;
begin
Result := SplineHorzTransform(img, img, leftSpline, splineType, backColor, offset);
end;
//------------------------------------------------------------------------------
function SplineHorzTransform(img, targetImg: TImage32; const leftSpline: TPathD;
splineType: TSplineType; backColor: TColor32; out offset: TPoint): Boolean;
var
i,j, len, w,h: integer;
x,y, xx, q, prevY: double;
xLimLo, xLimHi: double;
leftPath: TPathD;
distances: TArrayOfDouble;
rec: TRect;
pc: PColor32;
tmp: TArrayOfColor32;
backColoring, allowBackColoring: Boolean;
resampler: TResamplerFunction;
begin
offset := NullPoint;
if targetImg.Resampler = 0 then
resampler := nil else
resampler := GetResampler(targetImg.Resampler);
//convert the left spline control points into a flattened path
if splineType = stQuadratic then
leftPath := FlattenQSpline(leftSpline) else
leftPath := FlattenCSpline(leftSpline);
rec := GetBounds(leftPath);
//return false if the spline is invalid or there's no horizontal transformation
Result := Assigned(resampler) and not IsEmptyRect(rec);
if not Result then
begin
if targetImg <> img then targetImg.Assign(img);
Exit;
end;
offset := rec.TopLeft;
leftPath := InterpolatePathForY(leftPath);
len := Length(leftPath);
inc(rec.Right, img.Width);
RectWidthHeight(rec, w, h);
NewColor32Array(tmp, w * (h+1), False);
prevY := leftPath[0].Y;
allowBackColoring := GetAlpha(backColor) > 2;
backColor := backColor and $00FFFFFF;
distances := GetCumulativeDistances(leftPath);
q := img.Height / distances[High(distances)];;
xLimLo := -0.5;
xLimHi := img.Width + 0.5;
for i := 0 to len -1 do
begin
pc := @tmp[Round(leftPath[i].Y - rec.Top) * w];
backColoring := allowBackColoring and (prevY >= leftPath[i].Y);
prevY := leftPath[i].Y;
xx := leftPath[i].X;
y := Distances[i]*q;
for j := rec.left to rec.right -1 do
begin
x := j - xx;
if (x < xLimLo) or (x > xLimHi) then
// do nothing !
else if backColoring then
pc^ := BlendToAlpha(pc^, ReColor(resampler(img, x -0.5, y -0.5), backColor))
else
pc^ := BlendToAlpha(pc^, resampler(img, x -0.5, y -0.5));
inc(pc);
end;
end;
// tmp was creates with "w*(h+1)". We take advantage of the
// memory manager's inplace shrink.
SetLength(tmp, w * h);
targetImg.AssignPixelArray(tmp, w, h);
end;
//------------------------------------------------------------------------------
// Miscellaneous WeightedColor function
//------------------------------------------------------------------------------
function LimitByte(val: Cardinal): byte; {$IFDEF INLINE} inline; {$ENDIF}
begin
if val > 255 then result := 255
else result := val;
end;
//------------------------------------------------------------------------------
// TWeightedColor
//------------------------------------------------------------------------------
procedure TWeightedColor.Reset;
{$IFDEF CPUX64}
var
Zero: Int64;
{$ENDIF CPUX64}
begin
{$IFDEF CPUX64}
Zero := 0;
fAddCount := Zero;
fAlphaTot := Zero;
fColorTotR := Zero;
fColorTotG := Zero;
fColorTotB := Zero;
{$ELSE}
fAddCount := 0;
fAlphaTot := 0;
fColorTotR := 0;
fColorTotG := 0;
fColorTotB := 0;
{$ENDIF CPUX64}
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Reset(c: TColor32; w: Integer);
var
a: Cardinal;
begin
fAddCount := w;
a := w * Byte(c shr 24);
if a = 0 then
begin
fAlphaTot := 0;
fColorTotB := 0;
fColorTotG := 0;
fColorTotR := 0;
end else
begin
fAlphaTot := a;
fColorTotB := (a * Byte(c));
fColorTotG := (a * Byte(c shr 8));
fColorTotR := (a * Byte(c shr 16));
end;
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.AddWeight(w: Integer);
begin
inc(fAddCount, w);
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Add(c: TColor32; w: Integer);
var
a: Int64;
begin
inc(fAddCount, w);
a := Byte(c shr 24);
if a <> 0 then
begin
a := a * Cardinal(w);
inc(fAlphaTot, a);
inc(fColorTotB, (a * Byte(c)));
inc(fColorTotG, (a * Byte(c shr 8)));
inc(fColorTotR, (a * Byte(c shr 16)));
end;
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Add(c: TColor32);
// Optimized for w=1
var
a: Int64;
begin
inc(fAddCount);
a := Byte(c shr 24);
if a = 0 then Exit;
inc(fAlphaTot, a);
inc(fColorTotB, (a * Byte(c)));
inc(fColorTotG, (a * Byte(c shr 8)));
inc(fColorTotR, (a * Byte(c shr 16)));
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Add(const other: TWeightedColor);
begin
inc(fAddCount, other.fAddCount);
inc(fAlphaTot, other.fAlphaTot);
inc(fColorTotR, other.fColorTotR);
inc(fColorTotG, other.fColorTotG);
inc(fColorTotB, other.fColorTotB);
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Subtract(c: TColor32; w: Integer);
var
a: Int64;
begin
dec(fAddCount, w);
a := w * Byte(c shr 24);
if a = 0 then Exit;
dec(fAlphaTot, a);
dec(fColorTotB, (a * Byte(c)));
dec(fColorTotG, (a * Byte(c shr 8)));
dec(fColorTotR, (a * Byte(c shr 16)));
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Subtract(c: TColor32);
// Optimized for w=1
var
a: Int64;
begin
dec(fAddCount);
a := Byte(c shr 24);
if a = 0 then Exit;
dec(fAlphaTot, a);
dec(fColorTotB, (a * Byte(c)));
dec(fColorTotG, (a * Byte(c shr 8)));
dec(fColorTotR, (a * Byte(c shr 16)));
end;
//------------------------------------------------------------------------------
procedure TWeightedColor.Subtract(const other: TWeightedColor);
begin
dec(fAddCount, other.fAddCount);
dec(fAlphaTot, other.fAlphaTot);
dec(fColorTotR, other.fColorTotR);
dec(fColorTotG, other.fColorTotG);
dec(fColorTotB, other.fColorTotB);
end;
//------------------------------------------------------------------------------
function TWeightedColor.AddSubtract(addC, subC: TColor32): Boolean;
var
a: Int64;
begin
// add+subtract => fAddCount stays the same
// skip identical colors
Result := False;
if addC = subC then Exit;
a := Byte(addC shr 24);
if a > 0 then
begin
inc(fAlphaTot, a);
inc(fColorTotB, (a * Byte(addC)));
inc(fColorTotG, (a * Byte(addC shr 8)));
inc(fColorTotR, (a * Byte(addC shr 16)));
Result := True;
end;
a := Byte(subC shr 24);
if a > 0 then
begin
dec(fAlphaTot, a);
dec(fColorTotB, (a * Byte(subC)));
dec(fColorTotG, (a * Byte(subC shr 8)));
dec(fColorTotR, (a * Byte(subC shr 16)));
Result := True;
end;
end;
//------------------------------------------------------------------------------
function TWeightedColor.AddNoneSubtract(c: TColor32): Boolean;
var
a: Int64;
begin
// add+subtract => fAddCount stays the same
a := Byte(c shr 24);
if a > 0 then
begin
dec(fAlphaTot, a);
dec(fColorTotB, (a * Byte(c)));
dec(fColorTotG, (a * Byte(c shr 8)));
dec(fColorTotR, (a * Byte(c shr 16)));
Result := True;
end
else
Result := False;
end;
//------------------------------------------------------------------------------
function TWeightedColor.GetColor: TColor32;
var
oneDivAlphaTot: double;
alpha: Integer;
begin
result := clNone32;
if (fAlphaTot <= 0) or (fAddCount <= 0) then
Exit;
{$IFDEF CPUX86}
if fAlphaTot and $FFFFFFFF80000000 = 0 then // small, so can avoid _lldiv call
alpha := (Cardinal(fAlphaTot) + (Cardinal(fAddCount) shr 1)) div
Cardinal(fAddCount)
else
{$ENDIF CPUX86}
alpha := (fAlphaTot + (Cardinal(fAddCount) shr 1)) div Cardinal(fAddCount);
result := TColor32(Min(255, alpha)) shl 24;
// alpha weighting has been applied to color channels, so div by fAlphaTot
if fAlphaTot < DivOneByXTableSize then // use precalculated 1/X values
oneDivAlphaTot := DivOneByXTable[fAlphaTot] else
oneDivAlphaTot := 1/(fAlphaTot);
// 1. Skip zero calculations.
// 2. LimitByte(Integer): Values can't be less than 0, so don't use ClampByte.
// 3. x86: Round expects the value in the st(0)/xmm1 FPU register.
// Thus we need to do the calculation and Round call in one expression.
// Otherwise the compiler will use a temporary double variable on
// the stack that will cause unnecessary store and load operations.
if fColorTotB > 0 then
result := result or LimitByte(System.Round(fColorTotB * oneDivAlphaTot));
if fColorTotG > 0 then
result := result or LimitByte(System.Round(fColorTotG * oneDivAlphaTot)) shl 8;
if fColorTotR > 0 then
result := result or LimitByte(System.Round(fColorTotR * oneDivAlphaTot)) shl 16;
end;
//------------------------------------------------------------------------------
// Initialization
//------------------------------------------------------------------------------
procedure MakeDivOneByXTable;
var
i: Integer;
begin
DivOneByXTable[0] := 0; // NaN
for i := 1 to High(DivOneByXTable) do
DivOneByXTable[i] := 1/i;
end;
initialization
MakeDivOneByXTable;
end.
Reference in a new issue View git blame Copy permalink