Algorithm GRI  Robust Geometric Resampling Algorithm (RGR)  Smoothing-Resampling Algorithm (SRALG)  SRALG-Based Interpolation Algorithm  Free Dot Gain Compensation Calculator  Free Ink Balance Online Tool  Free Utilities About the Author

SRALG-Based Interpolation Algorithm

If you wish to use the materials on this page, you are obliged to preserve the link to it.

The algorithm offered below solves the problem of finding a set of intermediate values according to a given set of original values. It is effective when there is a need in the set of values mentioned above without a requirement for a defined interpolation function. That is indeed the case in a large number of computational problems.

In terms of signal processing, the proposed algorithm is, in fact, an algorithm of resampling. In many cases it offers an acceptable result with only a small amount of iterations, thus cutting expenses on computation time. It is based on the SRALG I proposed, designed for smoothing. The core idea is utilizing SRALG's innate feature (preservation of points lying in the middle of a line segment defined by its neighbors) to preserve the original set of values, yet allowing it to run any amount of iterations, till the desired smoothness, precision or any other desired criteria is achieved.

Thus, the essence of this algorithm is in fact the pre-iterative construction made, which defines the shape of the interpolating function. Certainly enough, using SRALG's innate feature it is possible to offer different methods of pre-iteration construction to produce a different kinds of curves at the end of the iteration cycle. Still, the construction offered below proves effective enough in many cases.

Pre-Iteration construction principles:

  • Each of the original points is surrounded by a pair of "neighbors", so that it lies in the middle of line segment defined by them.
  • The "neighbors" are chosen so that they lie on bisectors of angles between lines defined by consequent original points. It ensures that resultant curve will appear natural.
  • The distance between the "neighbors" and the original point are defined by an initial parameter in order to have control over the deviation magnitude from the original set of points.
pre-iteration illustration

D.Grilikhes, PhD
March 24, 2013



MATLAB realization of the pre-iteration construction

function Y = preiteration( X,t )

%PREITERATION prepares the pre-iterational set of points

% for SRALG-based interpolation algorithm.

% preiteration(X,t) uses original data points as array 

% whose each row is a pooint (x,y) and t - the parameter

% ( it should be between 0 and 1 )

% defining how close to original nodes would be the additional points.


% D.Grilikhes 03-24-2013

% Copyright 2013 Dmitry Grilikhes


Y = zeros((size(X,1)-2)*3+2,2); % memory allocating


% start first segment processing 

Y (1,1:2) = X(1,1:2);

mid = middle(X(1,1:2),X(2,1:2));

slope0=-1/coef(X(1,1:2),X(2,1:2));

if ~isinf(slope0) && ~all(isinf(inters1(X(2,1:2),X(3,1:2),mid,slope0)))

    vertex = inters1(X(2,1:2),X(3,1:2),mid,slope0);

    if vertex(1,1)<= mid(1,1) && vertex(1,1)>=X(2,1)

        vertex = inters2(X(2,1:2),X(3,1:2),mid);

    end

else

    vertex = inters2(X(2,1:2),X(3,1:2),mid);

end

Y=addtores(shift(incenter(X(2,1:2),mid,vertex),X(2,1:2),t),Y,2);

 % end first segment processing


% start main part processing

m=3;

for n = 2:(size(X,1)-2)

