'=========================================================================== ' Subject: CUBIC SPLINE INTERPOLATOR Date: 11-01-99 (20:03) ' Author: Adam O'Neill Code: PBDLL ' Origin: adamaki@hotmail.com Packet: PBDLL.ABC '=========================================================================== '====================================================================== ' Splint.bas Cubic spline algorithm ' ' This algorithm fits a cubic equation to a temporal or spatial 2D ' (x,y) dataset and interpolates data in between the measured points. ' ' Although the highest possible sample density is always desired, ' usually time(cost) and/or memory are the limiting factors. ' ' Splines are also useful for digital to analog conversion, ' ie. making a smooth curve out of discrete points. ' ' This spline algorithm requires: ' 1. Noise (spurious signal) to be minimal - any extreme noise will ' introduce anomalies into the curve, and it may even become ' unstable. ' 2. Input data "should" be evenly sampled - this is for best results, ' although as shown here it seems to work OK otherwise ' ' Translated from VAX-VMS Fortran by Adam O'Neill ' Perth, Western Australia, Oct 1999 ' ' Original code from School of Physical Sciences, ' Curtin University of Technology, Perth, Western Australia ' ' Graphic window framework is from Dave Navarro's ' ftp://ftp.powerbasic.com/pub/pbdll32/GD.ZIP ' ' Notes: ' 1. There are better algorithms for this - please share them around! ' They also include a variable "tension factor" for closer fitting. ' 2. No error trapping is included ' 3. It "sometimes" misbehaves - when the interpolation interval is ' exactly half the sample interval, or greater than the minimum ' sample interval. ' 4. See how it deviates from the expected curve in areas of sparse ' data, such as towards the far right in the example data here. ' '====================================================================== $DIM ALL $COMPILE EXE $INCLUDE "win32api.inc" %AbsMinPts = 1 ' Arrays start at index 1 %AbsMaxPts = 32767 ' Absolute allowable no. of points %MaxCols = 2 ' One x column and one y column %Radius = 2 ' Radius of the data points on the graph %BorderThick = 3 ' Thickness of border %GraphThick = 1 ' Thickness of curve %XMargin = 10 ' As a percent of the window width %YMargin = 10 ' %GMargin = 2 ' As a percent of the window width '====================================================================== GLOBAL gX() AS SINGLE ' The x data GLOBAL gY() AS SINGLE ' Corresponding y data GLOBAL gXMin AS SINGLE ' Min and max of the arrays GLOBAL gXMax AS SINGLE GLOBAL gYMin AS SINGLE GLOBAL gYMax AS SINGLE GLOBAL gMinPts AS INTEGER ' Array minimum index GLOBAL gMaxPts AS INTEGER ' Array maximum index (i.e. data counter) GLOBAL gsTemp AS STRING GLOBAL gxInt AS SINGLE '====================================================================== SUB SPLINT (x() AS SINGLE,_ y() AS SINGLE,_ maxpts AS INTEGER,_ xint AS SINGLE) EXPORT '---------- 'IMPORT: x() array of x values ' y() array of corresponding y values ' maxpts number of data points in each the arrays ' xint the interpolation interval 'EXPORT: x() array of interpolated x values ' y() array of corresponding interpolated y values ' maxpts number of data points in each the interpolated arrays '---------- DIM xi(1:100*MaxPts) AS SINGLE 'Splined x and y arrays DIM yi(1:100*MaxPts) AS SINGLE DIM b(1:MaxPts) AS SINGLE 'Temporary arrays DIM c(1:MaxPts) AS SINGLE DIM d(1:MaxPts) AS SINGLE DIM j AS INTEGER 'Points in input array DIM n AS INTEGER 'Points in the splined array DIM jtemp AS INTEGER DIM i AS INTEGER DIM ib AS INTEGER DIM ii AS INTEGER DIM jj AS INTEGER DIM kk AS INTEGER DIM k AS INTEGER DIM t AS SINGLE DIM u AS SINGLE DIM dx AS SINGLE '========== ' Start the SPLINE subroutine j = maxpts jtemp = j - 1 IF j < 2 THEN GOTO 299 END IF IF j < 3 THEN GOTO 250 END IF '---------- ' Set up tri-diagonal system ' b=diagonal, d=offdiagonal, c=RHS d(1) = x(2) - x(1) c(2) = (y(2) - y(1)) / d(1) FOR i = 2 TO jtemp d(i) = x(i + 1) - x(i) b(i) = 2.0 * (d(i - 1) + d(i)) c(i + 1) = (y(i + 1) - y(i)) / d(i) c(i) = c(i + 1) - c(i) NEXT i '---------- ' End conditions, 3rd der. at x(1) and x(n) ' obtained from divided differences b(1) = -d(1) b(j) = -d(j - 1) c(1) = 0.0 c(j) = 0.0 IF j = 3 THEN GOTO 215 END IF c(1) = c(3) / (x(4) - x(2)) - c(2) / (x(3) - x(1)) c(j) = c(j - 1) / (x(j) - x(j - 2)) - c(j - 2) / (x(j - 1) - x(j - 3)) c(1) = (c(1) * (d(1) ^ 2)) / (x(4) - x(1)) c(j) = -c(j) * (d(j - 1) ^ 2) / (x(j) - x(j - 3)) '---------- ' Forward elimination 215 FOR i = 2 TO j t = d(i - 1) / b(i - 1) b(i) = b(i) - t * d(i - 1) c(i) = c(i) - t * c(i - 1) NEXT i '---------- ' Back substitution c(j) = c(j) / b(j) FOR ib = 1 TO jtemp i = j - ib c(i) = (c(i) - d(i) * c(i + 1)) / b(i) NEXT ib '---------- ' Compute the polynomial coefficients b(j) = (y(j) - y(jtemp)) / d(jtemp) + d(jtemp) * (c(jtemp) + 2 * c(j)) FOR i = 1 TO jtemp b(i) = (y(i + 1) - y(i)) / d(i) - d(i) * (c(i + 1) + 2 * c(i)) d(i) = (c(i + 1) - c(i)) / d(i) c(i) = 3.0 * c(i) NEXT i c(j) = 3.0 * c(j) d(j) = d(j - 1) GOTO 299 250 b(1) = (y(2) - y(1)) / (x(2) - x(1)) c(1) = 0.0 d(1) = 0.0 b(2) = b(1) c(2) = 0.0 d(2) = 0.0 299 xi(1) = x(1) yi(1) = y(1) n = 2 30 xi(n) = xi(n - 1) + xint IF xi(n) > x(j) THEN GOTO 40 END IF u = xi(n) '========== ' Now do the SEVAL function ii = 1 IF ii >= j THEN ii = 1 END IF IF u < x(ii) THEN GOTO 310 END IF IF u <= x(ii + 1) THEN GOTO 330 END IF '---------- ' Binary search 310 ii = 1 jj = j + 1 320 kk = INT( (ii + jj) / 2 ) IF u < x(kk) THEN jj = kk END IF IF u >= x(kk) THEN ii = kk END IF '---------- ' NOTE: I don't know what purpose this next block serves, ' and it sometimes causes infinite loop errors. Maybe in very noisy data??? IF u > x(ii + 1) THEN GOTO 320 END IF '---------- ' Evaluate spline 330 dx = u - x(ii) yi(n) = y(ii) + dx * (b(ii) + dx * (c(ii) + dx * d(ii))) '========== n = n + 1 GOTO 30 '---------- ' Wind n back by one and output the SPLINTERED array 40 n = n - 1 maxpts = n 'Overwrites the inputs REDIM x(1 TO maxpts) REDIM y(1 TO maxpts) FOR k = 1 TO maxpts x(k) = xi(k) y(k) = yi(k) NEXT k 99 END SUB '====================================================================== SUB OriginalXY () EXPORT LOCAL i AS INTEGER LOCAL xTemp AS STRING LOCAL yTemp AS STRING gMinPts = %AbsMinPts : gMaxPts = %AbsMaxPts DIM gX(gMinPts : gMaxPts) DIM gY(gMinPts : gMaxPts) '---------- gMaxPts = gMinPts i = gMinPts WHILE i <= DATACOUNT xTemp = READ$(i) yTemp = READ$(i+1) gX(gMaxPts) = VAL(TRIM$(xTemp)) gY(gMaxPts) = VAL(TRIM$(yTemp)) IF gMaxPts = gMinPts THEN 'Set the max and min gXMin = gX(gMaxPts) : gXMax = gX(gMaxPts) 'Set the data Min and Max gYMin = gY(gMaxPts) : gYMax = gY(gMaxPts) ELSE SELECT CASE gY(gMaxPts) 'Find the max and min of the y array CASE gYMax gYMax = gY(gMaxPts) END SELECT END IF i = i + %MaxCols ' DATA counter INCR gMaxPts ' X,Y array index counter WEND DECR gMaxPts REDIM PRESERVE gX(gMinPts : gMaxPts) REDIM PRESERVE gY(gMinPts : gMaxPts) '---------- ARRAY SORT gX(), TAGARRAY gY() 'In ascending order gXMin = gX(gMinPts) gXMax = gX(gMaxPts) DATA 0, 1299 DATA 5, 1202 DATA 10, 1066 DATA 15, 896 DATA 20, 701 DATA 25, 488 DATA 30, 266 DATA 35, 45 DATA 40, -166 DATA 45, -359 DATA 50, -525 DATA 60, -756 DATA 70, -828 DATA 80, -739 DATA 90, -515 DATA 100, -206 DATA 110, 125 DATA 120, 416 DATA 130, 613 DATA 140, 683 DATA 150, 614 DATA 165, 293 DATA 180, -157 DATA 195, -542 DATA 210, -692 DATA 225, -522 DATA 240, -75 DATA 255, 492 DATA 275, 1072 DATA 295, 1098 DATA 315, 541 DATA 335, -240 DATA 355, -744 DATA 375, -674 DATA 400, 34 DATA 425, 680 DATA 450, 630 DATA 500, -517 DATA 550, 367 END SUB '====================================================================== SUB NewXY () EXPORT LOCAL i AS INTEGER i = gMinPts gXMin = gX(i) : gXMax = gX(i) 'Set the data Min and Max gYMin = gY(i) : gYMax = gY(i) INCR i WHILE i <= gMaxPts SELECT CASE gY(i) 'Find the maximum and minimum of the y array CASE gYMax gYMax = gY(i) END SELECT INCR i WEND '---------- ARRAY SORT gX(), TAGARRAY gY() 'In ascending order gXMin = gX(1) gXMax = gX(gMaxPts) END SUB '====================================================================== SUB PLOT(BYVAL hWnd AS LONG) LOCAL LpPaint AS PaintStruct LOCAL hDC AS LONG LOCAL hMemDC AS LONG LOCAL hBmp AS LONG LOCAL hBrush AS LONG LOCAL hPen AS LONG LOCAL hWndShow AS SINGLE LOCAL WinX AS LONG LOCAL WinY AS LONG LOCAL i AS LONG ' Since it is used in passing to API LOCAL Colr AS LONG LOCAL j AS LONG LOCAL r AS RECT ' For Window dimensions LOCAL X0 AS LONG ' Temporary origin to draw FROM LOCAL Y0 AS LONG LOCAL X1 AS LONG ' Point to draw TO LOCAL Y1 AS LONG LOCAL XMarg AS LONG LOCAL YMarg AS LONG LOCAL GMarg AS LONG LOCAL XWide AS LONG LOCAL YWide AS LONG LOCAL XUnit AS SINGLE LOCAL YUnit AS SINGLE '---------- 'Prepare window and background GetClientRect hWnd, r WinX = r.nRight - r.nLeft 'Get windows dimensions WinY = r.nBottom - r.nTop hDC = BeginPaint(hWnd, LpPaint) hMemDC = CreateCompatibleDC(hDC) hBmp = CreateCompatibleBitmap(hDC, WinX, WinY) SelectObject hMemDC, hBmp hBrush = CreateSolidBrush(%White) FillRect hMemDC, r, hBrush DeleteObject hBrush '---------- 'Initializations XMarg = WinX/%XMargin 'Margins are 10% of the Window YMarg = WinY/%YMargin GMarg = WinY*(%GMargin/100) XWide = WinX - 2*XMarg YWide = WinY - 2*YMarg XUnit = (gXMax-gXMin) / CSNG(XWide-(2*GMarg)) YUnit = (gYMax-gYMin) / CSNG(YWide-(2*GMarg)) '---------- 'For the border hPen = CreatePen(%PS_SOLID, %BorderThick, %Black) SelectObject hMemDC, hPen MoveToEx hMemDC, XMarg, YMarg, BYVAL %NULL LineTo hMemDC, XMarg+XWide, YMarg MoveToEx hMemDC, XMarg+XWide, YMarg, BYVAL %NULL LineTo hMemDC, XMarg+XWide, YMarg+YWide MoveToEx hMemDC, XMarg+XWide, YMarg+YWide, BYVAL %NULL LineTo hMemDC, XMarg, YMarg+YWide MoveToEx hMemDC, XMarg, YMarg+YWide, BYVAL %NULL LineTo hMemDC, XMarg, YMarg '---------- 'For the graph Colr = %Black hPen = CreatePen(%PS_SOLID, %GraphThick, Colr) SelectObject hMemDC, hPen XMarg = XMarg+GMarg 'i.e. the inside margins between the border YMarg = YMarg-GMarg ' and the actual graph X0 = XMarg + CLNG ( (gX(gMinPts)-gXMin) / XUnit ) 'First point Y0 = Ymarg + (YWide - CLNG ( (gY(gMinPts)-gYMin) / YUnit ) ) Ellipse hMemDC, X0-%Radius, Y0-%Radius, X0+%Radius, Y0+%Radius FOR j = gMinPts TO gMaxPts 'And the points thereafter X1 = XMarg + CLNG ( (gX(j)-gXMin) / XUnit ) 'Point TO which to draw Y1 = YMarg + (YWide - CLNG ( (gY(j)-gYMin) / YUnit ) ) ' becomes next origin MoveToEx hMemDC, X0, Y0, BYVAL %NULL 'Move to start LineTo hMemDC, X1, Y1 ' and draw line from it to next Ellipse hMemDC, X1-%Radius, Y1-%Radius, X1+%Radius, Y1+%Radius X0 = X1 'Point FROM which to draw Y0 = Y1 'Origin point NEXT j '---------- 'Clean up BitBlt hDC, r.nLeft, r.nTop, WinX, WinY, hMemDC, 0, 0, %SRCCOPY EndPaint hWndShow, LpPaint DeleteObject hPen END SUB '====================================================================== FUNCTION WinMain (BYVAL hInstance AS LONG, _ BYVAL hPrevInstance AS LONG, _ lpCmdLine AS ASCIIZ PTR, _ BYVAL iCmdShow AS LONG) AS LONG LOCAL Msg AS tagMsg LOCAL wndclass AS WndClassEx LOCAL szAppName AS ASCIIZ * 80 LOCAL hWnd AS LONG szAppName = "GDIWINDOW" wndclass.cbSize = SIZEOF(WndClass) wndclass.style = %CS_HREDRAW OR %CS_VREDRAW OR %CS_DBLCLKS wndclass.lpfnWndProc = CODEPTR( Spline_WndProc ) wndclass.cbClsExtra = 0 wndclass.cbWndExtra = 0 wndclass.hInstance = hInstance wndclass.hIcon = LoadIcon( hInstance, BYVAL %IDI_APPLICATION ) wndclass.hCursor = LoadCursor( %NULL, BYVAL %IDC_ARROW ) wndclass.hbrBackground = GetStockObject( %BLACK_BRUSH ) wndclass.lpszMenuName = %NULL wndclass.lpszClassName = VARPTR( szAppName ) wndclass.hIconSm = LoadIcon( hInstance, BYVAL %IDI_APPLICATION ) RegisterClassEx wndclass ' Create a window using the registered class hWnd = CreateWindow("GDIWINDOW",_ ' window class name "Double-click to interpolate, Right-click for original data", _ ' window caption %WS_OVERLAPPEDWINDOW, _ ' window style 0, _ ' initial x position 0, _ ' initial y position GetSystemMetrics(%SM_CXSCREEN),_ ' initial x size GetSystemMetrics(%SM_CYSCREEN)/2,_ ' initial y size %NULL, _ ' parent window handle %NULL, _ ' window menu handle hInstance, _ ' program instance handle BYVAL %NULL) ' creation parameters ShowWindow hWnd, iCmdShow UpdateWindow hWnd WHILE GetMessage(Msg, %NULL, 0, 0) TranslateMessage Msg DispatchMessage Msg WEND FUNCTION = msg.wParam 'FUNCTION = 0 means failure - stop END FUNCTION '====================================================================== FUNCTION Spline_WndProc (BYVAL hWnd AS LONG, BYVAL wMsg AS LONG, _ BYVAL wParam AS LONG, BYVAL lParam AS LONG) EXPORT AS LONG LOCAL hDC AS LONG LOCAL LpPaint AS PaintStruct LOCAL tRect AS Rect STATIC First AS LONG SELECT CASE wMsg CASE %WM_NCCREATE CALL OriginalXY() CASE %WM_CREATE CASE %WM_SIZE InvalidateRect hWnd, BYVAL %NULL, %TRUE CASE %WM_LBUTTONDOWN InvalidateRect hWnd, BYVAL %NULL, %TRUE CASE %WM_RBUTTONDOWN 'Show original array CALL OriginalXY() InvalidateRect hWnd, BYVAL %NULL, %TRUE CASE %WM_LBUTTONDBLCLK 'Interpolate the XY arrays gxInt = 2 CALL Splint(gX(), gY(), gMaxPts, gxInt) CALL NewXY() InvalidateRect hWnd, BYVAL %NULL, %TRUE CASE %WM_ERASEBKGND IF First = 0 THEN First = 1 ELSE FUNCTION = 1 EXIT FUNCTION END IF CASE %WM_PAINT Plot hWnd CASE %WM_SYSCOMMAND IF wParam = %SC_CLOSE THEN DestroyWindow hWnd EXIT FUNCTION END IF CASE %WM_DESTROY PostQuitMessage 0 FUNCTION = 0 EXIT FUNCTION END SELECT FUNCTION = DefWindowProc(hWnd, wMsg, wParam, lParam) 'VITAL - must give function a meaning END FUNCTION '======================================================================