/* @file  pmFootprintFindAtPoint.c
 * find footprints in a small image relative to a reference point
 *
 * @author RHL, Princeton & IfA; EAM, IfA
 *
 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
 * @date $Date: 2009-01-27 06:39:38 $
 * Copyright 2006 Institute for Astronomy, University of Hawaii
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <strings.h>
#include <pslib.h>
#include "pmSpan.h"
#include "pmFootprint.h"
#include "pmPeaks.h"

/*
 * A data structure to hold the starting point for a search for pixels above threshold,
 * used by pmFootprintsFindAtPoint
 *
 * We don't want to find this span again --- it's already part of the footprint ---
 * so we set appropriate mask bits
 *
 * EAM : these function were confusingly using "startspan" and "spartspan"  
 * I've rationalized them all to 'startspan'
 */

//
// An enum for what we should do with a Startspan
//
typedef enum {PM_SSPAN_DOWN = 0,	// scan down from this span
	      PM_SSPAN_UP,		// scan up from this span
	      PM_SSPAN_RESTART,		// restart scanning from this span
	      PM_SSPAN_DONE		// this span is processed
} PM_SSPAN_DIR;				// How to continue searching
//
// An enum for mask's pixel values.  We're looking for pixels that are above threshold, and
// we keep extra book-keeping information in the PM_SSPAN_STOP plane.  It's simpler to be
// able to check for 
//
enum {
    PM_SSPAN_INITIAL = 0x0,		// initial state of pixels.
    PM_SSPAN_DETECTED = 0x1,		// we've seen this pixel
    PM_SSPAN_STOP = 0x2			// you may stop searching when you see this pixel
};
//
// The struct that remembers how to [re-]start scanning the image for pixels
//
typedef struct {
    const pmSpan *span;			// save the pixel range
    PM_SSPAN_DIR direction;		// How to continue searching
    bool stop;				// should we stop searching?
} Startspan;

static void startspanFree(Startspan *sspan) {
    psFree((void *)sspan->span);
}

static Startspan *
StartspanAlloc(const pmSpan *span,	// The span in question
	       psImage *mask,		// Pixels that we've already detected
	       const PM_SSPAN_DIR dir	// Should we continue searching towards the top of the image?
    ) {
    Startspan *sspan = psAlloc(sizeof(Startspan));
    psMemSetDeallocator(sspan, (psFreeFunc)startspanFree);

    sspan->span = psMemIncrRefCounter((void *)span);
    sspan->direction = dir;
    sspan->stop = false;
    
    if (mask != NULL) {			// remember that we've detected these pixels
	psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];

	for (int i = 0; i <= span->x1 - span->x0; i++) {
	    mpix[i] |= PM_SSPAN_DETECTED;
	    if (mpix[i] & PM_SSPAN_STOP) {
		sspan->stop = true;
	    }
	}
    }
    
    return sspan;
}

//
// Add a new Startspan to an array of Startspans.  Iff we see a stop bit, return true
//
static bool add_startspan(psArray *startspans, // the saved Startspans
			  const pmSpan *sp, // the span in question
			  psImage *mask, // mask of detected/stop pixels
			  const PM_SSPAN_DIR dir) { // the desired direction to search
    if (dir == PM_SSPAN_RESTART) {
	if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
	    add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
	    return true;
	}
    } else {
	Startspan *sspan = StartspanAlloc(sp, mask, dir);
	if (sspan->stop) {		// we detected a stop bit
	    psFree(sspan);		// don't allocate new span

	    return true;
	} else {
	    psArrayAdd(startspans, 1, sspan);
	    psFree(sspan);		// as it's now owned by startspans
	}
    }

    return false;
}

/************************************************************************************************************/
/*
 * Search the image for pixels above threshold, starting at a single Startspan.
 * We search the array looking for one to process; it'd be better to move the
 * ones that we're done with to the end, but it probably isn't worth it for
 * the anticipated uses of this routine.
 *
 * This is the guts of pmFootprintsFindAtPoint
 */
