testinterpolation2dunit.pas

maths lib with source

unit testinterpolation2dunit;
interface
uses Math, Ap, Sysutils, spline3, spline2d;

function Test2DInterpolation(Silent : Boolean):Boolean;
function testinterpolation2dunit_test_silent():Boolean;
function testinterpolation2dunit_test():Boolean;

implementation

procedure LConst(const C : TReal1DArray;
     const LX : TReal1DArray;
     const LY : TReal1DArray;
     M : Integer;
     N : Integer;
     LStep : Double;
     var LC : Double;
     var LCX : Double;
     var LCY : Double;
     var LCXY : Double);forward;
procedure TwoDNumDer(const C : TReal1DArray;
     X : Double;
     Y : Double;
     H : Double;
     var F : Double;
     var FX : Double;
     var FY : Double;
     var FXY : Double);forward;
function TestUnpack(const C : TReal1DArray;
     const LX : TReal1DArray;
     const LY : TReal1DArray):Boolean;forward;
function TestLinTrans(const C : TReal1DArray;
     AX : Double;
     BX : Double;
     AY : Double;
     BY : Double):Boolean;forward;


function Test2DInterpolation(Silent : Boolean):Boolean;
var
    WasErrors : Boolean;
    BLErrors : Boolean;
    BCErrors : Boolean;
    DSErrors : Boolean;
    CPErrors : Boolean;
    UPErrors : Boolean;
    LTErrors : Boolean;
    SYErrors : Boolean;
    RLErrors : Boolean;
    RCErrors : Boolean;
    Pass : Integer;
    PassCount : Integer;
    JobType : Integer;
    LStep : Double;
    H : Double;
    X : TReal1DArray;
    Y : TReal1DArray;
    C : TReal1DArray;
    C2 : TReal1DArray;
    LX : TReal1DArray;
    LY : TReal1DArray;
    F : TReal2DArray;
    FR : TReal2DArray;
    FT : TReal2DArray;
    AX : Double;
    AY : Double;
    BX : Double;
    BY : Double;
    I : Integer;
    J : Integer;
    K : Integer;
    N : Integer;
    M : Integer;
    N2 : Integer;
    M2 : Integer;
    Err : Double;
    T : Double;
    T1 : Double;
    T2 : Double;
    L1 : Double;
    L1X : Double;
    L1Y : Double;
    L1XY : Double;
    L2 : Double;
    L2X : Double;
    L2Y : Double;
    L2XY : Double;
    FM : Double;
    F1 : Double;
    F2 : Double;
    F3 : Double;
    F4 : Double;
    V1 : Double;
    V1X : Double;
    V1Y : Double;
    V1XY : Double;
    V2 : Double;
    V2X : Double;
    V2Y : Double;
    V2XY : Double;
    MF : Double;
