flst_bilinear.c

image processing including fourier,wavelet,segmentation etc.

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

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

extern int fsaddles();


/* 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

#define NOT_A_SADDLE FLT_MAX /* Value must be coherent with fsaddles.c */
#define is_a_saddle(x) ((x) != NOT_A_SADDLE)

/* To deal with another type of image, just change these macros, and the call
to fsaddles */
#define PIXEL_T float
#define IMAGE_T Fimage
#define IMAGE_ALLOC_FUNC mw_change_fimage

/* A 'local configuration' of the pixel is the logical OR of these values,
stored in one byte. The bit is set if the neighbor is in the 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 connected components of the complement created (sometimes negative,
since cc's can be deleted), when the pixel is appended to the region. This
is the change in the number of 4-connected components of the complement.
Configurations are stored on one byte. */
static char tabPatterns[256] =
{ 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0,
  0, 1, 1, 1, 1, 2, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0,
  0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1,
  0, 1, 1, 1, 1, 2, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0,
  0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1,
  1, 2, 2, 2, 2, 3, 2, 2, 1, 2, 2, 2, 1, 2, 1, 1,
  0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 1, 1,
  0, 1, 1, 1, 1, 2, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0,
  0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 2, 1, 1, 1, 1, 0,
  1, 1, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 0,
  1, 1, 2, 1, 2, 2, 2, 1, 2, 2, 3, 2, 2, 2, 2, 1,
  1, 1, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 0,
  0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 2, 1, 1, 1, 1, 0,
  1, 1, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 0,
  0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 2, 1, 1, 1, 1, 0,
  0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, -1};

/* A pixel or a saddle point */
typedef struct {
  short int x, y; /* For saddle point, north-east corner of dual pixel */
  char bSaddle; /* Pixel (0) or saddle point (1)? */
} PixelOrSaddle;

/* 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 {
  PixelOrSaddle point; /* Neighbor pixel or saddle point */
  PIXEL_T 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? */
  PIXEL_T 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;

/* Global variables. Ugly, but avoids weighing the code */
static int iWidth, iHeight;
static int iMinArea, iMaxArea, iAreaImage, iHalfAreaImage, iPerimeterImage;
static int iExploration; /* Used to avoid reinitializing images */
static PixelOrSaddle* tabPointsInShape;
static int** tabtabVisitedNeighbors; /* Exterior boundary */
static int** tabtabVisitedNeighborSaddles; /* Adjacent saddle points */
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) (tabNeighbors[k].value > tabNeighbors[l].value)
#define ORDER_MIN(k,l) (tabNeighbors[k].value < tabNeighbors[l].value)
#define SWAP(k,l) tabNeighbors[0] = tabNeighbors[k]; \
                  tabNeighbors[k] = tabNeighbors[l]; \
                  tabNeighbors[l] = tabNeighbors[0];
/* Put the last neighbor at a position so that we fix the heap */
void fix_up(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor* tabNeighbors = 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;
  }
}

/* Put the first neighbor at a position so that we fix the heap */
void fix_down(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor* tabNeighbors = 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_MAX(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_MIN(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, bSaddle)
Neighborhood* pNeighborhood;
short int x, y;
PIXEL_T value;
char bSaddle;
{
  Neighbor* pNewNeighbor;

  /* 0) Tag the pixel as 'visited neighbor' */
  if(bSaddle)
    tabtabVisitedNeighborSaddles[y][x] = iExploration;
  else
    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 to the heap and update it */
  pNewNeighbor = &pNeighborhood->tabPoints[++pNeighborhood->iNbPoints];
  pNewNeighbor->point.x = x; /* Initialise the new neighbor point */
  pNewNeighbor->point.y = y;
  pNewNeighbor->point.bSaddle = bSaddle;
  pNewNeighbor->value = value;
  fix_up(pNeighborhood); /* Update the heap of neighbors */
}

/* Remove  neighbor at the top of the heap, i.e., the root */
void remove_neighbor(pNeighborhood)
Neighborhood* pNeighborhood;
{
  Neighbor* pTop = &pNeighborhood->tabPoints[1];
  PIXEL_T 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, ptabtabVisitedSaddles)
int ***ptabtabVisitedPixels, ***ptabtabVisitedSaddles;
{
  int i;

  ARRAY_2D_ALLOC(*ptabtabVisitedPixels, iWidth, iHeight, int);
  ARRAY_2D_ALLOC(tabtabVisitedNeighbors, iWidth, iHeight, int);
  ARRAY_2D_ALLOC(*ptabtabVisitedSaddles, iWidth-1, iHeight-1, int);
  ARRAY_2D_ALLOC(tabtabVisitedNeighborSaddles, iWidth-1, iHeight-1, int);
}

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

  free(tabtabVisitedSaddles[0]);
  free(tabtabVisitedSaddles);

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

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

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

  *ptabtabPixelsOutput = (PIXEL_T**) malloc(iHeight * sizeof(PIXEL_T*));
  if(*ptabtabPixelsOutput == NULL)
    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)
PIXEL_T** tabtabPixelsOutput;
{
  free(tabtabPixelsOutput);
}

void init_region(iMaxArea)
int iMaxArea;
{
  tabPointsInShape = (PixelOrSaddle*)malloc(iMaxArea * sizeof(PixelOrSaddle));
  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)
PIXEL_T** ou; /* A 2-D array of the image */
short int x, y;
{
  PIXEL_T 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) &&
    n != 0;
}

