-
-CREATE OR REPLACE FUNCTION _CDB_Feature_Density(reloid REGCLASS, max_z integer)
+CREATE OR REPLACE FUNCTION _CDB_Feature_Density(reloid REGCLASS, nz integer)
 RETURNS FLOAT8
 AS $$
   DECLARE
     fd FLOAT8;
     min_features TEXT;
+    n integer = 4;
+    c FLOAT8;
   BEGIN
   -- TODO: for small total count or extents we could just:
   -- EXECUTE 'SELECT Count(*)/ST_Area(ST_Extent(the_geom_webmercator)) FROM ' || reloid::text || ';' INTO fd;

  -- the area of tiles at level Z: c*c*power(2, -2*z)
  -- with c = CDB_XYZ_Resolution(-8) (earth circumference)
  min_features = '500';

  -- TODO: compute min_z and seed tiles based on reloid extents
  -- so as to have at least N^2 tiles that cover the extents for some N
  -- For the time being we use the root tile as the seed.
  EXECUTE Format('
    WITH RECURSIVE t(x, y, z, e) AS (
      SELECT 0, 0, 0, (
        SELECT count(*) FROM %1$s
          WHERE the_geom_webmercator && CDB_XYZ_Extent(0, 0, 0)
      )
      UNION ALL
      SELECT x*2 + xx, y*2 + yy, z+1, (
        SELECT count(*) FROM %1$s
          WHERE the_geom_webmercator && CDB_XYZ_Extent(x*2 + xx, y*2 + yy, z+1)
      )
      FROM t, (VALUES (0, 0), (0, 1), (1, 1), (1, 0)) AS c(xx, yy)
      WHERE e > %2$s AND z < %3$s
    )
    SELECT MAX(e/ST_Area(CDB_XYZ_Extent(x,y,z))) FROM t where e > 0;
  ', reloid::text, min_features, max_z)
  INTO fd;
  RETURN fd;
  END

AS $$
  DECLARE
    lim FLOAT8 := 500; -- TODO: determine/parameterize this
    nz integer := 4;
    fd FLOAT8;
    c FLOAT8;
  BEGIN
    -- Compute fd as an estimation of the (maximum) number
    -- of features per unit of tile area (in webmercator squared meters)
    SELECT _CDB_Feature_Density(reloid, nz) INTO fd;
    -- lim maximum number of (desiderable) features per tile
    -- we have c = 2*Pi*R = CDB_XYZ_Resolution(-8) (earth circumference)
    -- ta(z): tile area = power(c*power(2,z), 2) = c*c*power(2,2*z)
