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:
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:
‘{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.:
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.
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;
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
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?
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)