/* Is pixel (x,y) a local maximum? */
char is_local_max(ou, x, y)
PIXEL_T** ou; /* A 2-D array of the image */
short int x, y;
{
  PIXEL_T 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) && 
    n != 0;
}

/* Set pixels and saddle points in `tabPoints' at level newGray */
void levelize(tabtabPixelsOutput, tabtabSaddleValues,
	      tabPoints, iNbPoints, newGray)
PIXEL_T** tabtabPixelsOutput;
float** tabtabSaddleValues;
PixelOrSaddle* tabPoints;
int iNbPoints;
PIXEL_T newGray;
{
  int i;
  for(i = iNbPoints - 1; i >= 0; i--)
    if(tabPoints[i].bSaddle)
      tabtabSaddleValues[tabPoints[i].y][tabPoints[i].x] = newGray;
    else
      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, tabtabVisitedSaddles, x, y)
int **tabtabVisitedPixels, **tabtabVisitedSaddles;
short int x, y;
{
  short int iMaxX = iWidth-1, iMaxY = iHeight-1;
  unsigned char cPattern = 0;

  if(x != 0) {
    if(tabtabVisitedPixels[y][x-1] == iExploration)
      cPattern = WEST;
  } else if(iAtBorder)
    cPattern = SOUTH_WEST | WEST | NORTH_WEST;
  
  if(x != iMaxX) {
    if(tabtabVisitedPixels[y][x+1] == iExploration)
      cPattern |= EAST;
  } else if(iAtBorder)
    cPattern |= SOUTH_EAST | EAST | NORTH_EAST;
  
  if(y != 0) {
     if(tabtabVisitedPixels[y-1][x] == iExploration)
       cPattern |= NORTH;
  } else if(iAtBorder)
    cPattern |= NORTH_EAST | NORTH | NORTH_WEST;

  if(y != iMaxY) {
     if(tabtabVisitedPixels[y+1][x] == iExploration)
       cPattern |= SOUTH;
  } else if(iAtBorder)
    cPattern |= SOUTH_EAST | SOUTH | SOUTH_WEST;

#define DIAGONAL_IN_REGION(pixel, x, y, diag) ( \
  (pixel == iExploration && \
        ((cPattern&(diag)) || tabtabVisitedSaddles[y][x] == iExploration)) || \
  ((cPattern&(diag))==(diag) && tabtabVisitedSaddles[y][x] == iExploration) \
 )
  if(x != 0 && y != 0 &&
     DIAGONAL_IN_REGION(tabtabVisitedPixels[y-1][x-1], x-1, y-1, NORTH|WEST))
    cPattern |= NORTH_WEST;
  if(x != 0 && y != iMaxY &&
     DIAGONAL_IN_REGION(tabtabVisitedPixels[y+1][x-1], x-1, y, SOUTH|WEST))
    cPattern |= SOUTH_WEST;
  if(x != iMaxX && y != 0 &&
     DIAGONAL_IN_REGION(tabtabVisitedPixels[y-1][x+1], x, y-1, NORTH|EAST))