testinterpolation2dunit.pas

maths lib with source

                                BuildBilinearSpline(X, Y, F, M, N, C);
                            end;
                            if JobType=1 then
                            begin
                                BicubicResampleCartesian(F, M, N, FR, M2, N2);
                                BuildBicubicSpline(X, Y, F, M, N, C);
                            end;
                            Err := 0;
                            MF := 0;
                            I:=0;
                            while I<=M2-1 do
                            begin
                                J:=0;
                                while J<=N2-1 do
                                begin
                                    V1 := SplineInterpolation2D(C, J/(N2-1), I/(M2-1));
                                    V2 := FR[I,J];
                                    Err := Max(Err, AbsReal(V1-V2));
                                    MF := Max(MF, AbsReal(V1));
                                    Inc(J);
                                end;
                                Inc(I);
                            end;
                            if JobType=0 then
                            begin
                                RLErrors := RLErrors or (Err/MF>10000*MachineEpsilon);
                            end;
                            if JobType=1 then
                            begin
                                RCErrors := RCErrors or (Err/MF>10000*MachineEpsilon);
                            end;
                            Inc(JobType);
                        end;
                        Inc(Pass);
                    end;
                    Inc(N2);
                end;
                Inc(M2);
            end;
            Inc(N);
        end;
        Inc(M);
    end;
    
    //
    // report
    //
    WasErrors := BLErrors or BCErrors or DSErrors or CPErrors or UPErrors or LTErrors or SYErrors or RLErrors or RCErrors;
    if  not Silent then
    begin
        Write(Format('TESTING 2D INTERPOLATION'#13#10'',[]));
        
        //
        // Normal tests
        //
        Write(Format('BILINEAR TEST:                           ',[]));
        if BLErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('BICUBIC TEST:                            ',[]));
        if BCErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('DIFFERENTIATION TEST:                    ',[]));
        if DSErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('COPY TEST:                               ',[]));
        if CPErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('UNPACK TEST:                             ',[]));
        if UPErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('LIN.TRANS. TEST:                         ',[]));
        if LTErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('SPECIAL SYMMETRY TEST:                   ',[]));
        if SYErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('BILINEAR RESAMPLING TEST:                ',[]));
        if RLErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('BICUBIC RESAMPLING TEST:                 ',[]));
        if RCErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        
        //
        // Summary
        //
        if WasErrors then
        begin
            Write(Format('TEST FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    
    //
    // end
    //
    Result :=  not WasErrors;
end;


(*************************************************************************
Lipschitz constants for spline inself, first and second derivatives.
*************************************************************************)
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);
var
    I : Integer;
    J : Integer;
    F1 : Double;
    F2 : Double;
    F3 : Double;
    F4 : Double;
    FX1 : Double;
    FX2 : Double;
    FX3 : Double;
    FX4 : Double;
    FY1 : Double;
    FY2 : Double;
    FY3 : Double;
    FY4 : Double;
    FXY1 : Double;
    FXY2 : Double;
    FXY3 : Double;
    FXY4 : Double;
    S2LStep : Double;
begin
    LC := 0;
    LCX := 0;
    LCY := 0;
    LCXY := 0;
    S2LStep := Sqrt(2)*LStep;
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            
            //
            // Calculate
            //
            TwoDNumDer(C, LX[J]-LStep/2, LY[I]-LStep/2, LStep/4, F1, FX1, FY1, FXY1);
            TwoDNumDer(C, LX[J]+LStep/2, LY[I]-LStep/2, LStep/4, F2, FX2, FY2, FXY2);
            TwoDNumDer(C, LX[J]+LStep/2, LY[I]+LStep/2, LStep/4, F3, FX3, FY3, FXY3);
            TwoDNumDer(C, LX[J]-LStep/2, LY[I]+LStep/2, LStep/4, F4, FX4, FY4, FXY4);
            
            //
            // Lipschitz constant for the function itself
            //
            LC := Max(LC, AbsReal((F1-F2)/LStep));
            LC := Max(LC, AbsReal((F2-F3)/LStep));
            LC := Max(LC, AbsReal((F3-F4)/LStep));
            LC := Max(LC, AbsReal((F4-F1)/LStep));
            LC := Max(LC, AbsReal((F1-F3)/S2LStep));
            LC := Max(LC, AbsReal((F2-F4)/S2LStep));
            
            //
            // Lipschitz constant for the first derivative
            //
            LCX := Max(LCX, AbsReal((FX1-FX2)/LStep));
            LCX := Max(LCX, AbsReal((FX2-FX3)/LStep));
            LCX := Max(LCX, AbsReal((FX3-FX4)/LStep));
            LCX := Max(LCX, AbsReal((FX4-FX1)/LStep));
            LCX := Max(LCX, AbsReal((FX1-FX3)/S2LStep));
            LCX := Max(LCX, AbsReal((FX2-FX4)/S2LStep));
            
            //
            // Lipschitz constant for the first derivative
            //
            LCY := Max(LCY, AbsReal((FY1-FY2)/LStep));
            LCY := Max(LCY, AbsReal((FY2-FY3)/LStep));
            LCY := Max(LCY, AbsReal((FY3-FY4)/LStep));
            LCY := Max(LCY, AbsReal((FY4-FY1)/LStep));
            LCY := Max(LCY, AbsReal((FY1-FY3)/S2LStep));
            LCY := Max(LCY, AbsReal((FY2-FY4)/S2LStep));
            
            //
            // Lipschitz constant for the cross-derivative
            //
            LCXY := Max(LCXY, AbsReal((FXY1-FXY2)/LStep));
            LCXY := Max(LCXY, AbsReal((FXY2-FXY3)/LStep));
            LCXY := Max(LCXY, AbsReal((FXY3-FXY4)/LStep));
            LCXY := Max(LCXY, AbsReal((FXY4-FXY1)/LStep));
            LCXY := Max(LCXY, AbsReal((FXY1-FXY3)/S2LStep));
            LCXY := Max(LCXY, AbsReal((FXY2-FXY4)/S2LStep));
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Numerical differentiation.
*************************************************************************)
procedure TwoDNumDer(const C : TReal1DArray;
     X : Double;
     Y : Double;
     H : Double;
     var F : Double;
     var FX : Double;
     var FY : Double;
     var FXY : Double);