[Y,m]=addtores(X(n,1:2),Y,m);

    vertex = inters(X(n,1:2),X(n-1,1:2),X(n+1,1:2),X(n+2,1:2));

    if X(n,1)<vertex(1,1) && vertex(1,1)<X(n+1,1) % the simplest case

        tmp1=shift(incenter(X(n,1:2),vertex,X(n+1,1:2)),X(n,1:2),t);

        if d(tmp1,X(n,1:2))~=0

            if d(tmp1,X(n,1:2))>d(Y(m-2,1:2),X(n,1:2))

                tmp1=osp(Y(m-2,1:2),X(n,1:2));

            else

                Y(m-2,1:2)=osp(tmp1,X(n,1:2));

            end

        end

        [Y,m]=addtores(tmp1,Y,m);

        [Y,m]=addtores(shift(incenter(X(n,1:2),vertex,X(n+1,1:2)),X(n+1,1:2),t),Y,m);

    else

        if  all(isinf(vertex)) % neighbor segments are parallel

            slope1=-1/coef(X(n,1:2),X(n+1,1:2));

        else

            slope1=-1/coef(vertex,mid);

        end

        mid = middle(X(n,1:2),X(n+1,1:2));

        if ~isinf(slope1)

        vertex1 = inters1(X(n,1:2),X(n-1,1:2),mid,slope1);

           if vertex1(1,1) > X(n,1) && vertex1(1,1)<mid(1,1)

                 vertex2 = inters1(X(n+1,1:2),X(n+2,1:2),mid,slope1);

           else

                vertex1 = inters2(X(n,1:2),X(n-1,1:2),mid);

                vertex2 = inters2(X(n+1,1:2),X(n+2,1:2),mid);

           end

        else

           vertex1 = inters2(X(n,1:2),X(n-1,1:2),mid);

           vertex2 = inters2(X(n+1,1:2),X(n+2,1:2),mid); 

        end

        tmp2=shift(incenter(X(n,1:2),mid,vertex1),X(n,1:2),t);

        if d(tmp2,X(n,1:2))~=0

            if d(tmp2,X(n,1:2))>d(Y(m-2,1:2),X(n,1:2))

                tmp2=osp(Y(m-2,1:2),X(n,1:2));

            else

                Y(m-2,1:2)=osp(tmp2,X(n,1:2));

            end

        end

        [Y,m]=addtores(tmp2,Y,m);

        [Y,m]=addtores(shift(incenter(X(n+1,1:2),mid,vertex2),X(n+1,1:2),t),Y,m);

    end

end

% end main part processing

% start last segment processing

Y (end-2,1:2) = X(end-1,1:2);

mid = middle(X(end,1:2),X(end-1,1:2));

slope2=-1/coef(X(end-1,1:2),X(end,1:2));

if ~isinf(slope2) && ~all(isinf(inters1(X(end-2,1:2),X(end-1,1:2),mid,slope2)))

    vertex = inters1(X(end-2,1:2),X(end-1,1:2),mid,slope2);

    if vertex(1,1) >= mid(1,1) && vertex(1,1)<=X(end-2,1)

         vertex = inters2(X(end-2,1:2),X(end-1,1:2),mid);

    end

else

    vertex = inters2(X(end-2,1:2),X(end-1,1:2),mid);

end

        tmp3=shift(incenter(X(end-1,1:2),mid,vertex),X(end-1,1:2),t);

        if d(tmp3,X(end-1,1:2))~=0

            if d(tmp3,X(end-1,1:2))>d(Y(end-3,1:2) ,X(end-1,1:2))

                tmp3=osp(Y(end-3,1:2),X(end-1,1:2));

            else

                Y(end-3,1:2)=osp(tmp3,X(end-1,1:2));

            end

        end

 Y=addtores(tmp3,Y,size(Y,1)-1);

 Y (end,1:2) = X(end,1:2);

% end last segment processing


end


function k= coef(A,B)

%coef(a,b) calculates the slope k a line defined by two points A and B


k = (B(1,2)-A(1,2))/(B(1,1)-A(1,1));

end


function di=d(p1,p2)

%d(p1,p2) calculates the distance between two points p1 and p2.


di=sqrt((p1(1,1)-p2(1,1))^2+(p1(1,2)-p2(1,2))^2);

end


function center = incenter( A,B,C )

%INCENTER finds incenter of a triangle defined by A, B and C points.


x = (d(B,C)*A(1,1)+d(A,C)*B(1,1)+d(B,A)*C(1,1))/(d(B,C)+d(A,C)+d(B,A));

