-- Function: rdat_infer_common_filename(name, name, name, name, text)

-- DROP FUNCTION rdat_infer_common_filename(name, name, name, name, text);

CREATE OR REPLACE FUNCTION rdat_infer_common_filename(table_name name, schema_name name, raster_column name, filename_column name, filename text)
  RETURNS raster AS
$BODY$DECLARE
   num_tile_geometries integer := 0 ;
   interior_tile RECORD ;
   dummy RECORD ; 
   theTable regclass := (schema_name || '.' || table_name)::regclass;
   num_results int ; 
   gt_probe RECORD ; 
   scaleX real ; 
   scaleY real ; 
   skewX  real ; 
   skewY  real ; 
   srid int ; 
   thresh_ScaleX real := 1e-6 ; 
   thresh_ScaleY real := 1e-6 ; 
   thresh_SkewX real := 1e-6 ; 
   thresh_SkewY real := 1e-6 ;
   arbitrary_tile raster ; 
   misaligned_tiles int ; 
   dataset_raster raster ; 
   dataset_envelope RECORD ; 
BEGIN 

--
-- Check to ensure that the filename exists in the table.
--
EXECUTE 'SELECT 1 FROM '
   || theTable 
   || ' WHERE '
   || filename_column || '= $1 LIMIT 1'
  INTO dummy 
  USING filename ; 
GET DIAGNOSTICS num_results = ROW_COUNT ; 
IF num_results = 0 THEN
   RAISE EXCEPTION 'Cannot find filename %.', filename ; 
END IF ;

--
-- Do a basic check to ensure that the tile geometry 
-- (integer width and height) are the same. 
--
EXECUTE 'SELECT count(*) FROM '
   || '(SELECT DISTINCT ST_Width(' || raster_column 
   || '), ST_Height(' || raster_column || ') '
   || ' FROM ' 
   || theTable 
   || ' WHERE ' || filename_column || ' = $1) as dummy ' 
   INTO num_tile_geometries
   USING filename ;
GET DIAGNOSTICS num_results = ROW_COUNT ; 

IF num_results = 0 THEN
  RAISE EXCEPTION 'Image % not found.', filename ; 
END IF;

IF num_tile_geometries > 1 THEN 
  RAISE EXCEPTION 'Image % does not appear to be regularly tiled.', filename ;
END IF;

--
-- So the image was found and there's only one tile geometry.
-- Get the tile geometry.
--
EXECUTE 'SELECT ST_Width(' || raster_column 
   || ') as width, ST_Height(' || raster_column 
   || ') as height '
   || 'FROM ' 
   || theTable 
   || ' WHERE ' || filename_column || ' = $1 ' 
   || ' LIMIT 1'
  INTO interior_tile
  USING filename ;

--
-- Now check to ensure that the geotransform parameters are 
-- all the same. Since these are floating point, we check that 
-- max() - min() for each parameter does not exceed a specified 
-- threshold.
--

-- ScaleX
EXECUTE 'SELECT '
   || 'max(ST_ScaleX(' || raster_column || ')) as scalex_max, '
   || 'min(ST_ScaleX(' || raster_column || ')) as scalex_min '
   || 'FROM '
   || theTable
   || ' WHERE ' || filename_column || ' = $1 ' 
  INTO  gt_probe
  USING filename ;
IF (gt_probe.scalex_max - gt_probe.scalex_min) > thresh_ScaleX THEN 
  RAISE EXCEPTION 'All tiles must have the same x scale geotransform parameter';
END IF ; 
scaleX := (gt_probe.scalex_max + gt_probe.scalex_min) / 2.0 ; 
    
-- ScaleY
EXECUTE 'SELECT '
   || 'max(ST_ScaleY(' || raster_column || ')) as scaley_max, '
   || 'min(ST_ScaleY(' || raster_column || ')) as scaley_min '
   || 'FROM '
   || theTable
   || ' WHERE ' || filename_column || ' = $1 ' 
  INTO  gt_probe
  USING filename ;
