flst.c

image processing including fourier,wavelet,segmentation etc.

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

  if((y == iMaxY && iAtBorder) ||
     (y < iMaxY && tabtabVisitedPixels[y+1][x] == iExploration))
    cPattern |= SOUTH;

  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;
float 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)
Point_plane 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--)
    {
      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)
Point_plane 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;
    pShape = tabConnections[iIndex].shape;
    if(pShape != NULL) {
      level = tabConnections[iIndex].level;
      pParent = pSmallestShape;
      while(pParent->value != level) {
	assert(pParent != &pGlobalTree->the_shapes[pGlobalTree->nb_shapes-1]);
	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)
Point_plane pPoint;
float level;
Connection* tabConnections;
{
  int iIndex;
  Shape pSibling, pShape = &pGlobalTree->the_shapes[pGlobalTree->nb_shapes-1];

  iIndex = pPoint->y*iWidth + pPoint->x;
  if(tabConnections[iIndex].shape == NULL) {
    tabConnections[iIndex].shape = pShape;
    tabConnections[iIndex].level = level;
  } else {
    assert(tabConnections[iIndex].level == level);
    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)
float** 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]);
  if(x < iWidth-1  && NEIGHBOR_NOT_STORED(x+1,y))
    add_neighbor(pNeighborhood, x+1, y, ou[y][x+1]);
  if(y > 0         && NEIGHBOR_NOT_STORED(x,y-1))
    add_neighbor(pNeighborhood, x, y-1, ou[y-1][x]);
  if(y < iHeight-1 && NEIGHBOR_NOT_STORED(x,y+1))
    add_neighbor(pNeighborhood, x, y+1, ou[y+1][x]);
}

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

/* Add the points in the neighborhood of gray level currentGrayLevel to the
region tabPointsInShape and return 1 if a new shape is (maybe) detected.
This "maybe" is linked to `pIgnoreHoles', indicating if we can count the
holes. New points are added to `tabPointsInShape' from position `pCurrentArea'.
This value is changed at exit in case of success. */
char add_iso_level(tabPointsInShape, pCurrentArea,
		   currentGrayLevel, pNeighborhood, ou,
		   tabtabVisitedPixels, p8Connected, pIgnoreHoles)
Point_plane tabPointsInShape;
int* pCurrentArea;
float currentGrayLevel;
Neighborhood* pNeighborhood;
float** ou; 
int** tabtabVisitedPixels;
char* p8Connected;
char* pIgnoreHoles;
{
  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 */
    x = pNeighbor->point.x;
    y = pNeighbor->point.y;
    tabPointsInShape[iCurrentArea].x = x;
    tabPointsInShape[iCurrentArea++].y = y;
    if(! *pIgnoreHoles) {
      cPattern = configuration(tabtabVisitedPixels, x, y);
      iNbHoles += tabPatterns[*p8Connected][cPattern];
    }
    if(x == 0 || x == iWidth-1 || y == 0 || y == iHeight-1)
      ++ iAtBorder;
    tabtabVisitedPixels[y][x] = iExploration;
    /* 2) Store new neighbors */
    store_4neighbors(ou, x, y, pNeighborhood);
    if(pNeighborhood->type == MAX) {
      if(! *p8Connected)
	*pIgnoreHoles = *p8Connected = (char)1;
      store_8neighbors(ou, x, y, pNeighborhood);
    }
    remove_neighbor(pNeighborhood);
  } while(iCurrentArea <= iMaxArea &&
	  pNeighbor->value == currentGrayLevel &&
	  pNeighborhood->type != INVALID);

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

/* Extract the terminal branch containing the point (x,y) */
void find_terminal_branch(ou, tabtabVisitedPixels, x, y,
			  b8Connected, pNeighborhood, tabConnections)