static bool do_startspan(pmFootprint *fp, // the footprint that we're building
			 const psImage *img, // the psImage we're working on
			 psImage *mask, // the associated masks
			 const float threshold,	// Threshold
			 psArray *startspans) {	// specify which span to process next
    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
    psImageMaskType *maskRow = NULL;		//  masks's row pointer
    
    const int row0 = img->row0;
    const int col0 = img->col0;
    const int numRows = img->numRows;
    const int numCols = img->numCols;
    
    /********************************************************************************************************/
    
    Startspan *sspan = NULL;
    for (int i = 0; i < startspans->n; i++) {
	sspan = startspans->data[i];
	if (sspan->direction != PM_SSPAN_DONE) {
	    break;
	}
	if (sspan->stop) {
	    break;
	}
    }
    if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process
	return false;
    }
    if (sspan->stop) {			// they don't want any more spans processed
	return false;
    }
    /*
     * Work
     */
    const PM_SSPAN_DIR dir = sspan->direction;
    /*
     * Set initial span to the startspan
     */
    int x0 = sspan->span->x0 - col0, x1 = sspan->span->x1 - col0;
    /*
     * Go through image identifying objects
     */
    int nx0, nx1 = -1;			// new values of x0, x1
    const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row
    bool stop = false;			// should I stop searching for spans?

    for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) {
	imgRowF32 = img->data.F32[i];	// only one of
	imgRowS32 = img->data.S32[i];	//      these is valid!
	maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i];
	//
	// Search left from the pixel diagonally to the left of (i - di, x0). If there's
	// a connected span there it may need to grow up and/or down, so push it onto
	// the stack for later consideration
	//
	nx0 = -1;
	for (int j = x0 - 1; j >= -1; j--) {
	    double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
	    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
		if (j < x0 - 1) {	// we found some pixels above threshold
		    nx0 = j + 1;
		}
		break;
	    }
	}

	if (nx0 < 0) {			// no span to the left
	    nx1 = x0 - 1;		// we're going to resume searching at nx1 + 1
	} else {
	    //
	    // Search right in leftmost span
	    //
	    //nx1 = 0;			// make gcc happy
	    for (int j = nx0 + 1; j <= numCols; j++) {
		double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
		if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
		    nx1 = j - 1;
		    break;
		}
	    }
	    
	    const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
	    
	    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
		stop = true;
		break;
	    }
	}
	//
	// Now look for spans connected to the old span.  The first of these we'll
	// simply process, but others will have to be deferred for later consideration.
	//
	// In fact, if the span overhangs to the right we'll have to defer the overhang
	// until later too, as it too can grow in both directions
	//
	// Note that column numCols exists virtually, and always ends the last span; this
	// is why we claim below that sx1 is always set
	//
	bool first = false;		// is this the first new span detected?
	for (int j = nx1 + 1; j <= x1 + 1; j++) {
	    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
	    if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
		int sx0 = j++;		// span that we're working on is sx0:sx1
		int sx1 = -1;		// We know that if we got here, we'll also set sx1
		for (; j <= numCols; j++) {
		    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
		    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
			sx1 = j;
			break;
		    }
		}
		assert (sx1 >= 0);

		const pmSpan *sp;
		if (first) {
		    if (sx1 <= x1) {
			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
			    stop = true;
			    break;
			}
		    } else {		// overhangs to right
			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
			    stop = true;
			    break;
			}
			sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
			if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
			    stop = true;
			    break;
			}
		    }
		    first = false;
		} else {
		    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
		    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
			stop = true;
			break;
		    }
		}
	    }
	}

	if (stop || first == false) {	// we're done
	    break;
	}

	x0 = nx0; x1 = nx1;
    }
    /*
     * Cleanup
     */

    sspan->direction = PM_SSPAN_DONE;
    return stop ? false : true;
}

/*
 * Go through an image, starting at (row, col) and assembling all the pixels
 * that are connected to that point (in a chess kings-move sort of way) into
 * a pmFootprint.
 *
 * This is much slower than pmFootprintsFind if you want to find lots of
 * footprints, but if you only want a small region about a given point it
 * can be much faster
 *
 * N.b. The returned pmFootprint is not in "normal form"; that is the pmSpans
 * are not sorted by increasing y, x0, x1.  If this matters to you, call
 * pmFootprintNormalize()
 */
pmFootprint *
pmFootprintsFindAtPoint(const psImage *img,	// image to search
		       const float threshold,	// Threshold
		       const psArray *peaks, // array of peaks; finding one terminates search for footprint
		       int row, int col) { // starting position (in img's parent's coordinate system)
   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;
/*
 * Is point in image, and above threshold?
 */
   row -= row0; col -= col0;
   if (row < 0 || row >= numRows ||
       col < 0 || col >= numCols) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",
		row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
       return NULL;
   }

   double pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col];
   if (pixVal < threshold) {
       return pmFootprintAlloc(0, img);
   }
   
   pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img);
/*
 * We need a mask for two purposes; to indicate which pixels are already detected,
 * and to store the "stop" pixels --- those that, once reached, should stop us
 * looking for the rest of the pmFootprint.  These are generally set from peaks.
 */
   psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
   P_PSIMAGE_SET_ROW0(mask, row0);
   P_PSIMAGE_SET_COL0(mask, col0);
   psImageInit(mask, PM_SSPAN_INITIAL);
   //
   // Set stop bits from peaks list
   //
   assert (peaks == NULL || peaks->n == 0 || psMemCheckPeak(peaks->data[0]));
   if (peaks != NULL) {
       for (int i = 0; i < peaks->n; i++) {
	   pmPeak *peak = peaks->data[i];
	   mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
       }
   }
/*
 * Find starting span passing through (row, col)
 */
   psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search
   
   imgRowF32 = img->data.F32[row];	// only one of
   imgRowS32 = img->data.S32[row];	//      these is valid!
   psImageMaskType *maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[row];
   {
       int i;
       for (i = col; i >= 0; i--) {
	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
	       break;
	   }
       }
       int i0 = i;
       for (i = col; i < numCols; i++) {
	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
	       break;
	   }
       }
       int i1 = i;
       const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1);

       (void)add_startspan(startspans, sp, mask, PM_SSPAN_RESTART);
   }
   /*
    * Now workout from those Startspans, searching for pixels above threshold
    */
   while (do_startspan(fp, img, mask, threshold, startspans)) continue;
   /*
    * Cleanup
    */
   psFree(mask);
   psFree(startspans);			// restores the image pixel

   return fp;				// pmFootprint really
}