begin
    F := SplineInterpolation2D(C, X, Y);
    FX := (SplineInterpolation2D(C, X+H, Y)-SplineInterpolation2D(C, X-H, Y))/(2*H);
    FY := (SplineInterpolation2D(C, X, Y+H)-SplineInterpolation2D(C, X, Y-H))/(2*H);
    FXY := (SplineInterpolation2D(C, X+H, Y+H)-SplineInterpolation2D(C, X-H, Y+H)-SplineInterpolation2D(C, X+H, Y-H)+SplineInterpolation2D(C, X-H, Y-H))/Sqr(2*H);
end;


(*************************************************************************
Unpack test
*************************************************************************)
function TestUnpack(const C : TReal1DArray;
     const LX : TReal1DArray;
     const LY : TReal1DArray):Boolean;
var
    I : Integer;
    J : Integer;
    N : Integer;
    M : Integer;
    CI : Integer;
    CJ : Integer;
    P : Integer;
    Err : Double;
    TX : Double;
    TY : Double;
    V1 : Double;
    V2 : Double;
    Pass : Integer;
    PassCount : Integer;
    Tbl : TReal2DArray;
begin
    PassCount := 20;
    Err := 0;
    SplineUnpack2D(C, M, N, Tbl);
    I:=0;
    while I<=M-2 do
    begin
        J:=0;
        while J<=N-2 do
        begin
            Pass:=1;
            while Pass<=PassCount do
            begin
                P := (N-1)*I+J;
                TX := (0.001+0.999*RandomReal)*(Tbl[P,1]-Tbl[P,0]);
                TY := (0.001+0.999*RandomReal)*(Tbl[P,3]-Tbl[P,2]);
                
                //
                // Interpolation properties
                //
                V1 := 0;
                CI:=0;
                while CI<=3 do
                begin
                    CJ:=0;
                    while CJ<=3 do
                    begin
                        V1 := V1+Tbl[P,4+CI*4+CJ]*Power(TX, CI)*Power(TY, CJ);
                        Inc(CJ);
                    end;
                    Inc(CI);
                end;
                V2 := SplineInterpolation2D(C, Tbl[P,0]+TX, Tbl[P,2]+TY);
                Err := Max(Err, AbsReal(V1-V2));
                
                //
                // Grid correctness
                //
                Err := Max(Err, AbsReal(LX[2*J]-Tbl[P,0]));
                Err := Max(Err, AbsReal(LX[2*(J+1)]-Tbl[P,1]));
                Err := Max(Err, AbsReal(LY[2*I]-Tbl[P,2]));
                Err := Max(Err, AbsReal(LY[2*(I+1)]-Tbl[P,3]));
                Inc(Pass);
            end;
            Inc(J);
        end;
        Inc(I);
    end;
    Result := Err<10000*MachineEpsilon;