float **ou;
int** tabtabVisitedPixels;
short int x, y;
char b8Connected;
Neighborhood* pNeighborhood;
Connection* tabConnections;
{
  float level;
  int iArea, iLastShapeArea;
  char bUnknownHoles;
  Shape pSmallestShape, pLastShape;

 restart_branch:
  ++ iExploration;
  iArea = iAtBorder = 0;
  bUnknownHoles = 0;
  pSmallestShape = pLastShape = NULL;
  level = ou[y][x];
  reinit_neighborhood(pNeighborhood, b8Connected ? MAX: MIN);
  add_neighbor(pNeighborhood, x, y, level);
  while(add_iso_level(tabPointsInShape, &iArea,
		      level, pNeighborhood, ou, tabtabVisitedPixels,
		      &b8Connected, &bUnknownHoles) != 0) {
    if(bUnknownHoles) {
      assert(iArea != 0);
      if(pLastShape != NULL) {
	iArea = pLastShape->area;
	connect(tabPointsInShape, iArea, tabConnections, pSmallestShape);
	new_connection(&tabPointsInShape[iArea], level, tabConnections);
      }
      levelize(ou, tabPointsInShape, iArea, level);
      goto restart_branch;
    }
    if(iMinArea <= iArea) { /* Store new shape? */
      iLastShapeArea = (pLastShape == NULL) ? 0 : pLastShape->area;
      pLastShape = new_shape(iArea, level, !b8Connected, pLastShape);
      if(pSmallestShape == NULL) pSmallestShape = pLastShape;
      update_smallest_shapes(tabPointsInShape+iLastShapeArea,
			     iArea-iLastShapeArea);
    }
    if(iAtBorder && iArea == iHalfAreaImage)
      break;
    bUnknownHoles = (char)(b8Connected && pNeighborhood->type == AMBIGUOUS);
    if(bUnknownHoles) b8Connected = 0;
    level = pNeighborhood->tabPoints[1].value;
  }
  if(pLastShape != NULL) {
    connect(tabPointsInShape, iArea, tabConnections, pSmallestShape);
    if(iAtBorder && iArea == iHalfAreaImage)
      insert_children(pGlobalTree->the_shapes, pLastShape);
    else if(iArea != 0)
      new_connection(&tabPointsInShape[iArea], level, tabConnections);
  }
  levelize(ou, tabPointsInShape, iArea, level);
}

/* Scan the image, calling a procedure to extract terminal branch at each
   (not yet visited) local extremum */
void scan(tabtabPixelsOutput, tabtabVisitedPixels,pNeighborhood,tabConnections)
float **tabtabPixelsOutput;
int** tabtabVisitedPixels;
Neighborhood* pNeighborhood;
Connection* tabConnections;
{
  short int i, j;
  char b8Connected = 0;
  int iExplorationInit;

  iExplorationInit = iExploration;
  for(i = 0; i < iHeight; i++)
    for(j = 0; j < iWidth; j++)
      if(tabtabVisitedPixels[i][j] <= iExplorationInit &&
	 (is_local_min(tabtabPixelsOutput, j, i, (char)0) ||
	  (is_local_max(tabtabPixelsOutput,j,i,(char)1) &&(b8Connected=1)==1)))
	{
	  find_terminal_branch(tabtabPixelsOutput, tabtabVisitedPixels, j, i,
			       b8Connected, pNeighborhood, tabConnections);
	  b8Connected = 0;
	}
}

/* ------------------------------------------------------------------------
   --------- 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 put in the tree. pMinArea==NULL means 1.
Output: *pTree is filled */
void flst(pMinArea, pImageInput, pTree)
int* pMinArea;
Fimage pImageInput;
Shapes pTree;
{
  float **tabtabPixelsOutput; /* Array accessing pixels of output image */
  Neighborhood neighborhood; /* The neighborhood of the current region */
  int** tabtabVisitedPixels;  /* Image of last visit for each pixel */
  Connection* tabConnections;
  int i;

  pGlobalTree = pTree;
  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]);
  pGlobalTree->interpolation = 0;
  tabConnections = (Connection*) malloc(iAreaImage * sizeof(Connection));
  if(tabConnections == NULL)
    mwerror(FATAL, 1, "Image of largest shapes: allocation error\n");
  for(i = iAreaImage-1; i >= 0; i--)
    tabConnections[i].shape = NULL;

  init_output_image(pImageInput->gray, &tabtabPixelsOutput);
  init_image_of_visited_pixels(&tabtabVisitedPixels);
  init_neighborhood(&neighborhood, iAreaImage);
  init_region(iAreaImage);

  iExploration = 0;
  /* Loop with increasing iMaxArea: speed optimization. Result correct
  with only one call to `scan' and iMaxArea = iAreaImage */
  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,&neighborhood,tabConnections);
  } while(iMaxArea+1 < iAreaImage);

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

  free(tabConnections);
  free_image_of_visited_pixels(tabtabVisitedPixels);
  free_region();
  free_neighborhood(&neighborhood);
  free_output_image(tabtabPixelsOutput);
}