flst.c

image processing including fourier,wavelet,segmentation etc.

/*----------------------------- MegaWave Module -----------------------------*/
/* mwcommand 
  name = {flst}; 
  version = {"2.0"}; 
  author = {"Pascal Monasse"}; 
  function = {"Fast Level Sets Transform of an image"}; 
  usage = { 
    'a': min_area-> pMinArea "argument of the grain filter",
    image-> pImageInput  "Input fimage",
    tree<- pTree "The tree of shapes"
    }; 
*/ 

#include <assert.h>
#include "mw.h"

/* Optimization parameters. Probably should depend on image size, but these
values seem good for most images. */
#define INIT_MAX_AREA 10
#define STEP_MAX_AREA 8

/* A 'local configuration' of the pixel is the logical 'or' of these values,
stored in one byte. Corresponding bit is set if the neighbor is in region */
#define EAST         1
#define NORTH_EAST   2
#define NORTH        4
#define NORTH_WEST   8
#define WEST        16
#define SOUTH_WEST  32
#define SOUTH       64
#define SOUTH_EAST 128

/* Gives for each configuration of the neighborhood around the pixel the number
of new cc's of the complement created (sometimes negative, since
cc's can be deleted), when the pixel is appended to the region.
Configurations are stored on one byte. 
tabPatterns[0]: set in 4-connectedness and complement in 8-connectedness.
tabPatterns[1]: set in 8-connectedness and complement in 4-connectedness. */
static char tabPatterns[2][256] =
{{  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1,
    1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
    0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0,
    1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0,
    0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0,
    0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,-1},
 {  0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
    0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
    0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1,
    0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
    0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1,
    0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1,
    0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1,
    1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1,
    0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0,
    0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1}};

/* The structure representing the neighborhood of a region. It is a priority
queue organized as a heap (a binary tree, where each node has a key larger,
resp. smaller, than its two children), stored in an array.
Index 1 is the root, its children are at index 2 and 3.
Therefore, the parent of node at k is at k/2 and its children at 2*k and 2*k+1.
Index 0 is only a temporary buffer to swap two nodes. */
enum TypeOfTree {AMBIGUOUS, MIN, MAX, INVALID};
typedef struct {
  struct point_plane point; /* Neighbor pixel */
  float value; /* Its gray level */
} Neighbor;
typedef struct
{
  Neighbor* tabPoints; /* The array of neighbors, organized as a binary tree */
  int iNbPoints; /* The size of the previous arrays */
  enum TypeOfTree type; /* max- or min- oriented heap? */
  float otherBound; /* Min gray level if max-oriented, max if min-oriented */
} Neighborhood;

/* Structure to find connections between shapes. This is used when a monotone
section is extracted. The goal is to find the parent of its largest shape. The
gray level of the parent is known, so as a point of the parent, since we use
in fact an image of connections */
typedef struct {
  Shape shape; /* Largest shape of the monotone section */
  float level; /* Connection level */
} Connection;

/* Well, here are global variables. Ugly, but avoids to weigh the code,
since they are used almost everywhere */
static int iWidth, iHeight;
static int iMinArea, iMaxArea, iAreaImage, iHalfAreaImage, iPerimeterImage;
static int iExploration;  /* Used to avoid reinitializing images */
static Point_plane tabPointsInShape;
static int** tabtabVisitedNeighbors; /* Exterior boundary */
static Shapes pGlobalTree;
static int iAtBorder; /* #points meeting at border of image */

#define ARRAY_2D_ALLOC(array, iWidth, iHeight, type) \
(array) = (type**) malloc((iHeight)*sizeof(type*)); \
if((array) == NULL) mwerror(FATAL, 1, "ARRAY_2D_ALLOC --> memory"); \
(array)[0] = (type*) calloc((iWidth)*(iHeight), sizeof(type)); \
if((array)[0] == NULL) mwerror(FATAL, 1, "ARRAY_2D_ALLOC --> memory"); \
for(i = 1; i < (iHeight); i++) \
  (array)[i] = (array)[i-1]+(iWidth);

/* ------------------------------------------------------------------------
   --------- Functions to manage the neighborhood structure ---------------
   ------------------------------------------------------------------------ */

/* Reinitialise the neighborhood, so that it will be used for a new region */
void reinit_neighborhood(pNeighborhood, type)
Neighborhood* pNeighborhood;
enum TypeOfTree type;
{
  pNeighborhood->iNbPoints = 0;
  pNeighborhood->type = type;
}

/* To allocate the structure representing the neighborhood of a region */
void init_neighborhood(pNeighborhood, iMaxArea)
Neighborhood* pNeighborhood;
int iMaxArea;
{
  iMaxArea = 4*(iMaxArea+1);
  if(iMaxArea > iWidth*iHeight)
    iMaxArea = iWidth*iHeight;

  pNeighborhood->tabPoints = (Neighbor*) malloc((iMaxArea+1)*sizeof(Neighbor));
  if(pNeighborhood->tabPoints == NULL)
    mwerror(FATAL, 1, "init_neighborhood --> neighbors allocation error\n");
  reinit_neighborhood(pNeighborhood, AMBIGUOUS);
}

