PLSQL for Curve Smoothing

In an engineering project, I have been assigned to a task constructing smooth curve for given sets of sample points. This is accomplished by implementing Cubic Splines and Thomas algorithm to solve linear algebraic equation. So, I write PLSQL stored procedure as follows:

CREATE OR REPLACE FUNCTION spline(
IN x_arr numeric[], IN y_arr numeric[], IN resolution numeric,
OUT xs numeric, OUT ys numeric)
RETURNS SETOF record AS
$BODY$
DECLARE
  x numeric[];
  y numeric[];
  e numeric[] DEFAULT ‘{}’::numeric[];
  f numeric[] DEFAULT ‘{}’::numeric[];
  g numeric[] DEFAULT ‘{}’::numeric[];
  r numeric[] DEFAULT ‘{}’::numeric[];
  f2 numeric[] DEFAULT ‘{}’::numeric[];
  count integer;
  i integer;
  j integer;
  idx integer;
  tmp numeric;
BEGIN
  x:=x_arr;
  y:=y_arr;
  count:=ARRAY_UPPER(x,1);
  IF count<3 THEN
    RAISE EXCEPTION ‘sample size must be >=3’;
  END IF;
  IF ARRAY_UPPER(y,1)<>count THEN
    RAISE EXCEPTION ‘x and y array size must be identical’;
  END IF;

——–sort x, ascendingly:
  FOR i IN 1..count-1 LOOP
    idx:=i+1;
    FOR j IN i+2..count LOOP
      IF x[j]<x[idx] THEN
        idx:=j;
      END IF;
    END LOOP;
    IF x[idx]=x[i] THEN
      RAISE EXCEPTION ‘x values must be unique’;
    ELSIF x[idx]<x[i] THEN –swap
      tmp=x[i];
      x[i]=x[idx];
      x[idx]=tmp;
      tmp=y[i];
      y[i]=y[idx];
      y[idx]=tmp;
    END IF;
  END LOOP;

——–prepare variables for Thomas algorithm
  FOR i IN 1..count LOOP
    IF i=1 OR i=count THEN
      e:= e || 0.0;
      f:= f || 0.0;
      g:= g || 0.0;
      r:= r || 0.0;
      f2:=f2 || 0.0;
      CONTINUE;
    END IF;
    e:= e || x[i] – x[i-1];
    f:= f || 2*(x[i+1]-x[i-1]);
    g:= g || x[i+1] – x[i];
    r:= r || ( (6/(x[i+1]-x[i])) * (y[i+1]-y[i]) )
      + ( (6/(x[i]-x[i-1])) * (y[i-1]-y[i]) );
    f2:=f2 || 0.0;
  END LOOP;

——–decomposition
  FOR i IN 2..count LOOP
    IF f[i-1]<>0 THEN
      e[i]:=e[i]/f[i-1];
    END IF;
    f[i]:=f[i]-(e[i]*g[i-1]);
  END LOOP;

——–forward substitution
  FOR i IN 2..count LOOP
    r[i]:=r[i]-(e[i]*r[i-1]);
  END LOOP;

——–backward substitution
  FOR i IN REVERSE count-1..1 LOOP
    IF f[i]<>0 THEN
      f2[i]:=(r[i] – (g[i]*f2[i+1])) / f[i];
    END IF;
  END LOOP;

——–do interpolation
  FOR i IN 2..count LOOP
    xs:=x[i-1]-resolution;
    LOOP
      xs:=xs+resolution;
      EXIT WHEN i=count AND xs>x[i];
      EXIT WHEN i=x[i];
      ys:= ( (f2[i-1] * (x[i]-xs)^3)
        + (f2[i] * (xs-x[i-1])^3)
        ) / (6*(x[i]-x[i-1]));
      ys:= ys + (
        (y[i-1]/(x[i]-x[i-1])) –
        (f2[i-1]*(x[i]-x[i-1])/6)
        ) * (x[i]-xs);
      ys:= ys + (
        (y[i]/(x[i]-x[i-1])) –
        (f2[i]*(x[i]-x[i-1])/6)
        ) * (xs-x[i-1]);
      RETURN NEXT;
      EXIT WHEN xs>=x[i];
    END LOOP;
  END LOOP;

  RETURN;
END;
$BODY$
LANGUAGE ‘plpgsql’ IMMUTABLE STRICT;

The stored procedure takes 3 input arguments and 2 outs. x_arr and y_arr are sample points’ x and y. With resolution, you can determine the resolution of smoothed curve’s x.

Output arguments: xs and ys contain the set of x and y of smoothed curve, respectively.

Let us use the function. At first, we have sample points:

#	x	y
1	3.3	1.2
2	4.2	-1.1
3	6.8	1.5
4	8.2	1.2
5	10.2	2.5

Call the function, supply x and y array as displayed on above table, set resolution to 0.1:

SELECT * FROM spline(
  ‘{3.3,4.2,6.8,8.2,10.2}’::numeric[],
   ‘{1.2,-1.1,1.5,1.2,2.5}’::numeric[],
  0.1);

We will have set of records containing xs and ys.:

xs  ys
3.3  1.199999999999999970
3.4  0.8860175045157231811
3.5  0.5764170295261004851
3.6  0.2755805955257859750
3.7  -0.0121097769905662565
3.8  -0.28227206752830211620
… and so on

Let us put them in graph:

The blue-y1 and red-y2 are sample and interpolation points as returned by the spline function, respectively.