mirror of
https://gitlab.com/freepascal.org/lazarus/lazarus.git
synced 2025-06-06 09:58:16 +02:00
479 lines
11 KiB
ObjectPascal
479 lines
11 KiB
ObjectPascal
//----------------------------------------------------------------------------
|
|
// Anti-Grain Geometry - Version 2.4 (Public License)
|
|
// Copyright (C) 2002-2005 Maxim Shemanarev (http://www.antigrain.com)
|
|
//
|
|
// Anti-Grain Geometry - Version 2.4 Release Milano 3 (AggPas 2.4 RM3)
|
|
// Pascal Port By: Milan Marusinec alias Milano
|
|
// milan@marusinec.sk
|
|
// http://www.aggpas.org
|
|
// Copyright (c) 2005-2006
|
|
//
|
|
// Permission to copy, use, modify, sell and distribute this software
|
|
// is granted provided this copyright notice appears in all copies.
|
|
// This software is provided "as is" without express or implied
|
|
// warranty, and with no claim as to its suitability for any purpose.
|
|
//
|
|
//----------------------------------------------------------------------------
|
|
// Contact: mcseem@antigrain.com
|
|
// mcseemagg@yahoo.com
|
|
// http://www.antigrain.com
|
|
//
|
|
//----------------------------------------------------------------------------
|
|
//
|
|
// class bspline
|
|
//
|
|
// [Pascal Port History] -----------------------------------------------------
|
|
//
|
|
// 23.06.2006-Milano: ptrcomp adjustments
|
|
// 19.01.2006-Milano: Complete Unit Port
|
|
// 18.01.2006-Milano: Unit port establishment
|
|
//
|
|
{ agg_bspline.pas }
|
|
unit
|
|
agg_bspline ;
|
|
|
|
INTERFACE
|
|
|
|
{$I agg_mode.inc }
|
|
|
|
uses
|
|
agg_basics ;
|
|
|
|
{ TYPES DEFINITION }
|
|
type
|
|
//----------------------------------------------------------------bspline
|
|
// A very simple class of Bi-cubic Spline interpolation.
|
|
// First call init(num, x[], y[]) where num - number of source points,
|
|
// x, y - arrays of X and Y values respectively. Here Y must be a function
|
|
// of X. It means that all the X-coordinates must be arranged in the ascending
|
|
// order.
|
|
// Then call get(x) that calculates a value Y for the respective X.
|
|
// The class supports extrapolation, i.e. you can call get(x) where x is
|
|
// outside the given with init() X-range. Extrapolation is a simple linear
|
|
// function.
|
|
bspline = object
|
|
m_max ,
|
|
m_num : int;
|
|
|
|
m_x ,
|
|
m_y ,
|
|
m_am : double_ptr;
|
|
|
|
m_last_idx : int;
|
|
|
|
constructor Construct; overload;
|
|
constructor Construct(num : int ); overload;
|
|
constructor Construct(num : int; x ,y : double_ptr ); overload;
|
|
destructor Destruct;
|
|
|
|
procedure init(max : int ); overload;
|
|
procedure init(num : int; x ,y : double_ptr ); overload;
|
|
|
|
procedure add_point(x ,y : double );
|
|
procedure prepare;
|
|
|
|
function get (x : double ) : double;
|
|
function get_stateful(x : double ) : double;
|
|
|
|
procedure bsearch(n : int; x : double_ptr; x0 : double; i : int_ptr );
|
|
|
|
function extrapolation_left (x : double ) : double;
|
|
function extrapolation_right(x : double ) : double;
|
|
|
|
function interpolation(x : double; i : int ) : double;
|
|
|
|
end;
|
|
|
|
{ GLOBAL PROCEDURES }
|
|
|
|
|
|
IMPLEMENTATION
|
|
{ LOCAL VARIABLES & CONSTANTS }
|
|
{ UNIT IMPLEMENTATION }
|
|
{ CONSTRUCT }
|
|
constructor bspline.Construct;
|
|
begin
|
|
m_max:=0;
|
|
m_num:=0;
|
|
|
|
m_x :=NIL;
|
|
m_y :=NIL;
|
|
m_am:=NIL;
|
|
|
|
m_last_idx:=-1;
|
|
|
|
end;
|
|
|
|
{ CONSTRUCT }
|
|
constructor bspline.Construct(num : int );
|
|
begin
|
|
Construct;
|
|
|
|
init(num );
|
|
|
|
end;
|
|
|
|
{ CONSTRUCT }
|
|
constructor bspline.Construct(num : int; x ,y : double_ptr );
|
|
begin
|
|
Construct;
|
|
|
|
init(num ,x ,y );
|
|
|
|
end;
|
|
|
|
{ DESTRUCT }
|
|
destructor bspline.Destruct;
|
|
begin
|
|
agg_freemem(pointer(m_am ) ,m_max * 3 * sizeof(double ) );
|
|
|
|
end;
|
|
|
|
{ INIT }
|
|
procedure bspline.init(max : int );
|
|
begin
|
|
if (max > 2 ) and
|
|
(max > m_max ) then
|
|
begin
|
|
agg_freemem(pointer(m_am ) ,m_max * 3 * sizeof(double ) );
|
|
agg_getmem (pointer(m_am ) ,max * 3 * sizeof(double ) );
|
|
|
|
m_max:=max;
|
|
|
|
m_x:=double_ptr(ptrcomp(m_am ) + m_max * sizeof(double ) );
|
|
m_y:=double_ptr(ptrcomp(m_am ) + m_max * 2 * sizeof(double ) );
|
|
|
|
end;
|
|
|
|
m_num :=0;
|
|
m_last_idx:=-1;
|
|
|
|
end;
|
|
|
|
{ INIT }
|
|
procedure bspline.init(num : int; x ,y : double_ptr );
|
|
var
|
|
i : int;
|
|
|
|
begin
|
|
if num > 2 then
|
|
begin
|
|
init(num );
|
|
|
|
for i:=0 to num - 1 do
|
|
begin
|
|
add_point(x^ ,y^ );
|
|
|
|
inc(ptrcomp(x ) ,sizeof(double ) );
|
|
inc(ptrcomp(y ) ,sizeof(double ) );
|
|
|
|
end;
|
|
|
|
prepare;
|
|
|
|
end;
|
|
|
|
m_last_idx:=-1;
|
|
|
|
end;
|
|
|
|
{ ADD_POINT }
|
|
procedure bspline.add_point;
|
|
begin
|
|
if m_num < m_max then
|
|
begin
|
|
double_ptr(ptrcomp(m_x ) + m_num * sizeof(double ) )^:=x;
|
|
double_ptr(ptrcomp(m_y ) + m_num * sizeof(double ) )^:=y;
|
|
|
|
inc(m_num );
|
|
|
|
end;
|
|
|
|
end;
|
|
|
|
{ PREPARE }
|
|
procedure bspline.prepare;
|
|
var
|
|
i ,k ,n1 ,sz : int;
|
|
temp ,r ,s ,al : double_ptr;
|
|
h ,p ,d ,f ,e : double;
|
|
|
|
begin
|
|
if m_num > 2 then
|
|
begin
|
|
for k:=0 to m_num - 1 do
|
|
double_ptr(ptrcomp(m_am ) + k * sizeof(double ) )^:=0;
|
|
|
|
n1:=3 * m_num;
|
|
sz:=n1;
|
|
|
|
agg_getmem(pointer(al ) ,n1 * sizeof(double ) );
|
|
|
|
temp:=al;
|
|
|
|
for k:=0 to n1 - 1 do
|
|
double_ptr(ptrcomp(temp ) + k * sizeof(double ) )^:=0;
|
|
|
|
r:=double_ptr(ptrcomp(temp ) + m_num * sizeof(double ) );
|
|
s:=double_ptr(ptrcomp(temp ) + m_num * 2 * sizeof(double ) );
|
|
|
|
n1:=m_num - 1;
|
|
d :=double_ptr(ptrcomp(m_x ) + sizeof(double ) )^ - m_x^;
|
|
e :=(double_ptr(ptrcomp(m_y ) + sizeof(double ) )^ - m_y^ ) / d;
|
|
|
|
k:=1;
|
|
|
|
while k < n1 do
|
|
begin
|
|
h:=d;
|
|
d:=double_ptr(ptrcomp(m_x ) + (k + 1 ) * sizeof(double ) )^ - double_ptr(ptrcomp(m_x ) + k * sizeof(double ) )^;
|
|
f:=e;
|
|
e:=(double_ptr(ptrcomp(m_y ) + (k + 1 ) * sizeof(double ) )^ - double_ptr(ptrcomp(m_y ) + k * sizeof(double ) )^ ) / d;
|
|
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^:=d / (d + h );
|
|
double_ptr(ptrcomp(r ) + k * sizeof(double ) )^ :=1.0 - double_ptr(ptrcomp(al ) + k * sizeof(double ) )^;
|
|
double_ptr(ptrcomp(s ) + k * sizeof(double ) )^ :=6.0 * (e - f ) / (h + d );
|
|
|
|
inc(k );
|
|
|
|
end;
|
|
|
|
k:=1;
|
|
|
|
while k < n1 do
|
|
begin
|
|
p:=1.0 / (double_ptr(ptrcomp(r ) + k * sizeof(double ) )^ * double_ptr(ptrcomp(al ) + (k - 1 ) * sizeof(double ) )^ + 2.0 );
|
|
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^:=
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^ * -p;
|
|
|
|
double_ptr(ptrcomp(s ) + k * sizeof(double ) )^ :=
|
|
(double_ptr(ptrcomp(s ) + k * sizeof(double ) )^ -
|
|
double_ptr(ptrcomp(r ) + k * sizeof(double ) )^ *
|
|
double_ptr(ptrcomp(s ) + (k - 1 ) * sizeof(double ) )^ ) * p;
|
|
|
|
inc(k );
|
|
|
|
end;
|
|
|
|
double_ptr(ptrcomp(m_am ) + n1 * sizeof(double ) )^:=0.0;
|
|
|
|
double_ptr(ptrcomp(al ) + (n1 - 1 ) * sizeof(double ) )^:=
|
|
double_ptr(ptrcomp(s ) + (n1 - 1 ) * sizeof(double ) )^;
|
|
|
|
double_ptr(ptrcomp(m_am ) + (n1 - 1 ) * sizeof(double ) )^:=
|
|
double_ptr(ptrcomp(al ) + (n1 - 1 ) * sizeof(double ) )^;
|
|
|
|
k:=n1 - 2;
|
|
i:=0;
|
|
|
|
while i < m_num - 2 do
|
|
begin
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^:=
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^ *
|
|
double_ptr(ptrcomp(al ) + (k + 1 ) * sizeof(double ) )^ +
|
|
double_ptr(ptrcomp(s ) + k * sizeof(double ) )^;
|
|
|
|
double_ptr(ptrcomp(m_am ) + k * sizeof(double ) )^:=
|
|
double_ptr(ptrcomp(al ) + k * sizeof(double ) )^;
|
|
|
|
inc(i );
|
|
dec(k );
|
|
|
|
end;
|
|
|
|
agg_freemem(pointer(al ) ,sz * sizeof(double ) );
|
|
|
|
end;
|
|
|
|
m_last_idx:=-1;
|
|
|
|
end;
|
|
|
|
{ GET }
|
|
function bspline.get;
|
|
var
|
|
i : int;
|
|
|
|
begin
|
|
if m_num > 2 then
|
|
begin
|
|
// Extrapolation on the left
|
|
if x < m_x^ then
|
|
begin
|
|
result:=extrapolation_left(x );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
// Extrapolation on the right
|
|
if x >= double_ptr(ptrcomp(m_x ) + (m_num - 1 ) * sizeof(double ) )^ then
|
|
begin
|
|
result:=extrapolation_right(x );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
// Interpolation
|
|
bsearch(m_num ,m_x ,x ,@i );
|
|
|
|
result:=interpolation(x ,i );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
result:=0.0;
|
|
|
|
end;
|
|
|
|
{ GET_STATEFUL }
|
|
function bspline.get_stateful;
|
|
begin
|
|
if m_num > 2 then
|
|
begin
|
|
// Extrapolation on the left
|
|
if x < m_x^ then
|
|
begin
|
|
result:=extrapolation_left(x );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
// Extrapolation on the right
|
|
if x >= double_ptr(ptrcomp(m_x ) + (m_num - 1 ) * sizeof(double ) )^ then
|
|
begin
|
|
result:=extrapolation_right(x );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
if m_last_idx >= 0 then
|
|
begin
|
|
// Check if x is not in current range
|
|
if (x < double_ptr(ptrcomp(m_x ) + m_last_idx * sizeof(double ) )^ ) or
|
|
(x > double_ptr(ptrcomp(m_x ) + (m_last_idx + 1 ) * sizeof(double ) )^ ) then
|
|
// Check if x between next points (most probably)
|
|
if (m_last_idx < m_num - 2 ) and
|
|
(x >= double_ptr(ptrcomp(m_x ) + (m_last_idx + 1 ) * sizeof(double ) )^ ) and
|
|
(x <= double_ptr(ptrcomp(m_x ) + (m_last_idx + 2 ) * sizeof(double ) )^ ) then
|
|
inc(m_last_idx )
|
|
else
|
|
if (m_last_idx > 0 ) and
|
|
(x >= double_ptr(ptrcomp(m_x ) + (m_last_idx - 1 ) * sizeof(double ) )^ ) and
|
|
(x <= double_ptr(ptrcomp(m_x ) + m_last_idx * sizeof(double ) )^ ) then
|
|
// x is between pevious points
|
|
dec(m_last_idx )
|
|
else
|
|
// Else perform full search
|
|
bsearch(m_num ,m_x ,x ,@m_last_idx );
|
|
|
|
result:=interpolation(x ,m_last_idx );
|
|
|
|
exit;
|
|
|
|
end
|
|
else
|
|
begin
|
|
// Interpolation
|
|
bsearch(m_num ,m_x ,x ,@m_last_idx );
|
|
|
|
result:=interpolation(x ,m_last_idx );
|
|
|
|
exit;
|
|
|
|
end;
|
|
|
|
end;
|
|
|
|
result:=0.0;
|
|
|
|
end;
|
|
|
|
{ BSEARCH }
|
|
procedure bspline.bsearch;
|
|
var
|
|
j ,k : int;
|
|
|
|
begin
|
|
j :=n - 1;
|
|
i^:=0;
|
|
|
|
while j - i^ > 1 do
|
|
begin
|
|
k:=shr_int32(i^ + j ,1 );
|
|
|
|
if x0 < double_ptr(ptrcomp(x ) + k * sizeof(double ) )^ then
|
|
j:=k
|
|
else
|
|
i^:=k;
|
|
|
|
end;
|
|
|
|
end;
|
|
|
|
{ EXTRAPOLATION_LEFT }
|
|
function bspline.extrapolation_left;
|
|
var
|
|
d : double;
|
|
|
|
begin
|
|
d:=double_ptr(ptrcomp(m_x ) + sizeof(double ) )^ - m_x^;
|
|
|
|
result:=
|
|
(-d * double_ptr(ptrcomp(m_am ) + sizeof(double ) )^ / 6 +
|
|
(double_ptr(ptrcomp(m_y ) + sizeof(double ) )^ - m_y^ ) / d ) *
|
|
(x - m_x^ ) + m_y^;
|
|
|
|
end;
|
|
|
|
{ EXTRAPOLATION_RIGHT }
|
|
function bspline.extrapolation_right;
|
|
var
|
|
d : double;
|
|
|
|
begin
|
|
d:=
|
|
double_ptr(ptrcomp(m_x ) + (m_num - 1 ) * sizeof(double ) )^ -
|
|
double_ptr(ptrcomp(m_x ) + (m_num - 2 ) * sizeof(double ) )^;
|
|
|
|
result:=
|
|
(d * double_ptr(ptrcomp(m_am ) + (m_num - 2 ) * sizeof(double ) )^ / 6 +
|
|
(double_ptr(ptrcomp(m_y ) + (m_num - 1 ) * sizeof(double ) )^ -
|
|
double_ptr(ptrcomp(m_y ) + (m_num - 2 ) * sizeof(double ) )^ ) / d ) *
|
|
(x - double_ptr(ptrcomp(m_x ) + (m_num - 1 ) * sizeof(double ) )^ ) +
|
|
double_ptr(ptrcomp(m_y ) + (m_num - 1 ) * sizeof(double ) )^;
|
|
|
|
end;
|
|
|
|
{ INTERPOLATION }
|
|
function bspline.interpolation;
|
|
var
|
|
j : int;
|
|
|
|
d ,h ,r ,p : double;
|
|
|
|
begin
|
|
j:=i + 1;
|
|
d:=double_ptr(ptrcomp(m_x ) + i * sizeof(double ) )^ - double_ptr(ptrcomp(m_x ) + j * sizeof(double ) )^;
|
|
h:=x - double_ptr(ptrcomp(m_x ) + j * sizeof(double ) )^;
|
|
r:=double_ptr(ptrcomp(m_x ) + i * sizeof(double ) )^ - x;
|
|
p:=d * d / 6.0;
|
|
|
|
result:=
|
|
(double_ptr(ptrcomp(m_am ) + j * sizeof(double ) )^ * r * r * r +
|
|
double_ptr(ptrcomp(m_am ) + i * sizeof(double ) )^ * h * h * h ) / 6.0 / d +
|
|
((double_ptr(ptrcomp(m_y ) + j * sizeof(double ) )^ -
|
|
double_ptr(ptrcomp(m_am ) + j * sizeof(double ) )^ * p ) * r +
|
|
(double_ptr(ptrcomp(m_y ) + i * sizeof(double ) )^ -
|
|
double_ptr(ptrcomp(m_am ) + i * sizeof(double ) )^ * p ) * h ) / d;
|
|
|
|
end;
|
|
|
|
END.
|
|
|