flst_bilinear.c

image processing including fourier,wavelet,segmentation etc.

    cPattern |= NORTH_EAST;
  if(x != iMaxX && y != iMaxY &&
     DIAGONAL_IN_REGION(tabtabVisitedPixels[y+1][x+1], x, y, SOUTH|EAST))
    cPattern |= SOUTH_EAST;

  return cPattern;
}

/* Insert a new shape and its siblings in the tree, with parent pParent */
void insert_children(pParent, pNewChildToInsert)
Shape pParent, pNewChildToInsert;
{
  Shape pSibling = pNewChildToInsert;
  while(pSibling->next_sibling != NULL) {
    pSibling->parent = pParent;
    pSibling = pSibling->next_sibling;
  }
  pSibling->parent = pParent;
  pSibling->next_sibling = pParent->child;
  pParent->child = pNewChildToInsert;
}

Shape new_shape(iCurrentArea, currentGrayLevel, bOfInferiorType, pChild)
int iCurrentArea;
PIXEL_T currentGrayLevel;
char bOfInferiorType;
Shape pChild; /* Supposed to have no sibling. Can be NULL */
{
  Shape pNewShape = &pGlobalTree->the_shapes[pGlobalTree->nb_shapes++];

  pNewShape->inferior_type = bOfInferiorType;
  pNewShape->value = currentGrayLevel;
  pNewShape->open = (char)(iAtBorder != 0);
  pNewShape->area = iCurrentArea;
  pNewShape->removed = 0;
  pNewShape->pixels = NULL;
  pNewShape->boundary = pNewShape->data = NULL;
  pNewShape->data_size = 0;
  /* Make links */
  pNewShape->parent = NULL;
  pNewShape->next_sibling = NULL;
  pNewShape->child = pChild;
  if(pChild != NULL)
    pChild->parent = pNewShape;
  return pNewShape;
}

/* Knowing that the last extracted shape contains the points, update,
for each one, the smallest shape containing it */
void update_smallest_shapes(tabPoints, iNbPoints)
PixelOrSaddle* tabPoints;
int iNbPoints;
{
  int i, iIndex;
  Shape pNewShape, pRoot = &pGlobalTree->the_shapes[0];

  pNewShape = &pGlobalTree->the_shapes[pGlobalTree->nb_shapes-1];
  for(i = iNbPoints - 1; i >= 0; i--)
    if(! tabPoints[i].bSaddle) {
      iIndex = tabPoints[i].y * iWidth + tabPoints[i].x;
      if(pGlobalTree->smallest_shape[iIndex] == pRoot)
	pGlobalTree->smallest_shape[iIndex] = pNewShape;
    }
}

/* Find children of the last constructed monotone section, which is composed
of the interval between pSmallestShape and the last extracted shape. That is,
find shapes in other monotone sections whose parent is inside this interval */
void connect(tabPoints, iNbPoints, tabConnections, pSmallestShape)
PixelOrSaddle* tabPoints;
int iNbPoints;
Connection* tabConnections;
Shape pSmallestShape;
{
  int i, iIndex;
  Shape pShape, pParent;
  float level;

  for(i = iNbPoints-1; i >= 0; i--) {
    iIndex = (tabPoints[i].y * iWidth + tabPoints[i].x) << 1;
    if(tabPoints[i].bSaddle) ++ iIndex;
    pShape = tabConnections[iIndex].shape;
    if(pShape != NULL) {
      level = tabConnections[iIndex].level;
      pParent = pSmallestShape;
      while(pParent->value != level)
	pParent = pParent->parent;
      insert_children(pParent, pShape);
      tabConnections[iIndex].shape = NULL;
    }
  }
}

