-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-- 
-- $Id: cleanGeometry.sql 2008-04-24 10:30Z Dr. Horst Duester $
--
-- cleanGeometry - remove self- and ring-selfintersections from 
--                 input Polygon geometries 
-- http://www.sogis.ch
-- Copyright 2008 SO!GIS Koordination, Kanton Solothurn, Switzerland
-- Version 1.0
-- contact: horst dot duester at bd dot so dot ch
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--  
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


CREATE OR REPLACE FUNCTION cleanGeometry(geometry)
  RETURNS geometry AS
$BODY$DECLARE
  inGeom ALIAS for $1;
  outGeom geometry;
  tmpLinestring geometry;
  tmpPolygon geometry;
  Pk geometry;
  Pi geometry;
  nGeometries integer;
  nGeom integer;
  intRingPi geometry;
  nIntRings integer;
  extRingPk geometry;
  isHole boolean;

Begin
  
  outGeom := NULL;
  
-- Clean Process for Polygon 
  IF (GeometryType(inGeom) = 'POLYGON' OR GeometryType(inGeom) = 'MULTIPOLYGON') THEN

-- Only process if geometry is not valid, 
-- otherwise put out without change
    if not isValid(inGeom) THEN
    
-- create nodes at all self-intersecting lines by union the polygon boundaries
-- with the startingpoint of the boundary.  
      tmpLinestring := st_union(st_multi(st_boundary(inGeom)),st_pointn(boundary(inGeom),1));
      
-- Recreate the polygons with polygonize
-- => all POLYGON, and all holes are POLYGON now!      
      tmpPolygon = polygonize(tmpLinestring);       
      
-- From how many POLYGON parts the feature is build from
      nGeometries = st_numGeometries(tmpPolygon);
      
-- For all POLYGON parts k = Pk      
      for k in 1..nGeometries LOOP

-- Take the Pk for further investigations      
        Pk := st_GeometryN(tmpPolygon,k);
        
-- Detect the exterior ring of Pk        
        extRingPk = st_exteriorRing(Pk);        
        isHole := FALSE;

-- For all POLYGON i = Pi 
-- does have Pi an inner ring which is an exterior ring of Pk?
-- if yes, Pk is a hole could be ignored for further process
        for i in 1..nGeometries LOOP
          Pi := st_GeometryN(tmpPolygon,i);
          
-- Pi consists of how many inner rings?
          nIntRings := st_NumInteriorRings(Pi);
                                           
-- For all inner ringe of Pi                  
          for n in 1..nIntRings LOOP
            intRingPi := st_interiorRingN(Pi,n);            
                        
-- Check, whether the inner ring of Pi is equal the exterior ring of Pk
-- use ST_SNAPTOGRID due to rounding differences          
            if equals(st_snaptogrid(intRingPi,0.00001),st_snaptogrid(extRingPk,0.00001)) THEN
              isHole := TRUE;
            END IF;
          END LOOP;
        END LOOP;    
        
        if NOT isHole THEN    

-- For the first Geometry
          if outGeom ISNULL THEN
            outGeom := Pk;
          ELSE  
            outGeom := st_Union(outGeom,Pk);
          END IF;  
        END IF;

      END LOOP;  
      RETURN st_multi(outGeom);
    else    
      RETURN st_multi(inGeom);
    END IF;


------------------------------------------------------------------------------
-- Clean Process for LINESTRINGS, self-intersecting parts of linestrings 
-- will be divided into multiparts of the mentioned linestring 
------------------------------------------------------------------------------
  ELSIF (GeometryType(inGeom) = 'LINESTRING' OR GeometryType(inGeom) = 'MULTILINESTRING') THEN
    
-- create nodes at all self-intersecting lines by union the linestrings
-- with the startingpoint of the linestring.  
    outGeom := st_union(st_multi(inGeom),st_pointn(inGeom,1));
    RETURN outGeom;
  ELSE 
    RAISE EXCEPTION 'The input type % is not supported',GeometryType(inGeom);
  END IF;	  
End;$BODY$
  LANGUAGE 'plpgsql' VOLATILE;