IF (gt_probe.scaley_max - gt_probe.scaley_min) > thresh_ScaleY THEN 
  RAISE EXCEPTION 'All tiles must have the same y scale geotransform parameter';
END IF ; 
scaleY := (gt_probe.scaley_max + gt_probe.scaley_min) / 2.0 ; 

-- SkewX
EXECUTE 'SELECT '
   || 'max(ST_SkewX(' || raster_column || ')) as skewx_max, '
   || 'min(ST_SkewX(' || raster_column || ')) as skewx_min '
   || 'FROM '
   || theTable
   || ' WHERE ' || filename_column || ' = $1 ' 
  INTO  gt_probe
  USING filename ;
IF (gt_probe.skewx_max - gt_probe.skewx_min) > thresh_SkewX THEN 
  RAISE EXCEPTION 'All tiles must have the same x skew geotransform parameter';
END IF ; 
skewX := (gt_probe.skewx_max + gt_probe.skewx_min) / 2.0 ; 

-- SkewY
EXECUTE 'SELECT '
   || 'max(ST_SkewY(' || raster_column || ')) as skewy_max, '
   || 'min(ST_SkewY(' || raster_column || ')) as skewy_min '
   || 'FROM '
   || theTable
   || ' WHERE ' || filename_column || ' = $1 ' 
  INTO  gt_probe
  USING filename ;
IF (gt_probe.skewy_max - gt_probe.skewy_min) > thresh_SkewY THEN 
  RAISE EXCEPTION 'All tiles must have the same x scale geotransform parameter';
END IF ; 
skewY := (gt_probe.skewy_max + gt_probe.skewy_min) / 2.0 ; 

RAISE NOTICE 'Tiles are %x% pixels, (%,%,%,%)',
         interior_tile.width, interior_tile.height,
         scaleX, scaleY, skewX, skewY ; 

--
-- All of the candidate tiles must be in the same SRID.
-- 
EXECUTE 'SELECT count(*) FROM '
   || '(SELECT DISTINCT ST_Srid(' || raster_column 
   || ') FROM ' || theTable
   || ' WHERE ' || filename_column || ' = $1 ) as dummy'
   INTO num_results
   USING filename ; 
IF num_results <> 1 THEN 
   RAISE EXCEPTION 'All tiles must have the same SRID.' ; 
END IF;
EXECUTE 'SELECT ST_Srid(' || raster_column 
   || ') FROM ' || theTable
   || ' WHERE ' || filename_column || ' = $1 '
   || 'LIMIT 1'
   INTO srid
   USING filename ;  

--
-- Check to ensure that the offset vector 
-- (aka "UpperLeft") is always shifted by 
-- an integer number of tiles
--
   
EXECUTE 'SELECT ' || raster_column
   || ' FROM ' || theTable
   || ' WHERE ' || filename_column || ' = $1'
   || ' LIMIT 1'
   INTO arbitrary_tile
   USING filename ;

--
-- Convert the offset vector to each tile (UpperLeft) 
-- into zero-based pixel indices, then modulo the 
-- width/height of the tile. 
--
EXECUTE 'SELECT count(*) as aligned FROM ' || theTable
   || ' WHERE ' || filename_column || ' = $1 and '
   || ' ((ST_World2RasterCoordX($2, '
   || '  ST_UpperLeftX(' || raster_column 
   || '  ), ST_UpperLeftY(' || raster_column 
   || '  ))-1) % $3 <> 0) and ' 
   || ' ((ST_World2RasterCoordY($2, '
   || '  ST_UpperLeftX(' || raster_column 
   || '  ), ST_UpperLeftY(' || raster_column 
   || '  ))-1) % $4 <> 0)' 
   INTO misaligned_tiles
   USING filename, arbitrary_tile, interior_tile.width, interior_tile.height ; 
IF misaligned_tiles <> 0 THEN 
   RAISE EXCEPTION 'All tiles must be offset by an integer multiple of the tile size. %', misaligned_tiles ; 