/* Make a new connection structure at the given point */
void new_connection(pPoint, level, tabConnections)
PixelOrSaddle* pPoint;
float level;
Connection* tabConnections;
{
  int iIndex;
  Shape pSibling, pShape = &pGlobalTree->the_shapes[pGlobalTree->nb_shapes-1];

  iIndex = (pPoint->y*iWidth + pPoint->x) << 1;
  if(pPoint->bSaddle) ++ iIndex;
  if(tabConnections[iIndex].shape == NULL) {
    tabConnections[iIndex].shape = pShape;
    tabConnections[iIndex].level = level;
  } else {
    pSibling = tabConnections[iIndex].shape;
    while(pSibling->next_sibling != NULL)
      pSibling = pSibling->next_sibling;
    pSibling->next_sibling = pShape;
  }
}

/* Is the neighbor pixel already stored for this exploration? */
#define NEIGHBOR_NOT_STORED(x,y) (tabtabVisitedNeighbors[y][x] < iExploration)

/* Store the 4-neighbors of pixel (x,y) */
void store_4neighbors(ou, x, y, pNeighborhood)
PIXEL_T** ou;
short int x, y;
Neighborhood* pNeighborhood;
{
  if(x > 0         && NEIGHBOR_NOT_STORED(x-1,y))
    add_neighbor(pNeighborhood, x-1, y, ou[y][x-1], (char)0);
  if(x < iWidth-1  && NEIGHBOR_NOT_STORED(x+1,y))
    add_neighbor(pNeighborhood, x+1, y, ou[y][x+1], (char)0);
  if(y > 0         && NEIGHBOR_NOT_STORED(x,y-1))
    add_neighbor(pNeighborhood, x, y-1, ou[y-1][x], (char)0);
  if(y < iHeight-1 && NEIGHBOR_NOT_STORED(x,y+1))
    add_neighbor(pNeighborhood, x, y+1, ou[y+1][x], (char)0);
}

/* Store the neighbors of the saddle point (x,y) */
void store_neighbors_to_saddle(ou, x, y, pNeighborhood)
PIXEL_T** ou;
short int x, y;
Neighborhood* pNeighborhood;
{
  if(NEIGHBOR_NOT_STORED(x,y))
    add_neighbor(pNeighborhood, x, y, ou[y][x], (char)0);
  if(NEIGHBOR_NOT_STORED(x+1,y))
    add_neighbor(pNeighborhood, x+1, y, ou[y][x+1], (char)0);
  if(NEIGHBOR_NOT_STORED(x,y+1))
    add_neighbor(pNeighborhood, x, y+1, ou[y+1][x], (char)0);
  if(NEIGHBOR_NOT_STORED(x+1,y+1))
    add_neighbor(pNeighborhood, x+1, y+1, ou[y+1][x+1], (char)0);  
}

/* Is the neighbor saddle point already stored? */
#define SADDLE_NOT_STORED(x,y) \
  (is_a_saddle(ou[y][x]) && tabtabVisitedNeighborSaddles[y][x] < iExploration)

/* Store the saddle points being neighbors of pixel (x,y) */
void store_saddle_neighbors(ou, x, y, pNeighborhood)
float** ou;
short int x, y;
Neighborhood* pNeighborhood;
{
  if(x > 0) {
    if(y > 0         && SADDLE_NOT_STORED(x-1,y-1))
      add_neighbor(pNeighborhood, x-1, y-1, ou[y-1][x-1], (char)1);
    if(y < iHeight-1 && SADDLE_NOT_STORED(x-1,y))
      add_neighbor(pNeighborhood, x-1, y, ou[y][x-1], (char)1);
  }
  if(x < iWidth-1) {
    if(y > 0         && SADDLE_NOT_STORED(x,y-1))
      add_neighbor(pNeighborhood, x, y-1, ou[y-1][x], (char)1);
    if(y < iHeight-1 && SADDLE_NOT_STORED(x,y))
      add_neighbor(pNeighborhood, x, y, ou[y][x], (char)1);
  }
}