begin
    WasErrors := False;
    PassCount := 10;
    H := 0.00001;
    LStep := 0.001;
    BLErrors := False;
    BCErrors := False;
    DSErrors := False;
    CPErrors := False;
    UPErrors := False;
    LTErrors := False;
    SYErrors := False;
    RLErrors := False;
    RCErrors := False;
    
    //
    // Test: bilinear, bicubic
    //
    N:=2;
    while N<=7 do
    begin
        M:=2;
        while M<=7 do
        begin
            SetLength(X, N-1+1);
            SetLength(Y, M-1+1);
            SetLength(LX, 2*N-2+1);
            SetLength(LY, 2*M-2+1);
            SetLength(F, M-1+1, N-1+1);
            SetLength(FT, N-1+1, M-1+1);
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Prepare task:
                // * X and Y stores grid
                // * F stores function values
                // * LX and LY stores twice dense grid (for Lipschitz testing)
                //
                AX := -1-RandomReal;
                BX := +1+RandomReal;
                AY := -1-RandomReal;
                BY := +1+RandomReal;
                J:=0;
                while J<=N-1 do
                begin
                    X[J] := 0.5*(BX+AX)-0.5*(BX-AX)*Cos(PI*(2*J+1)/(2*n));
                    if J=0 then
                    begin
                        X[J] := AX;
                    end;
                    if J=N-1 then
                    begin
                        X[J] := BX;
                    end;
                    LX[2*J] := X[J];
                    if J>0 then
                    begin
                        LX[2*J-1] := 0.5*(X[J]+X[J-1]);
                    end;
                    Inc(J);
                end;
                J:=0;
                while J<=N-1 do
                begin
                    K := RandomInteger(N);
                    if K<>J then
                    begin
                        T := X[J];
                        X[J] := X[K];
                        X[K] := T;
                    end;
                    Inc(J);
                end;
                I:=0;
                while I<=M-1 do
                begin
                    Y[I] := 0.5*(BY+AY)-0.5*(BY-AY)*Cos(PI*(2*i+1)/(2*M));
                    if I=0 then
                    begin
                        Y[I] := AY;
                    end;
                    if I=M-1 then
                    begin
                        Y[I] := BY;
                    end;
                    LY[2*I] := Y[I];
                    if I>0 then
                    begin
                        LY[2*I-1] := 0.5*(Y[I]+Y[I-1]);
                    end;
                    Inc(I);
                end;
                I:=0;
                while I<=M-1 do
                begin
                    K := RandomInteger(M);
                    if K<>I then
                    begin
                        T := Y[I];
                        Y[I] := Y[K];
                        Y[K] := T;
                    end;
                    Inc(I);
                end;
                I:=0;
                while I<=M-1 do
                begin
                    J:=0;
                    while J<=N-1 do
                    begin
                        F[I,J] := Exp(0.6*X[J])-Exp(-0.3*Y[I]+0.08*X[J])+2*Cos(Pi*(X[J]+1.2*Y[I]))+0.1*Cos(20*X[J]+15*Y[I]);
                        Inc(J);
                    end;
                    Inc(I);
                end;
                
                //
                // Test bilinear interpolation:
                // * interpolation at the nodes
                // * linearity
                // * continuity
                // * differentiation in the inner points
                //
                BuildBilinearSpline(X, Y, F, M, N, C);
                Err := 0;
                I:=0;
                while I<=M-1 do
                begin
                    J:=0;
                    while J<=N-1 do
                    begin
                        Err := Max(Err, AbsReal(F[I,J]-SplineInterpolation2D(C, X[J], Y[I])));
                        Inc(J);
                    end;
                    Inc(I);
                end;
                BLErrors := BLErrors or (Err>10000*MachineEpsilon);
                Err := 0;
                I:=0;
                while I<=M-2 do
                begin
                    J:=0;
                    while J<=N-2 do
                    begin
                        
                        //
                        // Test for linearity between grid points
                        // (test point - geometric center of the cell)
                        //
                        FM := SplineInterpolation2D(C, LX[2*J+1], LY[2*I+1]);
                        F1 := SplineInterpolation2D(C, LX[2*J], LY[2*I]);
                        F2 := SplineInterpolation2D(C, LX[2*J+2], LY[2*I]);
                        F3 := SplineInterpolation2D(C, LX[2*J+2], LY[2*I+2]);
                        F4 := SplineInterpolation2D(C, LX[2*J], LY[2*I+2]);
                        Err := Max(Err, AbsReal(0.25*(F1+F2+F3+F4)-FM));
                        Inc(J);
                    end;
                    Inc(I);
                end;
                BLErrors := BLErrors or (Err>10000*MachineEpsilon);
                LConst(C, LX, LY, M, N, LStep, L1, L1X, L1Y, L1XY);
                LConst(C, LX, LY, M, N, LStep/3, L2, L2X, L2Y, L2XY);
                BLErrors := BLErrors or (L2/L1>1.2);
                Err := 0;
                I:=0;
                while I<=M-2 do
                begin
                    J:=0;
                    while J<=N-2 do
                    begin
                        SplineDifferentiation2D(C, LX[2*J+1], LY[2*I+1], V1, V1X, V1Y, V1XY);
                        TwoDNumDer(C, LX[2*J+1], LY[2*I+1], H, V2, V2X, V2Y, V2XY);
                        Err := Max(Err, AbsReal(V1-V2));
                        Err := Max(Err, AbsReal(V1X-V2X));
                        Err := Max(Err, AbsReal(V1Y-V2Y));
                        Err := Max(Err, AbsReal(V1XY-V2XY));
                        Inc(J);
                    end;
                    Inc(I);
                end;
                DSErrors := DSErrors or (Err>1.0E-3);
                UPErrors := UPErrors or  not TestUnpack(C, LX, LY);
                LTErrors := LTErrors or  not TestLinTrans(C, AX, BX, AY, BY);
                
                //
                // Test bicubic interpolation.
                // * interpolation at the nodes
                // * smoothness
                // * differentiation
                //
                BuildBicubicSpline(X, Y, F, M, N, C);
                Err := 0;
                I:=0;
                while I<=M-1 do
                begin
                    J:=0;
                    while J<=N-1 do
                    begin
                        Err := Max(Err, AbsReal(F[I,J]-SplineInterpolation2D(C, X[J], Y[I])));
                        Inc(J);
                    end;
                    Inc(I);
                end;
                BCErrors := BCErrors or (Err>10000*MachineEpsilon);
                LConst(C, LX, LY, M, N, LStep, L1, L1X, L1Y, L1XY);
                LConst(C, LX, LY, M, N, LStep/3, L2, L2X, L2Y, L2XY);
                BCErrors := BCErrors or (L2/L1>1.2);
                BCErrors := BCErrors or (L2X/L1X>1.2);
                BCErrors := BCErrors or (L2Y/L1Y>1.2);
                if (L2XY>0.01) and (L1XY>0.01) then
                begin
                    
                    //
                    // Cross-derivative continuity is tested only when
                    // bigger than 0.01. When the task size is too
                    // small, the d2F/dXdY is nearly zero and Lipschitz
                    // constant ratio is ill-conditioned.
                    //
                    BCErrors := BCErrors or (L2XY/L1XY>1.2);
                end;
                Err := 0;
                I:=0;
                while I<=2*M-2 do
                begin
                    J:=0;
                    while J<=2*N-2 do
                    begin
                        SplineDifferentiation2D(C, LX[J], LY[I], V1, V1X, V1Y, V1XY);
                        TwoDNumDer(C, LX[J], LY[I], H, V2, V2X, V2Y, V2XY);
                        Err := Max(Err, AbsReal(V1-V2));
                        Err := Max(Err, AbsReal(V1X-V2X));
                        Err := Max(Err, AbsReal(V1Y-V2Y));
                        Err := Max(Err, AbsReal(V1XY-V2XY));
                        Inc(J);
                    end;
                    Inc(I);
                end;
                DSErrors := DSErrors or (Err>1.0E-3);
                UPErrors := UPErrors or  not TestUnpack(C, LX, LY);
                LTErrors := LTErrors or  not TestLinTrans(C, AX, BX, AY, BY);
                
                //
                // Copy test
                //
                if RandomReal>0.5 then
                begin
                    BuildBicubicSpline(X, Y, F, M, N, C);
                end
                else
                begin
                    BuildBilinearSpline(X, Y, F, M, N, C);
                end;
                Spline2DCopy(C, C2);
                Err := 0;
                I:=1;
                while I<=5 do
                begin
                    T1 := AX+(BX-AX)*RandomReal;
                    T2 := AY+(BY-AY)*RandomReal;
                    Err := Max(Err, AbsReal(SplineInterpolation2D(C, T1, T2)-SplineInterpolation2D(C2, T1, T2)));
                    Inc(I);
                end;
                CPErrors := CPErrors or (Err>10000*MachineEpsilon);
                
                //
                // Special symmetry test
                //
                Err := 0;
                JobType:=0;
                while JobType<=1 do
                begin
                    
                    //
                    // Prepare
                    //
                    I:=0;
                    while I<=M-1 do
                    begin
                        J:=0;
                        while J<=N-1 do
                        begin
                            FT[J,I] := F[I,J];
                            Inc(J);
                        end;
                        Inc(I);
                    end;
                    if JobType=0 then
                    begin
                        BuildBilinearSpline(X, Y, F, M, N, C);
                        BuildBilinearSpline(Y, X, FT, N, M, C2);
                    end
                    else
                    begin
                        BuildBicubicSpline(X, Y, F, M, N, C);
                        BuildBicubicSpline(Y, X, FT, N, M, C2);
                    end;
                    
                    //
                    // Test
                    //
                    I:=1;
                    while I<=10 do
                    begin
                        T1 := AX+(BX-AX)*RandomReal;
                        T2 := AY+(BY-AY)*RandomReal;
                        Err := Max(Err, AbsReal(SplineInterpolation2D(C, T1, T2)-SplineInterpolation2D(C2, T2, T1)));
                        Inc(I);
                    end;
                    Inc(JobType);
                end;
                SYErrors := SYErrors or (Err>10000*MachineEpsilon);
                Inc(Pass);
            end;
            Inc(M);
        end;
        Inc(N);
    end;
    
    //
    // Test resample
    //
    M:=2;
    while M<=6 do
    begin
        N:=2;
        while N<=6 do
        begin
            SetLength(F, M-1+1, N-1+1);
            SetLength(X, N-1+1);
            SetLength(Y, M-1+1);
            J:=0;
            while J<=N-1 do
            begin
                X[J] := J/(N-1);
                Inc(J);
            end;
            I:=0;
            while I<=M-1 do
            begin
                Y[I] := I/(M-1);
                Inc(I);
            end;
            I:=0;
            while I<=M-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    F[I,J] := Exp(0.6*X[J])-Exp(-0.3*Y[I]+0.08*X[J])+2*Cos(Pi*(X[J]+1.2*Y[I]))+0.1*Cos(20*X[J]+15*Y[I]);
                    Inc(J);
                end;
                Inc(I);
            end;
            M2:=2;
            while M2<=6 do
            begin
                N2:=2;
                while N2<=6 do
                begin
                    Pass:=1;
                    while Pass<=PassCount do
                    begin
                        JobType:=0;
                        while JobType<=1 do
                        begin
                            if JobType=0 then
                            begin
                                BilinearResampleCartesian(F, M, N, FR, M2, N2);