y = (d(B,C)*A(1,2)+d(A,C)*B(1,2)+d(B,A)*C(1,2))/(d(B,C)+d(A,C)+d(B,A));

center = [x y ];

end


function is= inters(f1,f2,g1,g2)

%inters(f1,f2,g1,g2) finds the intersection point of the lines f and g, each

%defined respectively by the points f1, f2 and g1, g2.


x = (g1(1,2)-f1(1,2)-coef(g1,g2)*g1(1,1)+coef(f1,f2)*f1(1,1))/(coef(f1,f2)-coef(g1,g2));

y =  coef(f1,f2)*(x-f1(1,1))+f1(1,2);

is = [ x y ];

end


function is= inters1(f1,f2,g1,k)

%inters(f1,f2,g1,k) finds the intersection point of the lines f and g, when

% f is defined by points f1, f2 and g is defined by a point g1 and a slope

% k.


x = (g1(1,2)-f1(1,2)-k*g1(1,1)+coef(f1,f2)*f1(1,1))/(coef(f1,f2)-k);

y =  coef(f1,f2)*(x-f1(1,1))+f1(1,2);

is = [ x y ];

end


function Y = inters2(f1,f2,g)

%inters2(f1,f2,g) finds the intersection point of a line defined by points

% f1, f2 and a vertical line through the point g.

    Y= [g(1,1) (g(1,1)-f2(1,1))*(f1(1,2)-f2(1,2))/(f1(1,1)-f2(1,1))+f2(1,2)]; 

end


function M= middle(P,N)

%middle(P,N) finds the middle point of the line segment [P,N]


M = (P+N)/2;

end


function C= shift( A,B,k )

%shift( A,B,k ) computes a point between A and B so that |C-A|/|B-A|=k

%  Where A and B are points with coordinates (x,y) and k - a given

%  coefficient. k is assumed to be between 0 and 1.


C = A*(1-k)+B*k;

end


function [Y,m]=addtores(p,Y,m)

%addtores(X,Y,m) adds the point p with coordinates (x,y) to the the matrix Y

% as row m. Returns Y and the number of the next row.

% If the point is the same as previous it does not add the point to the matrix and 

% removes the unused row.

    if ~isequal(p,Y(m-1,1:2))

        Y(m,1:2) = p;

        m=m+1;

        else

            Y(m,:)=[];

    end

end


function p3= osp(p1,p2)

%osp(P,N) finds the other side point

% so that p2 is the middle of the line segment

% defined by the result and and p1


p3 = 2*p2-p1;

end

MATLAB realization of SRALG-based interpolation algorithm.

function Y = sralgin( X,t,v,k )

%SRALGIN implements SRALG-based interpolation.

% sralgin( X,t,v,k ) does the folowing:

%  1. Create a pre-iterational set of points from a matrix X of the

% original data, where each row is a point (x,y), and t - a proximity

% parameter defining how close the additional points would be to the

% original data. 

% It is assumed that t is 0<t<1 where 0 is as far as the algorythm permits

% and 1 doesn't change the original curve.

%  2. Apply SRALG k times with smothness parameter v (0<v<1).

% The result is the resampling interpolational matrix Y, where each row is

% a point (x,y).


% D.Grilikhes 03-24-2013

% Copyright 2013 Dmitry Grilikhes 


Y=preiteration(X,t);

for i = 1:k

    Y = iterSRALG(Y,v);

end


end

* MATLAB realization of iterSRALG may be found there.

Illustrations

Fig.1 Comparison of SRALG-based interpolation and the traditional methods based on Hermit polynomials and third-degree splines

SRALG illustration

Fig.2 Pre-iterational state and comparison of SRALG-based interpolation and the traditional methods based on Hermit polynomials and third-degree splines

SRALG illustration

© Dmitry Grilikhes, 2013

If you wish to use the materials on this page, you are obliged to preserve the link to it.

www.000webhost.com