END IF ; 


--
-- Now calculate the envelope in (i,j) space of the 
-- overall raster dataset (incl all tiles) 
-- 

EXECUTE 'SELECT ' 
   || 'min(ST_World2RasterCoordX($2'
      || ',ST_UpperLeftX(' || raster_column
      || '),ST_UpperLeftY(' || raster_column
      || '))) as min_i, '
   || 'max(ST_World2RasterCoordX($2'
      || ',ST_UpperLeftX(' || raster_column
      || '),ST_UpperLeftY(' || raster_column
      || '))) as max_i, '
   || 'min(ST_World2RasterCoordY($2'
      || ',ST_UpperLeftX(' || raster_column
      || '),ST_UpperLeftY(' || raster_column
      || '))) as min_j, '
   || 'max(ST_World2RasterCoordY($2' 
      || ',ST_UpperLeftX(' || raster_column
      || '),ST_UpperLeftY(' || raster_column
      || '))) as max_j '
   || 'FROM ' || theTable
   || ' WHERE ' || filename_column || ' = $1 ' 
   INTO dataset_envelope
   USING filename, arbitrary_tile ; 

--
-- Calculate the offset vector to the origin of
-- the overall dataset (currently named "UpperLeft")
--
EXECUTE 'SELECT ST_Raster2WorldCoordX($1, $2, $3) as ulx, '
   || 'ST_Raster2WorldCoordY($1, $2, $3) as uly '
   INTO gt_probe 
   USING arbitrary_tile, dataset_envelope.min_i, dataset_envelope.min_j;

--
-- Now construct the dataset raster to serve as a 
-- tile index.
SELECT ST_MakeEmptyRaster(dataset_envelope.max_i - dataset_envelope.min_i + 1,
    dataset_envelope.max_j - dataset_envelope.min_j + 1,
    gt_probe.ulx, gt_probe.uly,
    scaleX, scaleY, skewX, skewY, srid) 
    INTO dataset_raster ; 
 
RETURN dataset_raster ;    

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_infer_common_filename(name, name, name, name, text)
  OWNER TO postgres;
COMMENT ON FUNCTION rdat_infer_common_filename(name, name, name, name, text) IS 'Attempts to reconstruct a single image description from a raster column paired with a filename column and a named numeric primary key column. The raster column, the filename column, and the primary key column must all exist in the same catalog/schema/table. The following constraints apply to tiles originating from the same image: 

+ The scale and skew parameters of the geotransform for each tile must be identical.
+ The tile size (width and height) is the same for every tile in the image, except for the "last tile in the row" and the "last tile in the column".
++ Every "interior tile" must be the same size.
++ Every "last tile in the row" must be the same size, and must be smaller than an interior tile.
++ Every "last tile in the column" must be the same size, and must be smaller than an interior tile.
+ Each tile must have the same number of bands.
+ Each tile''s bands must have the same data types.
+ The offset (x|y) parameters must be different by an integer number of interior tiles.';


-- Function: rdat_sync_dataset_using_filename(raster, name, name, name, name, name, text)

-- DROP FUNCTION rdat_sync_dataset_using_filename(raster, name, name, name, name, name, text);

CREATE OR REPLACE FUNCTION rdat_sync_dataset_using_filename(index raster, table_name name, schema_name name, raster_column name, key_column name, filename_column name, filename text)
  RETURNS raster AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  upperLeft geometry ; 
  newIndex raster ; 
  key bigint ; 
  tile raster ;
  
