# include "psphotInternal.h"
void *memset (void *s, int c, size_t n);

/**** find footprints in an image (fast on large scale images) ****/

// XXX EAM : why use WSPAN in here rather than pmSpan?
typedef struct {			/* run-length code for part of object*/
   int id;				/* ID for object */
   int y;				/* Row wherein WSPAN dwells */
   int x0, x1;				/* inclusive range of columns */
} WSPAN;

/*
 * comparison function for qsort; sort by ID then row
 */
static int
compar(const void *va, const void *vb)
{
   const WSPAN *a = va;
   const WSPAN *b = vb;

   if(a->id < b->id) {
      return(-1);
   } else if(a->id > b->id) {
      return(1);
   } else {
      return(a->y - b->y);
   }
}

/*
 * Follow a chain of aliases, returning the final resolved value.
 */
static int
resolve_alias(const int *aliases,	/* list of aliases */
	      int id)			/* alias to look up */
{
   int resolved = id;			/* resolved alias */

   while(id != aliases[id]) {
      resolved = id = aliases[id];
   }

   return(resolved);
}

/*
 * Go through an image, finding sets of connected pixels above threshold
 * and assembling them into pmFootprints;  the resulting set of objects
 * is returned as a psArray
 */
psArray *
pmFootprintsFind(const psImage *img,	// image to search
		 const float threshold,	// Threshold
		 const int npixMin)	// minimum number of pixels in an acceptable pmFootprint
{
   int i0;				/* initial value of i */
   int id;				/* object ID */
   int in_span;				/* object ID of current WSPAN */
   int nspan = 0;			/* number of spans */
   int nobj = 0;			/* number of objects found */
   int x0 = 0;			        /* unpacked from a WSPAN */
   int *tmp;				/* used in swapping idc/idp */

   assert(img != NULL);

   bool F32 = false;			// is this an F32 image?
   if (img->type.type == PS_TYPE_F32) {
       F32 = true;
   } else if (img->type.type == PS_TYPE_S32) {
       F32 = false;
   } else {				// N.b. You can't trivially add more cases here; F32 is just a bool
       psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
       return NULL;
   }
   psF32 *imgRowF32 = NULL;		// row pointer if F32
   psS32 *imgRowS32 = NULL;		//  "   "   "  "  !F32
   
   const int row0 = img->row0;
   const int col0 = img->col0;
   const int numRows = img->numRows;
   const int numCols = img->numCols;
/*
 * Storage for arrays that identify objects by ID. We want to be able to
 * refer to idp[-1] and idp[numCols], hence the (numCols + 2)
 */
   int *id_s = psAlloc(2*(numCols + 2)*sizeof(int));
   memset(id_s, '\0', 2*(numCols + 2)*sizeof(int)); assert(id_s[0] == 0);
   int *idc = id_s + 1;			// object IDs in current/
   int *idp = idc + (numCols + 2);	//                       previous row

   int size_aliases = 1 + numRows/20;	// size of aliases[] array
   int *aliases = psAlloc(size_aliases*sizeof(int)); // aliases for object IDs

   int size_spans = 1 + numRows/20;	// size of spans[] array
   WSPAN *spans = psAlloc(size_spans*sizeof(WSPAN)); // row:x0,x1 for objects
/*
 * Go through image identifying objects
 */
   for (int i = 0; i < numRows; i++) {
      int j;
      tmp = idc; idc = idp; idp = tmp;	/* swap ID pointers */
      memset(idc, '\0', numCols*sizeof(int));
      
      imgRowF32 = img->data.F32[i];	// only one of
      imgRowS32 = img->data.S32[i];	//      these is valid!

      in_span = 0;			/* not in a span */
      for (j = 0; j < numCols; j++) {
	 double pixVal = F32 ? imgRowF32[j] : imgRowS32[j];
	 if (pixVal < threshold) {
	    if (in_span) {
	       if(nspan >= size_spans) {
		  size_spans *= 2;
		  spans = psRealloc(spans, size_spans*sizeof(WSPAN));
	       }
	       spans[nspan].id = in_span;
	       spans[nspan].y = i;
	       spans[nspan].x0 = x0;
	       spans[nspan].x1 = j - 1;
	       
	       nspan++;

	       in_span = 0;
	    }
	 } else {			/* a pixel to fix */
	    if(idc[j - 1] != 0) {
	       id = idc[j - 1];
	    } else if(idp[j - 1] != 0) {
	       id = idp[j - 1];
	    } else if(idp[j] != 0) {
	       id = idp[j];
	    } else if(idp[j + 1] != 0) {
	       id = idp[j + 1];
	    } else {
	       id = ++nobj;

	       if(id >= size_aliases) {
		  size_aliases *= 2;
		  aliases = psRealloc(aliases, size_aliases*sizeof(int));
	       }
	       aliases[id] = id;
	    }

	    idc[j] = id;
	    if(!in_span) {
	       x0 = j; in_span = id;
	    }
/*
 * Do we need to merge ID numbers? If so, make suitable entries in aliases[]
 */
	    if(idp[j + 1] != 0 && idp[j + 1] != id) {
	       aliases[resolve_alias(aliases, idp[j + 1])] =
						    resolve_alias(aliases, id);
	       
	       idc[j] = id = idp[j + 1];
	    }
	 }
      }

      if(in_span) {
	 if(nspan >= size_spans) {
	    size_spans *= 2;
	    spans = psRealloc(spans, size_spans*sizeof(WSPAN));
	 }

	 assert(nspan < size_spans);	/* we checked for space above */
	 spans[nspan].id = in_span;
	 spans[nspan].y = i;
	 spans[nspan].x0 = x0;
	 spans[nspan].x1 = j - 1;
	 
	 nspan++;
      }
   }

   psFree(id_s);
/*
 * Resolve aliases; first alias chains, then the IDs in the spans
 */
   for (int i = 0; i < nspan; i++) {
      spans[i].id = resolve_alias(aliases, spans[i].id);
   }

   psFree(aliases);
/*
 * Sort spans by ID, so we can sweep through them once
 */
   if(nspan > 0) {
      qsort(spans, nspan, sizeof(WSPAN), compar);
   }
/*
 * Build pmFootprints from the spans
 */
   psArray *footprints = psArrayAlloc(nobj);
   int n = 0;			// number of pmFootprints

   if(nspan > 0) {
      id = spans[0].id;
      i0 = 0;
      for (int i = 0; i <= nspan; i++) {	/* nspan + 1 to catch last object */
	 if(i == nspan || spans[i].id != id) {
	    pmFootprint *fp = pmFootprintAlloc(i - i0, img);
	    
	    for(; i0 < i; i0++) {
		pmFootprintAddSpan(fp, spans[i0].y + row0,
				   spans[i0].x0 + col0, spans[i0].x1 + col0);
	    }

	    if (fp->npix < npixMin) {
	       psFree(fp);
	    } else {
	       footprints->data[n++] = fp;
	    }
	 }
	 
	 id = spans[i].id;
      }
   }

   footprints = psArrayRealloc(footprints, n);
/*
 * clean up
 */
   psFree(spans);

   return footprints;
}