/* Does the addition of saddle point (x,y) change the Euler number? */
char saddle_change(tabtabVisitedPixels, x, y)
int** tabtabVisitedPixels;
short int x, y;
{
  char cDiag1, cDiag2, cDiag3;

  cDiag1 = (tabtabVisitedPixels[y][x] == iExploration);
  if(cDiag1)
    cDiag2 = (tabtabVisitedPixels[y+1][x+1] == iExploration);
  else
    cDiag2 = (tabtabVisitedPixels[y+1][x+1] != iExploration);

  if(cDiag2) { /* Pixels of 1st diag inside (cDiag1) or outside (!cDiag1)? */
    cDiag2 = (tabtabVisitedPixels[y+1][x] == iExploration);
    if(cDiag2)
      cDiag3 = (tabtabVisitedPixels[y][x+1] == iExploration);
    else
      cDiag3 = (tabtabVisitedPixels[y][x+1] != iExploration);
    /* cDiag3: Pixels of 2nd diag inside (cDiag2) or outside (!cDiag2)? */
    if(cDiag3 && cDiag1 != cDiag2)
      return (char)1;
  }
  return 0;
}

/* Add the points in the neighborhood of gray level `currentGrayLevel' to the
region `tabPointsInShape' and return 1 if a new shape is detected. New points
are added from position `pCurrentArea'. This value is changed at exit in case
of success. */
char add_iso_level(tabPointsInShape, pCurrentArea, pNbPixels,
		   currentGrayLevel, pNeighborhood, ou, tabtabVisitedPixels,
		   tabtabVisitedSaddles, tabtabSaddleValues)
PixelOrSaddle* tabPointsInShape;
int *pCurrentArea, *pNbPixels;
PIXEL_T currentGrayLevel, **ou;
Neighborhood* pNeighborhood;
int **tabtabVisitedPixels, **tabtabVisitedSaddles;
float** tabtabSaddleValues;
{
  short int x, y;
  Neighbor* pNeighbor;
  int iCurrentArea, iNbHoles;
  unsigned char cPattern;

  iNbHoles = 0;
  iCurrentArea = *pCurrentArea;
  pNeighbor = &pNeighborhood->tabPoints[1];
  do { /* 1) Neighbor is added to the region */
    tabPointsInShape[iCurrentArea].x = x = pNeighbor->point.x;
    tabPointsInShape[iCurrentArea].y = y = pNeighbor->point.y;
    tabPointsInShape[iCurrentArea++].bSaddle = pNeighbor->point.bSaddle;
    if(! pNeighbor->point.bSaddle) {
      ++ *pNbPixels;
      cPattern = configuration(tabtabVisitedPixels, tabtabVisitedSaddles, x,y);
      iNbHoles += tabPatterns[cPattern];
      if(x == 0 || x == iWidth-1 || y == 0 || y == iHeight-1)
	++ iAtBorder;
      tabtabVisitedPixels[y][x] = iExploration;
    } else {
      if(saddle_change(tabtabVisitedPixels, x, y))
	++ iNbHoles;
      tabtabVisitedSaddles[y][x] = iExploration;
    }
    /* 2) Store new neighbors */
    if(! pNeighbor->point.bSaddle) {
      store_4neighbors(ou, x, y, pNeighborhood);
      store_saddle_neighbors(tabtabSaddleValues, x, y, pNeighborhood);
    } else
      store_neighbors_to_saddle(ou, x, y, pNeighborhood);
    remove_neighbor(pNeighborhood);
  } while(*pNbPixels <= iMaxArea &&
	  pNeighbor->value == currentGrayLevel &&
	  pNeighborhood->type != INVALID);

  if(*pNbPixels <= iMaxArea &&
     iAtBorder != iPerimeterImage &&
     (! iAtBorder || *pNbPixels <= iHalfAreaImage) &&
     pNeighborhood->type != INVALID &&
     iNbHoles == 0) {
    *pCurrentArea = iCurrentArea;
    return (char)1;
  }
  return 0;
}

/* Extract the terminal branch containing the point (x,y) */
void find_terminal_branch(ou, tabtabVisitedPixels, tabtabVisitedSaddles,
			  tabtabSaddleValues, x, y,
			  pNeighborhood, type, tabConnections)