BEGIN
  -- Start over from scratch
  SELECT ST_MakeEmptyRaster(index)  INTO newIndex ; 
  SELECT ST_AddBand(newIndex, '32BSI', -1, -1) INTO newIndex;

  -- store the key to the "tile table" in the new index raster.
  FOR tile, key IN EXECUTE 'SELECT ' || raster_column || ',' || key_column
            || ' FROM ' || theTable 
            || ' WHERE ' || filename_column || ' = $1 '
           USING filename
  LOOP

    SELECT ST_SetSRID(ST_Point(ST_UpperLeftX(tile), ST_UpperLeftY(tile)),ST_Srid(index))
        INTO upperLeft ; 

    SELECT ST_SetVaue(newIndex, upperLeft, key::double precision) ; 

  END LOOP;    

  RETURN newIndex ;  
END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_sync_dataset_using_filename(raster, name, name, name, name, name, text)
  OWNER TO postgres;

-- Function: rdat_read_block(raster, integer, integer, name, name, name, name)

-- DROP FUNCTION rdat_read_block(raster, integer, integer, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_read_block(index raster, block_i integer, block_j integer, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS raster AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  keyval bigint ; 
  block raster ; 
  
BEGIN
  -- find if this block already exists.
  SELECT ST_Value(index, block_i, block_j)::bigint INTO keyval ; 

  block := NULL ; 
  IF keyval IS NOT NULL THEN 
    EXECUTE 'SELECT ' || raster_column 
        || ' FROM ' || theTable
        || ' WHERE ' || key_column || ' = $1'
    INTO block
    USING keyval ;  
  END IF ; 

  RETURN block ;
END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_read_block(raster, integer, integer, name, name, name, name)
  OWNER TO postgres;

-- Function: rdat_write_block(raster, raster, integer, integer, name, name, name, name)

-- DROP FUNCTION rdat_write_block(raster, raster, integer, integer, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_write_block(index raster, block raster, block_i integer, block_j integer, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS void AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  theSequence regclass ; 
  keyval bigint ; 
  
BEGIN
  -- find if this block already exists.
  SELECT ST_Value(index, block_i, block_j)::bigint INTO keyval ; 

  IF keyval IS NULL THEN 
    SELECT pg_get_serial_sequence(theTable, key_column) INTO theSequence ; 
    SELECT nextval(theSequence) INTO keyval ; 
    EXECUTE 'INSERT INTO ' || theTable 
      || ' (' || key_column || ',' || raster_column || ') VALUES ($1,$2)'
      USING keyval, block ;
  ELSE 
    EXECUTE 'UPDATE ' || theTable
      || ' SET ' || raster_column || ' = $1'
      || ' WHERE ' || key_column || ' = $2'
    USING block, keyval ;  
  END IF ; 
  
END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_write_block(raster, raster, integer, integer, name, name, name, name)
  OWNER TO postgres;

-- Function: rdat_setvalue(raster, geometry, integer, double precision, name, name, name, name)

-- DROP FUNCTION rdat_setvalue(raster, geometry, integer, double precision, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_setvalue(index raster, pt geometry, band integer, value double precision, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS void AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  keyval bigint ;
  
BEGIN
  -- find if this block exists.
  SELECT ST_Value(index, pt)::bigint INTO keyval ; 

  IF keyval IS NULL THEN
    -- to truly duplicate the behavior of ST_SetValue, we need to create
    -- the missing block here, then set the value.
    -- But to do that, we'll need lists of pixel types and nodata values
    -- for all the bands.
    RETURN ; 
  END IF ;  

  EXECUTE 'SELECT ST_SetValue(' || raster_column 
     || ', $2, $3) '
     || ' FROM ' || theTable
     || ' WHERE ' || key_column || ' = $1'
  INTO value
  USING keyval, band, pt ;  

END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_setvalue(raster, geometry, integer, double precision, name, name, name, name)
  OWNER TO postgres;

-- Function: rdat_setvalue(raster, integer, integer, integer, double precision, integer, integer, name, name, name, name)

-- DROP FUNCTION rdat_setvalue(raster, integer, integer, integer, double precision, integer, integer, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_setvalue(index raster, i integer, j integer, band integer, value double precision, tile_width integer, tile_height integer, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS void AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  keyval bigint ;
  block_i int ; 
  block_j int ; 
  inblock_i int ; 
  inblock_j int ; 
  block raster ; 
  
BEGIN
  block_i := ((i-1)/tile_width)+1 ; 
  block_j := ((j-1)/tile_height)+1 ; 
  inblock_i := ((i-1)%tile_width) + 1 ; 
  inblock_j := ((j-1)%tile_height) + 1; 

  -- find if this block exists.
  SELECT ST_Value(index, block_i, block_j)::bigint INTO keyval ; 

  IF keyval IS NULL THEN
    -- to truly duplicate the behavior of ST_SetValue, we need to create
    -- the missing block here, then set the value.
    -- But to do that, we'll need lists of pixel types and nodata values
    -- for all the bands.
    RETURN ; 
  END IF ;  

  EXECUTE 'SELECT ST_SetValue(' || raster_column 
     || ', $2, $3, $4) '
     || ' FROM ' || theTable
     || ' WHERE ' || key_column || ' = $1'
  INTO value
  USING keyval, band, inblock_i, inblock_j ;  

END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_setvalue(raster, integer, integer, integer, double precision, integer, integer, name, name, name, name)
  OWNER TO postgres;

-- Function: rdat_value(raster, integer, integer, integer, integer, integer, name, name, name, name)

-- DROP FUNCTION rdat_value(raster, integer, integer, integer, integer, integer, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_value(index raster, i integer, j integer, band integer, tile_width integer, tile_height integer, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS double precision AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  keyval bigint ;
  value double precision ;  
  block_i int ; 
  block_j int ; 
  inblock_i int ; 
  inblock_j int ; 
  block raster ; 
  
BEGIN
  block_i := ((i-1)/tile_width)+1 ; 
  block_j := ((j-1)/tile_height)+1 ; 
  inblock_i := ((i-1)%tile_width) + 1 ; 
  inblock_j := ((j-1)%tile_height) + 1; 

  -- find if this block exists.
  SELECT ST_Value(index, block_i, block_j)::bigint INTO keyval ; 

  IF keyval IS NULL THEN
    -- to truly duplicate the behavior of ST_Value, we need to return the 
    -- nodata value here. But to do that, we'll need a list of nodata values
    -- for all the bands.
    RETURN NULL ; 
  END IF ;  

  EXECUTE 'SELECT ST_Value(' || raster_column 
     || ', $2, $3, $4) '
     || ' FROM ' || theTable
     || ' WHERE ' || key_column || ' = $1'
  INTO value
  USING keyval, band, inblock_i, inblock_j ;  

  RETURN value ;
END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_value(raster, integer, integer, integer, integer, integer, name, name, name, name)
  OWNER TO postgres;

-- Function: rdat_value(raster, geometry, integer, name, name, name, name)

-- DROP FUNCTION rdat_value(raster, geometry, integer, name, name, name, name);

CREATE OR REPLACE FUNCTION rdat_value(index raster, pt geometry, band integer, table_name name, schema_name name, raster_column name, key_column name)
  RETURNS double precision AS
$BODY$
DECLARE
  theTable regclass := ( schema_name || '.' || table_name) :: regclass ;
  keyval bigint ;
  value double precision ;  
  block raster ; 
  
BEGIN
  -- find if this block exists.
  SELECT ST_Value(index, pt)::bigint INTO keyval ; 

  IF keyval IS NULL THEN
    -- to truly duplicate the behavior of ST_Value, we need to return the 
    -- nodata value here. But to do that, we'll need a list of nodata values
    -- for all the bands.
    RETURN NULL ; 
  END IF ;  

  EXECUTE 'SELECT ST_Value(' || raster_column 
     || ', $2, $3) '
     || ' FROM ' || theTable
     || ' WHERE ' || key_column || ' = $1'
  INTO value
  USING keyval, band, pt ;  

  RETURN value ;
END ;

$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
ALTER FUNCTION rdat_value(raster, geometry, integer, name, name, name, name)
  OWNER TO postgres;