end;


(*************************************************************************
LinTrans test
*************************************************************************)
function TestLinTrans(const C : TReal1DArray;
     AX : Double;
     BX : Double;
     AY : Double;
     BY : Double):Boolean;
var
    Err : Double;
    A1 : Double;
    A2 : Double;
    B1 : Double;
    B2 : Double;
    TX : Double;
    TY : Double;
    VX : Double;
    VY : Double;
    V1 : Double;
    V2 : Double;
    Pass : Integer;
    PassCount : Integer;
    XJob : Integer;
    YJob : Integer;
    C2 : TReal1DArray;
begin
    PassCount := 5;
    Err := 0;
    XJob:=0;
    while XJob<=1 do
    begin
        YJob:=0;
        while YJob<=1 do
        begin
            Pass:=1;
            while Pass<=PassCount do
            begin
                
                //
                // Prepare
                //
                repeat
                    A1 := 2*RandomReal-1;
                until A1<>0;
                A1 := A1*XJob;
                B1 := 2*RandomReal-1;
                repeat
                    A2 := 2*RandomReal-1;
                until A2<>0;
                A2 := A2*YJob;
                B2 := 2*RandomReal-1;
                
                //
                // Test XY
                //
                SplineCopy(C, C2);
                Spline2DLinTransXY(C2, A1, B1, A2, B2);
                TX := AX+RandomReal*(BX-AX);
                TY := AY+RandomReal*(BY-AY);
                if XJob=0 then
                begin
                    TX := B1;
                    VX := AX+RandomReal*(BX-AX);
                end
                else
                begin
                    VX := (TX-B1)/A1;
                end;
                if YJob=0 then
                begin
                    TY := B2;
                    VY := AY+RandomReal*(BY-AY);
                end
                else
                begin
                    VY := (TY-B2)/A2;
                end;
                V1 := SplineInterpolation2D(C, TX, TY);
                V2 := SplineInterpolation2D(C2, VX, VY);
                Err := Max(Err, AbsReal(V1-V2));
                
                //
                // Test F
                //
                SplineCopy(C, C2);
                Spline2DLinTransF(C2, A1, B1);
                TX := AX+RandomReal*(BX-AX);
                TY := AY+RandomReal*(BY-AY);
                V1 := SplineInterpolation2D(C, TX, TY);
                V2 := SplineInterpolation2D(C2, TX, TY);
                Err := Max(Err, AbsReal(A1*V1+B1-V2));
                Inc(Pass);
            end;
            Inc(YJob);
        end;
        Inc(XJob);
    end;
    Result := Err<10000*MachineEpsilon;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testinterpolation2dunit_test_silent():Boolean;
begin
    Result := Test2DInterpolation(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testinterpolation2dunit_test():Boolean;
begin
    Result := Test2DInterpolation(False);
end;


end.