PIXEL_T **ou;
int **tabtabVisitedPixels, **tabtabVisitedSaddles;
float** tabtabSaddleValues;
short int x, y;
Neighborhood* pNeighborhood;
enum TypeOfTree type;
Connection* tabConnections;
{
  PIXEL_T level;
  int iArea=0, iLastArea=0, iNbPixels=0;
  Neighbor* pNext = &pNeighborhood->tabPoints[1];
  Shape pSmallestShape=NULL, pLastShape=NULL;

  ++ iExploration;
  iAtBorder = 0; /* Shape does not meet the image border yet */
  level = ou[y][x];
  reinit_neighborhood(pNeighborhood, type);
  add_neighbor(pNeighborhood, x, y, level, (char)0);
  while(add_iso_level(tabPointsInShape, &iArea, &iNbPixels,
		      level, pNeighborhood, ou, tabtabVisitedPixels,
		      tabtabVisitedSaddles, tabtabSaddleValues) != 0) {
    if((type == MAX && level < pNext->value) ||
       (type == MIN && level > pNext->value)) { /* Type change */
      type = (level < pNext->value) ? MIN : MAX;
      if(pLastShape != NULL)
	connect(tabPointsInShape, iLastArea, tabConnections, pSmallestShape);
      pSmallestShape = NULL;
    }
    if(iMinArea <= iNbPixels) { /* Store new shape? */
      pLastShape = new_shape(iNbPixels, level, (char)(type==MIN), pLastShape);
      if(pSmallestShape == NULL) pSmallestShape = pLastShape;
      update_smallest_shapes(tabPointsInShape+iLastArea, iArea-iLastArea);
      iLastArea = iArea;
    }
    if(iAtBorder && iNbPixels == iHalfAreaImage)
      break;
    level = pNext->value;
  }
  if(pLastShape != NULL) {
    connect(tabPointsInShape, iArea, tabConnections, pSmallestShape);
    if(iAtBorder && iNbPixels == iHalfAreaImage)
      insert_children(pGlobalTree->the_shapes, pLastShape);
    else if(iArea != 0)
      new_connection(&tabPointsInShape[iArea], level, tabConnections);
  }
  levelize(ou, tabtabSaddleValues, tabPointsInShape, iArea, level);
}

/* Scan the image, calling a procedure to extract terminal branch at each
(not yet visited) local extremum */
void scan(tabtabPixelsOutput, tabtabVisitedPixels, tabtabVisitedSaddles,
	  tabtabSaddleValues, pNeighborhood, tabConnections)
PIXEL_T **tabtabPixelsOutput;
int **tabtabVisitedPixels, **tabtabVisitedSaddles;
float** tabtabSaddleValues;
Neighborhood* pNeighborhood;
Connection* tabConnections;
{
  short int i, j;
  enum TypeOfTree type;
  int iExplorationInit;

  iExplorationInit = iExploration;
  type = MIN;
  for(i = 0; i < iHeight; i++)
    for(j = 0; j < iWidth; j++)
      if(tabtabVisitedPixels[i][j] <= iExplorationInit &&
	 (is_local_min(tabtabPixelsOutput, j, i) ||
	  (is_local_max(tabtabPixelsOutput, j, i) && (type=MAX)==MAX))) {
	find_terminal_branch(tabtabPixelsOutput,
			     tabtabVisitedPixels, tabtabVisitedSaddles,
			     tabtabSaddleValues, j, i,
			     pNeighborhood, type, tabConnections);
	type = MIN;
      }
}

/* ------------------------------------------------------------------------
   --------- The main function --------------------------------------------
   ------------------------------------------------------------------------ */