/* Free the structure representing the neighborhood of a region */
void free_neighborhood(pNeighborhood)
Neighborhood* pNeighborhood;
{
  free(pNeighborhood->tabPoints);
}

#define ORDER_MAX(k,l) (tabPoints[k].value > tabPoints[l].value)
#define ORDER_MIN(k,l) (tabPoints[k].value < tabPoints[l].value)
#define SWAP(k,l) tabPoints[0] = tabPoints[k]; \
                  tabPoints[k] = tabPoints[l]; \
                  tabPoints[l] = tabPoints[0];
/* Put the last neighbor at a position so that we fix the heap */
void fix_up(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor *tabPoints = pNeighborhood->tabPoints;
  int k = pNeighborhood->iNbPoints, l;
  
  if(pNeighborhood->type == MAX)
    while(k > 1 && ORDER_MAX(k, l=k>>1))
      {
	SWAP(k, l);
	k = l;
      }
  else
    while(k > 1 && ORDER_MIN(k, l=k>>1))
      {
	SWAP(k, l);
	k = l;
      }
}

#define ORDER_MAX2(k,l) (tabPoints[k].value >= tabPoints[l].value)
#define ORDER_MIN2(k,l) (tabPoints[k].value <= tabPoints[l].value)
/* Put the first neighbor at a position so that we fix the heap */
void fix_down(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor *tabPoints = pNeighborhood->tabPoints;
  int N = pNeighborhood->iNbPoints, k = 1, l;

  if(pNeighborhood->type == MAX)
    while((l = k << 1) <= N)
      {
	if(l < N && ORDER_MAX(l+1,l)) ++l;
	if(ORDER_MAX2(k,l))
	  break;
	SWAP(k, l);
	k = l;
      }
  else
    while((l = k << 1) <= N)
      {
	if(l < N && ORDER_MIN(l+1,l)) ++l;
	if(ORDER_MIN2(k,l))
	  break;
	SWAP(k, l);
	k = l;
      }
}

/* Add the pixel (x,y), of gray-level VALUE, to the neighbor pixels */
void add_neighbor(pNeighborhood, x, y, value)
Neighborhood* pNeighborhood;
short int x, y;
float value;
{
  Neighbor* pNewNeighbor;

  /* 0) Tag the pixel as 'visited neighbor' */
  tabtabVisitedNeighbors[y][x] = iExploration; 
  /* 1) Update the fields TYPE and OTHERBOUND of PNEIGHBORHOOD */
  if(pNeighborhood->iNbPoints == 0)
    pNeighborhood->otherBound = value;
  else
    switch(pNeighborhood->type)
      {
      case MIN:
	if(value > pNeighborhood->otherBound)
	  pNeighborhood->otherBound = value;
	else if(value < pNeighborhood->tabPoints[1].value)
	  pNeighborhood->type = INVALID;
	break;
      case MAX:
	if(value < pNeighborhood->otherBound)
	  pNeighborhood->otherBound = value;
	else if(value > pNeighborhood->tabPoints[1].value)
	  pNeighborhood->type = INVALID;
	break;
      case AMBIGUOUS:
	if(value != pNeighborhood->tabPoints[1].value) {
	  pNeighborhood->type =
	    (value < pNeighborhood->tabPoints[1].value)? MAX: MIN;
	  pNeighborhood->otherBound = value;
	}
	break;
      }
  if(pNeighborhood->type == INVALID)
    return;
  /* 2) Add the point in the heap and update it */
  pNewNeighbor = &pNeighborhood->tabPoints[++pNeighborhood->iNbPoints];
  pNewNeighbor->point.x = x; /* Initialise the new neighbor point */
  pNewNeighbor->point.y = y;
  pNewNeighbor->value = value;
  fix_up(pNeighborhood); /* Update the heap of neighbors */
}

/* Remove the neighbor at the top of the heap, that is the root of the tree. */
void remove_neighbor(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor* pTop = &pNeighborhood->tabPoints[1];
  float value = pTop->value;

  if(pNeighborhood->type == INVALID)
    return;
  *pTop = pNeighborhood->tabPoints[pNeighborhood->iNbPoints--];
  if(pNeighborhood->iNbPoints == 0)
    return;
  fix_down(pNeighborhood);
  if(value != pTop->value && pTop->value == pNeighborhood->otherBound)
    pNeighborhood->type = AMBIGUOUS;
}

/* ------------------------------------------------------------------------
   --------- Allocations of structures used in the algorithm --------------
   ------------------------------------------------------------------------ */

/* Allocate image of the tags for visited pixels and the visited neighbors.
Do not be afraid about the parameters: pointers to 2-D arrays. */
void init_image_of_visited_pixels(ptabtabVisitedPixels)
int*** ptabtabVisitedPixels;
{
  int i;
  
  ARRAY_2D_ALLOC(*ptabtabVisitedPixels, iWidth, iHeight, int);
  ARRAY_2D_ALLOC(tabtabVisitedNeighbors, iWidth, iHeight, int);
}

