CREATE OR REPLACE FUNCTION ST_Extend(eje_ geometry, bound_ geometry)
RETURNS geometry 
AS $$
-- Use Case: I need to "extend" a linestring in another one. The new linestring must be
-- obtained such as their new extremes are the intersection of it with a
-- polygon that contains it. The original linestring is composed of only a
-- linear segment. 
-- version: alfa , by Julio A. Galindo, April 17, 2007: juliogalindoq at gmail.com
DECLARE
     b_ geometry = st_boundary(bound_);
     dist float;
     max_dist float = 0;
     n_points int;
     pto_1 geometry;
     pto_2 geometry;
     first_pto geometry;
     last_pto geometry;
     u_1 float;
     u_2 float;
     norm float;
     result text = 'LINESTRING(';
BEGIN
     IF ( ST_GeometryType(eje_)   NOT LIKE 'ST_LineString' OR 
          ST_GeometryType(bound_) NOT LIKE 'ST_Polygon' ) THEN 
        RETURN NULL; 
     END IF;
     -- First Search how far is the boundary: (worst case)
     pto_1 := ST_StartPoint(eje_);
     pto_2 := ST_EndPoint(eje_);
     FOR i IN 1..ST_NumPoints(b_)-1 LOOP
        dist := ST_distance(ST_PointN(b_,i),pto_1);
        IF dist > max_dist THEN 
           max_dist := dist; 
        END IF;
        dist := ST_Distance(ST_PointN(b_,i),pto_2);
        IF dist > max_dist THEN 
           max_dist := dist; 
        END IF;
     END LOOP;
     -- Now extent the linestring:
     pto_2 := ST_PointN(eje_,2);
     u_1 := ST_X(pto_2)-ST_X(pto_1);
     u_2 := ST_Y(pto_2)-ST_Y(pto_1);
     norm := sqrt(u_1^2 + u_2^2);
     first_pto := ST_MakePoint(ST_X(pto_1)-u_1/norm*dist,ST_Y(pto_1)-u_2/norm*dist);
     n_points := ST_nPoints(eje_);
     IF n_points > 2 THEN
        pto_1 := ST_PointN(eje_,n_points-1);
        pto_2 := ST_PointN(eje_,n_points);
        u_1   := ST_X(pto_2)-ST_X(pto_1);
        u_2   := ST_Y(pto_2)-ST_Y(pto_1);
        norm  := sqrt(u_1^2 + u_2^2);
     END IF;
     last_pto := ST_MakePoint(ST_X(pto_2)+u_1/norm*dist,ST_Y(pto_2)+u_2/norm*dist);
     result   := result || ST_X(first_pto) || ' ' || ST_Y(first_pto) || ',';
     FOR i IN 1..ST_NumPoints(eje_) LOOP
         result := result || ST_X(ST_PointN(eje_,i)) || ' ' || ST_Y(ST_PointN(eje_,i)) || ',';
     END LOOP;
     result := result || ST_X(last_pto) || ' ' || ST_Y(last_pto) || ')';
     -- Find the final Linestring:
     b_ := ST_intersection(ST_GeomFromText(result,ST_SRID(eje_)),bound_);
     RETURN b_;
END $$
LANGUAGE plpgsql
STABLE
RETURNS NULL ON NULL INPUT;

-- Test 1:
select ST_asText(ST_Extend(ST_geomFromText('LINESTRING(5 5, 12 12)',-1),ST_GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));

select ST_asText(ST_extend(ST_geomFromText('LINESTRING(5 5, 12 12, 13 12)',-1),ST_GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));