/* The "Fast Level Set Transform" gives the tree of interiors of level lines
(named 'shapes') representing the image.
Only shapes of area >= *pMinArea are in the tree. pMinArea==NULL means 1.
Output: *pTree is filled (pTree must point to an allocated tree). */
void flst_bilinear(pMinArea, pImageInput, pTree)
int* pMinArea;
IMAGE_T pImageInput;
Shapes pTree;
{
  PIXEL_T** tabtabPixelsOutput; /* Array accessing pixels of output image */
  struct fimage imageSaddles;
  float** tabtabSaddleValues;
  Neighborhood neighborhood; /* Neighborhood of current region */
  int** tabtabVisitedPixels; /* Image of last visit for each pixel */
  int** tabtabVisitedSaddles;
  int iNbSaddles;
  Connection* tabConnections;
  int i;

  iWidth = pImageInput->ncol;
  iHeight = pImageInput->nrow;
  iAreaImage = iWidth * iHeight; iHalfAreaImage = (iAreaImage+1) / 2;
  if(iWidth == 1) iPerimeterImage = iHeight;
  else if(iHeight == 1) iPerimeterImage = iWidth;
  else iPerimeterImage = (iWidth+iHeight-2)*2;
  iMinArea = (pMinArea != NULL) ? *pMinArea : 1;
  if(iMinArea > iAreaImage) mwerror(USAGE, 1, "min area > image");

  pGlobalTree = mw_change_shapes(pTree, iHeight, iWidth, pImageInput->gray[0]);
  if(pGlobalTree == NULL)
    mwerror(FATAL, 1, "shapes allocation error");
  pTree->interpolation = 1;
  /* tabConnection[2*k] <-> pixel(k), tabConnection[2*k+1] <-> saddle(k) */
  tabConnections = (Connection*) malloc(2*iAreaImage * sizeof(Connection));
  if(tabConnections == NULL)
    mwerror(FATAL, 1, "allocation of image of connections error\n");
  for(i = 2*iAreaImage-1; i >= 0; i--)
    tabConnections[i].shape = NULL;

  init_output_image(pImageInput->gray, &tabtabPixelsOutput);
  imageSaddles.nrow = imageSaddles.ncol = 0; imageSaddles.gray = NULL;
  imageSaddles.allocsize = 0;
  iNbSaddles = fsaddles(pImageInput, &imageSaddles);
  if((tabtabSaddleValues = mw_newtab_gray_fimage(&imageSaddles)) == NULL)
    mwerror(FATAL, 1, "Lines of saddle values: allocation error");
  init_image_of_visited_pixels(&tabtabVisitedPixels, &tabtabVisitedSaddles);
  init_neighborhood(&neighborhood, iAreaImage + iNbSaddles);
  init_region(iAreaImage + iNbSaddles);
  /* Reallocation, the kernel allocating for a 0-order interpolation */
  pTree->the_shapes = (Shape)
    realloc(pTree->the_shapes, (iNbSaddles+iAreaImage+1)*sizeof(struct shape));
  for(i = iAreaImage-1; i >= 0; i--)
    pGlobalTree->smallest_shape[i] = pGlobalTree->the_shapes;

  iExploration = 0;
  /* Loop with increasing iMaxArea: speed optimization. Result correct
  with only one call to `scan' and iMaxArea = iAreaImage-1 */
  iMaxArea = 0;
  do {
    if(iMaxArea == 0)
      iMaxArea = INIT_MAX_AREA;
    else
      iMaxArea <<= STEP_MAX_AREA;
    if(iMaxArea == 0 || iMaxArea >= iHalfAreaImage) /* iMaxArea==0: overflow */
      iMaxArea = iAreaImage-1;
    scan(tabtabPixelsOutput, tabtabVisitedPixels, tabtabVisitedSaddles,
	 tabtabSaddleValues, &neighborhood, tabConnections);
  } while(iMaxArea+1 < iAreaImage);

  /* Make connections with root */
  pTree->the_shapes[0].value = tabtabPixelsOutput[0][0];
  for(i = 2*iAreaImage-1; i >= 0; i--)
    if(tabConnections[i].shape != NULL) {
      assert(tabConnections[i].level == pTree->the_shapes[0].value);
      insert_children(pGlobalTree->the_shapes, tabConnections[i].shape);
    }
  free(tabConnections);

  free_image_of_visited_pixels(tabtabVisitedPixels, tabtabVisitedSaddles);
  free_region();
  free_neighborhood(&neighborhood);
  free_output_image(tabtabPixelsOutput);
  free(tabtabSaddleValues[0]); free(tabtabSaddleValues);
}