void free_image_of_visited_pixels(tabtabVisitedPixels)
int** tabtabVisitedPixels;
{
  free(tabtabVisitedPixels[0]);  /* Actually a 1-D array */
  free(tabtabVisitedPixels);

  free(tabtabVisitedNeighbors[0]);
  free(tabtabVisitedNeighbors);
}

/* Initialize the output image */
void init_output_image(tabPixelsIn, ptabtabPixelsOutput)
float *tabPixelsIn, ***ptabtabPixelsOutput;
{
  int i;

  *ptabtabPixelsOutput = (float**) malloc(iHeight * sizeof(float*));
  if(*ptabtabPixelsOutput == 0)
    mwerror(FATAL, 1, "init_output_image --> allocation error\n");
  for(i = 0; i < iHeight; i++)
    (*ptabtabPixelsOutput)[i] = tabPixelsIn + i * iWidth;
}

void free_output_image(tabtabPixelsOutput)
float** tabtabPixelsOutput;
{
  free(tabtabPixelsOutput);
}

void init_region(iMaxArea)
int iMaxArea;
{
  tabPointsInShape = (Point_plane) malloc(iMaxArea*sizeof(struct point_plane));
  if(tabPointsInShape == NULL)
    mwerror(FATAL, 1, "init_region --> impossible to allocate the array\n");
}

void free_region()
{
  free(tabPointsInShape);
}

/* ------------------------------------------------------------------------
   --------- The core extraction algorithm --------------------------------
   ------------------------------------------------------------------------ */

/* Is pixel (x, y) a local minimum? */
char is_local_min(ou, x, y, b8Connected)
float** ou; /* A 2-D array of the image */
short int x, y;
char b8Connected;
{
  float v;
  char n = 0;

  v = ou[y][x];
  return (x==iWidth-1 || (ou[y][x+1]>v && ++n) || ou[y][x+1]==v) &&
    (x==0 || (ou[y][x-1]>v && ++n) || ou[y][x-1]==v) &&
    (y==iHeight-1 || (ou[y+1][x]>v && ++n) || ou[y+1][x]==v) &&
    (y==0 || (ou[y-1][x]>v && ++n) || ou[y-1][x]==v) &&
    (b8Connected == 0 ||
     ((x==iWidth-1 || y==0 || (ou[y-1][x+1]>v  && ++n) || ou[y-1][x+1]==v) &&
      (x==iWidth-1||y==iHeight-1||(ou[y+1][x+1]>v && ++n)|| ou[y+1][x+1]==v) &&
      (x==0 || y==iHeight-1 || (ou[y+1][x-1]>v && ++n) || ou[y+1][x-1]==v) &&
      (x==0 || y==0 || (ou[y-1][x-1]>v && ++n) || ou[y-1][x-1]==v)))	&&
    n != 0;
}

/* Is pixel (x,y) a local maximum? */
char is_local_max(ou, x, y, b8Connected)
float** ou; /* A 2-D array of the image */
short int x, y;
char b8Connected;
{
  float v;
  char n = 0;

  v = ou[y][x];
  return (x==iWidth-1 || (ou[y][x+1]<v && ++n) || ou[y][x+1]==v) &&
    (x==0 || (ou[y][x-1]<v && ++n) || ou[y][x-1]==v) &&
    (y==iHeight-1 || (ou[y+1][x]<v && ++n) || ou[y+1][x]==v) &&
    (y==0 || (ou[y-1][x]<v && ++n) || ou[y-1][x]==v) && 
    (b8Connected == 0 ||
     ((x==iWidth-1 || y==0 || (ou[y-1][x+1]<v  && ++n) || ou[y-1][x+1]==v) &&
      (x==iWidth-1||y==iHeight-1||(ou[y+1][x+1]<v && ++n)|| ou[y+1][x+1]==v) &&
      (x==0 || y==iHeight-1 || (ou[y+1][x-1]<v && ++n) || ou[y+1][x-1]==v) &&
      (x==0 || y==0 || (ou[y-1][x-1]<v && ++n) || ou[y-1][x-1]==v))) &&
    n != 0;
}

/* Set pixels and saddle points in `tabPoints' at level newGray */
void levelize(tabtabPixelsOutput, tabPoints, iNbPoints, newGray)
float** tabtabPixelsOutput;
Point_plane tabPoints;
int iNbPoints;
float newGray;
{
  int i;
  for(i = iNbPoints - 1; i >= 0; i--)
    tabtabPixelsOutput[tabPoints[i].y][tabPoints[i].x] = newGray;
}

/* Return, coded in one byte, the local configuration around the pixel (x,y) */
unsigned char configuration(tabtabVisitedPixels, x, y)
int** tabtabVisitedPixels;
short int x, y;
{
  short int iMaxX = iWidth-1, iMaxY = iHeight-1;
  unsigned char cPattern = 0;

  if(x != 0) {