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:
curve

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

Advertisements
4 comments
  1. Günter Glosse said:

    Hi, please tell me, if this ok.
    The Script run over 1 hour. Thanks Günter

    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′;
    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;
    END LOOP;
    END LOOP;

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

  2. Gunter,

    I have noticed that you have mistyped at several places:
    1. Nearby lines “IF count=3” at eraly rows below BEGIN.
    2. At section –decomposition “IF f[i-1]0”
    3. At section –backward substitution “IF f[i]0”

    Other things: let me know your sample points and resolution so I can replay your case in my PC.

    Thanks

  3. Krish said:

    Hello Abdul,

    Thank you for this function. I am wondering if you would know how to convert this to PL/SQL (Oracle)? Or if you happen to have a PL/SQL version of the same?

  4. JK.Yang said:

    Hi. I convert this to PL/SQL (Oracle).
    I hope that this will help.

    — oracle array type
    create or replace type type_array as table of float;
    — spline conversion to oracle statement
    create or replace procedure spline(
    x_arr in type_array := type_array(3.3,4.2,6.8,8.2,10.2),
    y_arr in type_array := type_array(1.2,-1.1,1.5,1.2,2.5),
    resolution in float := 0.1)

    as
    x type_array := type_array();
    y type_array := type_array();
    e type_array := type_array();
    f type_array := type_array();
    g type_array := type_array();
    r type_array := type_array();
    f2 type_array := type_array();
    v_count float;
    i float;
    j float;
    idx float;
    tmp float;
    xs float;
    ys float;
    BEGIN
    x:=x_arr;
    y:=y_arr;
    v_count:=x.count;

    IF v_count=3′);
    END IF;
    IF y.count x.count THEN
    RAISE_APPLICATION_ERROR(-20020, ‘x and y array size must be identical’);
    END IF;

    — sort x, ascendingly:
    FOR i IN 1..v_count-1 LOOP
    idx:=i+1;
    FOR j IN i+2..v_count LOOP
    IF x(j)<x(idx) THEN
    idx:=j;
    END IF;
    END LOOP;
    IF x(idx)=x(i) THEN
    RAISE_APPLICATION_ERROR(-20030,'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..v_count LOOP
    e.extend;
    f.extend;
    g.extend;
    r.extend;
    f2.extend;
    IF i=1 OR i=v_count THEN
    e(i) := 0.0;
    f(i) := 0.0;
    g(i) := 0.0;
    r(i) := 0.0;
    f2(i) := 0.0;
    CONTINUE;
    END IF;
    e(i) := x(i) – x(i-1);
    f(i) := 2*(x(i+1)-x(i-1));
    g(i) := x(i+1) – x(i);
    r(i) := ( (6/(x(i+1)-x(i))) * (y(i+1)-y(i)) )
    + ( (6/(x(i)-x(i-1))) * (y(i-1)-y(i)) );
    f2(i) := 0.0;
    END LOOP;
    — decomposition
    FOR i IN 2..v_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..v_count LOOP
    r(i):=r(i)-(e(i)*r(i-1));
    END LOOP;

    — backward substitution
    FOR i IN REVERSE 1..v_count-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..v_count LOOP
    xs:=x(i-1)-resolution;
    LOOP
    xs:=xs+resolution;
    EXIT WHEN i=v_count AND xs>x(i);
    EXIT WHEN i=x(i);
    ys:= ( (f2(i-1) * power (x(i)-xs,3 ) )
    + (f2(i) * power(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;
    dbms_output.put_line (‘(xs,ys) = (‘ || xs ||’,’ || ys || ‘)’);
    EXIT WHEN xs>=x(i);
    END LOOP;
    END LOOP;

    END;

    — execution
    — default data is used. Check Definition of procedure
    exec spline();

    — result
    (xs,ys) = (3.3,1.2)
    (xs,ys) = (3.4,.886017504515723205775274835291278152328)
    (xs,ys) = (3.5,.5764170295261005044507373912690437765779)
    (xs,ys) = (3.6,.2755805955257859889265753886197843446386)
    (xs,ys) = (3.7,-.0121097769905662478970234519700126715761)
    (xs,ys) = (3.8,-.2822720675283021131198714098138598001583)
    (xs,ys) = (3.9,-.5305242555927675138417807642252695691978)
    (xs,ys) = (4,-.7524843206893083571625637945177545067892)
    (xs,ys) = (4.1,-.9437702423232705501820327800048271410326)
    (xs,ys) = (4.2,-1.1)
    (xs,ys) = (4.2,-1.1)
    (xs,ys) = (4.3,-1.21792938836992002242755873887690823838)
    (xs,ys) = (4.4,-1.29886546266376356812092630124967751715)
    (xs,ys) = (4.5,-1.34525309325734099644760099679255612388)
    (xs,ys) = (4.6,-1.35953715052646266677508113517979234614)
    (xs,ys) = (4.7,-1.3441625048469389384708650260856344715)
    (xs,ys) = (4.8,-1.30157402659458017090245097918433078751)
    (xs,ys) = (4.9,-1.23421658614519672343733730415012958175)
    (xs,ys) = (5,-1.14453505387459895544302231065727914179)
    (xs,ys) = (5.1,-1.0349743001585972262870043083800277552)
    (xs,ys) =
    (5.2,-.90797919537300189533678160699262370953)
    (xs,ys) = (5.3,-.76599460989362332195985251616931529236)
    (xs,ys) = (5.4,-.61146541409627186552371534558435079126)
    (xs,ys) = (5.5,-.44683647835675788539586840491197849378)
    (xs,ys) = (5.6,-.2745526730508917409438100038264466875)
    (xs,ys) = (5.7,-.09705886855448379153503845200200365998)
    (xs,ys) = (5.8,.08320006475665560346294794088710230119)
    (xs,ys) = (5.9,.26377925650671608468265086516662290848)
    (xs,ys) = (6,.4422338363198872927565720111623098743)
    (xs,ys) = (6.1,.61611893382035886831721306919991491107)
    (xs,ys) = (6.2,.78298967863232045199707572960518973126)
    (xs,ys) = (6.3,.94040120037996168442866168270388604728)
    (xs,ys) = (6.4,1.08590862868747220624447261882175557155)
    (xs,ys) = (6.5,1.21706709317904165807701022828455001654)
    (xs,ys) = (6.6,1.33143172347885968055877620141802109465)
    (xs,ys) = (6.7,1.42655764921111591432227222854792051834)
    (xs,ys) = (6.8,1.50000000000000000000000000000000000001)
    (xs,ys) = (6.8,1.5)
    (xs,ys) = (6.9,1.55015855192322326175045119570038334475)
    (xs,ys) = (7,1.5788116668725837578360774539766827282)
    (xs,ys) =
    (7.1,1.58858235319340123004532040275688241863)
    (xs,ys) = (7.2,1.58209361923099542016662166996896668431)
    (xs,ys) = (7.3,1.56196847333068606998842288354091979351)
    (xs,ys) = (7.4,1.53082992383779292129916567140072601451)
    (xs,ys) = (7.5,1.49130097909763571588729166147636961559)
    (xs,ys) = (7.6,1.446004647455534195541242481695834865)
    (xs,ys) = (7.7,1.39756393725680810204945975998710603103)
    (xs,ys) = (7.8,1.34860185684677717720038512427816738195)
    (xs,ys) = (7.9,1.30174141457076116278246020249700318604)
    (xs,ys) = (8,1.25960561877407980058412662257159771156)
    (xs,ys) = (8.1,1.22481747780205283239382601242993522679)
    (xs,ys) = (8.2,1.2)
    (xs,ys) = (8.2,1.2)
    (xs,ys) = (8.3,1.18723407839306030591314230667478132537)
    (xs,ys) = (8.4,1.18643214472564979553195502770728860068)
    (xs,ys) = (8.5,1.19696451542200377497719235181553624963)
    (xs,ys) = (8.6,1.21820150690635755036960846771753869589)
    (xs,ys) = (8.7,1.24951343560294642782995756413131036315)
    (xs,ys) = (8.8,1.29027061793600571347899382977486567511)
    (xs,ys) = (8.9,1.33984337032977071343747145336621905545)
    (xs,ys) =
    (9,1.39760200920847673382614462362338492785)
    (xs,ys) = (9.1,1.46291685099635908076576752926437771601)
    (xs,ys) = (9.2,1.53515821211765306037709435900721184361)
    (xs,ys) = (9.3,1.61369640899659397878087930156990173433)
    (xs,ys) = (9.4,1.69790175805741714209787654567046181187)
    (xs,ys) = (9.5,1.78714457572435785644884028002690649991)
    (xs,ys) = (9.6,1.88079517842165142795452469335725022214)
    (xs,ys) = (9.7,1.97822388257353316273568397437950740225)
    (xs,ys) = (9.8,2.07880100460423836691307231181169246393)
    (xs,ys) = (9.9,2.18189686093800234660744389437181983085)
    (xs,ys) = (10,2.28688176799906040793955291077790392671)
    (xs,ys) = (10.1,2.3931260422116478570301535497479591752)
    (xs,ys) = (10.2,2.5)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: