LabeledSkeletonImage.C

Go to the documentation of this file.

/*---------------------------------------------------------------------+
 | Library QgarLib, graphics analysis and recognition                  |
 | Copyright (C) 2002  Qgar Project, LORIA                             |
 |                                                                     |
 | This library is free software; you can redistribute it and/or       |
 | modify it under the terms of the GNU Lesser General Public          |
 | License version 2.1, as published by the Free Software Foundation.  |
 |                                                                     |
 | This library is distributed in the hope that it will be useful,     |
 | but WITHOUT ANY WARRANTY; without even the implied warranty of      |
 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                |
 | See the GNU Lesser General Public License for more details.         |
 |                                                                     |
 | The GNU Lesser General Public License is included in the file       |
 | LICENSE.LGPL, in the root directory of the Qgar packaging. See      |
 | http://www.gnu.org/licenses/lgpl.html for the terms of the licence. |
 | To receive a paper copy, write to the Free Software Foundation,     |
 | Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.       |
 |                                                                     |
 | Contact Project Qgar for any information:                           |
 |   LORIA - équipe Qgar                                               |
 |   B.P. 239, 54506 Vandoeuvre-lès-Nancy Cedex, France                |
 |   email: qgar-contact@loria.fr                                      |
 |   http://www.qgar.org/                                              |
 *---------------------------------------------------------------------*/


/**
 * @file   LabeledSkeletonImage.C
 * @brief  Implementation of class qgar::LabeledSkeletonImage.
 *
 *         See file LabeledSkeletonImage.H for the interface.
 *
 * @author <a href="mailto:qgar-develop@loria.fr?subject=Qgar fwd Karl Tombre">Adlane Habed and Karl Tombre</a>
 * @date   July 3, 2001  16:57
 * @since  Qgar 1.0
 */



// STD
#include <vector>
// QGAR
#include <qgarlib/Dist34BlackCCImage.H>
#include <qgarlib/GenImage.H>
#include <qgarlib/LabeledSkeletonImage.H>
#include <qgarlib/primitives.H>



using namespace std;



namespace qgar
{

// -------------------------------------------------------------------
// C O N S T R U C T O R
// -------------------------------------------------------------------


LabeledSkeletonImage::LabeledSkeletonImage(const Dist34BlackCCImage& img,
                                           int maxPruning)

  : GreyLevelImage(img)

{
  _pixMap2 = new GreyLevelImage::value_type[img.width() * img.height()];
  int pixsize = _width * _height; 
  GreyLevelImage::value_type* newp = _pixMap2; 
  int i = 0;
  int j = 0;    
  
  // Initialize _pixMap2
  for (i = 0 ; i < pixsize ; ++i, ++newp)
    {
      *newp = 255;
    }

  // Change labels from 3 to 1 and from 6 to 5
  // (see referenced article)
  GreyLevelImage::value_type* oldp = _pPixMap;
  for(i = 0 ; i < pixsize ; ++i, ++oldp)
    {
      if (*oldp == 3)
        {
          *oldp = 1;
        }
      if (*oldp == 6)
        {
          *oldp = 5;
        }
    }

  // Parallelwise detectable pixels
  parallelSkelet();
  
  // Change labels back from (1,5) to (3,6)
  oldp=_pPixMap;
  for(i = 0 ; i < pixsize ; ++i, ++oldp)
    {
      if (*oldp == 1)
        {
          *oldp = 3;
        }
      if (*oldp == 5)
        {
          *oldp = 6;
        }
    }

  // Sequentially detectable skeletal pixels
  sequentSkelet();

  // Now that they are tracked, mark also these pixels
  // as belonging to skeleton 
  newp = _pixMap2;
  for(i = 0 ; i < pixsize ; ++i, ++newp)
    {
      if (*newp == 1)
        {
          *newp = 0;
        }
    }

  // Hole filling and final thinning
  finalThinning();

  // Pruning if requested
  if (maxPruning > 0)
    {
      // End points detection
      // Create an array containing all the end points
      vector<Point> myEndPoints;
      GreyLevelImage::value_type n[10];

      for (i = 1 ; i < (_height - 1) ; ++i)
        {
          oldp = _pPixMap + (i * _width) + 1;
          newp = _pixMap2 + (i * _width) + 1;

          for (j = 1 ; j < (_width - 1) ; ++j, ++newp)
            {
              if(*newp == 0) // marked point
                {
                  loadBinaryNeighbors(newp, n);
                  int x4=crossingCnt(n);
                  if(x4 == 1)
                    {
                      Point endp(j, i);
                      myEndPoints.push_back(endp);
                    }
                }
            } // END for j
        } // END for i

      pruning(myEndPoints, maxPruning);
    }

  // Remove jaggedness
  removeJaggedness(); 
  
  // And now that we have marked everything, create the labeled skeleton
  // by setting to zero all non-skeleton pixels in original distance image
  oldp = _pPixMap;
  newp = _pixMap2;

  for (i = 0 ; i < pixsize ; ++i, ++oldp, ++newp)
    {
      if ((*oldp != 0) && (*newp != 0))
        {
          *oldp = 0;
        }
    }

  // Delete local images
  delete newp;
}


// -------------------------------------------------------------------
// D E S T R U C T O R
// -------------------------------------------------------------------


// VIRTUAL DESTRUCTOR

LabeledSkeletonImage::~LabeledSkeletonImage()
{
  if (_pixMap2 != 0)
    {
      delete [] _pixMap2;
    }
}


// -------------------------------------------------------------------
// A U X I L I A R I E S
// -------------------------------------------------------------------


// FIND PARALLELWISE DETECTABLE SKELETON PIXELS
// ============================================

void
LabeledSkeletonImage::parallelSkelet()
{
  GreyLevelImage::value_type lab[10];
  GreyLevelImage::value_type neighbors[10];
  GreyLevelImage::value_type suitbin[10];

  GreyLevelImage::value_type* newp;
  GreyLevelImage::value_type* oldp;
  
  for (int i = 1 ; i < (_height - 1) ; ++i)
    {
      oldp = _pPixMap + (i * _width) + 1;
      newp = _pixMap2 + (i * _width) + 1;
      
      for(int j = 1 ; j < (_width - 1) ; ++j, ++newp, ++oldp)
        {
          if (*oldp != 0) // Not background
            {  
              // Check if it is a maximal center

              loadLabeledNeighbors(oldp,lab);
            
              if ((lab[2] < (*oldp + 4)) && (lab[3] < (*oldp + 3)) &&
                  (lab[4] < (*oldp + 4)) && (lab[1] < (*oldp + 3)) &&
                  (lab[5] < (*oldp + 3)) && (lab[6] < (*oldp + 4)) &&
                  (lab[7] < (*oldp + 3)) && (lab[8] < (*oldp + 4)))
                {
                  // any maximal center having a triple of consecutive neighbors
                  // odd/even/odd which are all labeled more than it is not
                  // marked
                  if (((lab[1] > *oldp) && (lab[2] > *oldp) && (lab[3] > *oldp)) ||
                      ((lab[3] > *oldp) && (lab[4] > *oldp) && (lab[5] > *oldp)) ||
                      ((lab[5] > *oldp) && (lab[6] > *oldp) && (lab[7] > *oldp)) || 
                      ((lab[7] > *oldp) && (lab[8] > *oldp) && (lab[1] > *oldp)))
                    {
                      *newp = 255;
                    }
                  else
                    {
                      *newp = 0; 
                    }
                }
              else 
                { // Check if we have a saddle pixel
                  loadLabeledNeighbors(oldp,neighbors);
                  loadSuitCBinary(suitbin,neighbors);
                  int c8 = connexityCnt(suitbin);

                  // loadLabeledNeighbors(oldp,neighbors);
                  loadSuitXBinary(suitbin,neighbors);
                  int x4 = crossingCnt(suitbin);

                  if (*oldp != 1)
                    {
                      if ((c8 > 1) || (x4 > 1))
                        {
                          *newp = 0;
                        }
                      else 
                        {
                          *newp = 255;
                        }
                    }
                  else
                    {
                      if ((c8 > 1) || (x4 > 1) || 
                          (oddEvenOddLab3(lab,1)==3) ||
                          (oddEvenOddLab3(lab,3)==3) ||
                          (oddEvenOddLab3(lab,5)==3) ||
                          (oddEvenOddLab3(lab,7)==3))
                        {
                          *newp = 0;
                        }
                      else
                        {
                          *newp = 255;
                        }
                    }
                }
            }
          else
            {
              *newp = 255;
            }

        } // END for j

    } // END for i
}


// FIND SEQUENTIALLY DETECTABLE SKELETON PIXELS
// ============================================


void
LabeledSkeletonImage::sequentSkelet()
{
  GreyLevelImage::value_type lab[10];
  float grad[9];
  
  for (int i = 1 ; i < (_height - 1) ; ++i)
    {
      GreyLevelImage::value_type* oldp = _pPixMap + (i * _width) + 1;
      GreyLevelImage::value_type* newp = _pixMap2 + (i * _width) + 1;

      for (int j = 1 ; j < (_width - 1) ; ++j, ++newp, ++oldp)
        {
          if(*newp == 0)  // if marked as a maximal center by previous step
            {
              // Compute gradient in distance image
              loadLabeledNeighbors(oldp, lab);
              gradient(lab, grad);
            
              int indmax = maxGradientPos1(oldp, grad);
              int indmax2 = maxGradientPos2(grad, oldp, indmax);
            
              int encore = 1;
              while (encore)
                { 
                  GreyLevelImage::value_type* savoldp = oldp; 
                  GreyLevelImage::value_type* savnewp = newp;
                  while (indmax != 0)
                    { 
                      // Go in the direction of maximum gradient
                      switch(indmax)
                        {
                        case 1:
                          --oldp;
                          --newp;
                          *newp = 1; // mark
                          break;
                        case 2:
                          oldp -= _width + 1;
                          newp -= _width + 1;
                          *newp = 1;
                          break;
                        case 3:
                          oldp -= _width;
                          newp -= _width;
                          *newp = 1;
                          break;
                        case 4:
                          oldp -= _width - 1;
                          newp -= _width - 1;
                          *newp = 1;
                          break;
                        case 5:
                          ++oldp;
                          ++newp;
                          *newp = 1;
                          break;
                        case 6:
                          oldp += _width + 1;
                          newp += _width + 1;
                          *newp = 1;
                          break;
                        case 7:
                          oldp += _width;
                          newp += _width;
                          *newp = 1;
                          break;
                        case 8:
                          oldp += _width - 1;
                          newp += _width - 1;
                          *newp = 1;
                          break;
                        default:
                          break; 
                        } // END switch
                          
                      // Compute gradient at new position
                      loadLabeledNeighbors(oldp, lab);
                      gradient(lab, grad); 
                   
                      indmax = maxGradientPos1(oldp, grad);
                    } // END while indmax
                  
                  // Restore
                  oldp = savoldp;  
                  newp = savnewp;

                  if (indmax2 != 0) // there was another starting point
                    {
                      encore = 1;
                      indmax = indmax2;
                      indmax2 = 0;
                    }
                  else
                    {
                      encore = 0;
                    }
                } // END while encore
            }
        } // END for j
      
    } // END for i
}


// HOLE FILLING AND FINAL THINNING
// ===============================


void
LabeledSkeletonImage::finalThinning()
{
  GreyLevelImage::value_type lab[10];
    
  for (int i = 1 ; i < (_height - 1) ; ++i)
    {
      GreyLevelImage::value_type* newp = _pixMap2 + (i * _width) + 1;

      for (int j = 1 ; j < (_width - 1) ; ++j, ++newp)
        {
          if (*newp == 0) // marked: reduce to unit width
            {
              loadLabeledNeighbors(newp, lab);
              unitWidth(newp, lab);
            }
          else // not marked: detect holes
            {
              if ((*(newp - 1) == 0) && (*(newp - _width) == 0) &&
                  (*(newp + 1) == 0) && (*(newp + _width) == 0))
                {  // if all odd-neighbors are marked
                  // Fill the hole
                  *newp = 0;
                  
                  // check again already inspected odd-neighbors to see
                  // if we can remove marker
                  GreyLevelImage::value_type* savp = newp;
                  
                  --newp;
                  loadLabeledNeighbors(newp, lab);
                  unitWidth(newp, lab);
                  
                  newp = savp - _width;
                  loadLabeledNeighbors(newp, lab);
                  unitWidth(newp, lab);
                  
                  // Restore pointer to continue inspection
                  newp = savp;
                }
            }
        } // END for j
    } // END for i
}


// PRUNING
// =======

void
LabeledSkeletonImage::pruning(vector< GenPoint<int> >& myEndPoints,
                              int maxPruning)
{
  GreyLevelImage::value_type n[10];
  vector<Point>::iterator myIter = myEndPoints.begin();

  // Create a vector of pointers
  vector<GreyLevelImage::value_type*> pointsToBePruned;

  while (myIter != myEndPoints.end())
    {
      // Get the point
      Point endP = *myIter;
      // Put pointer on right position
      GreyLevelImage::value_type* startOldP =  _pPixMap + (endP.y() * _width) + endP.x();
      GreyLevelImage::value_type* startNewP =  _pixMap2 + (endP.y() * _width) + endP.x();

      GreyLevelImage::value_type* currOldP = startOldP;
      GreyLevelImage::value_type* currNewP = startNewP;

      // Where do we go from here?
      int nextPosition = 0;
      int d34 = 0; // (3,4) distance between start point and current point
      loadBinaryNeighbors(currNewP, n);

      for (int k = 1 ; (k < 9) && (nextPosition == 0) ; ++k)
        {
          if (n[k] != 0)
            {
              nextPosition = k;
            }
        }

      int iCameFrom = 0;

      // Go to next point
      switch (nextPosition)
        {
        case 1:
          --currOldP;
          --currNewP;
          d34 += 3;
          iCameFrom = 5;
          break;
        case 2:
          currOldP -= _width + 1;
          currNewP -= _width + 1;
          d34 += 4;
          iCameFrom = 6;
          break;
        case 3:
          currOldP -= _width;
          currNewP -= _width;
          d34 += 3;
          iCameFrom = 7;
          break;
        case 4:
          currOldP -= _width - 1;
          currNewP -= _width - 1;
          d34 += 4;
          iCameFrom = 8;
          break;
        case 5:
          ++currOldP;
          ++currNewP;
          d34 += 3;
          iCameFrom = 1;
          break;
        case 6:
          currOldP += _width + 1;
          currNewP += _width + 1;
          iCameFrom = 2;
          d34 += 4;
          break;
        case 7:
          currOldP += _width;
          currNewP += _width;
          iCameFrom = 3;
          d34 += 3;
          break;
        case 8:
          currOldP += _width - 1;
          currNewP += _width - 1;
          iCameFrom = 4;
          d34 += 4;
          break;
        default:
          break; 
        } // END switch
          
      // if pruning condition does hold, start a pruning process
      if (((*currOldP > *startOldP) || // labeled more than startig point
           (((*currOldP + 2)/3) == ((*startOldP +2)/3))) && // or same layer
          ((*startOldP - *currOldP + d34) <= 3 * maxPruning))
        {
          // Potential pruning branch
          vector<GreyLevelImage::value_type*> possiblePointsToBePruned;

          // The starting point can be pruned
          possiblePointsToBePruned.push_back(startNewP);

          // Can we go on?
          int canContinue = 1;
          // Can we prune?
          int canPrune = 1;

          while (canContinue)
            {
              // compute degree of current point
              int mydeg = degree(currNewP);
              if (mydeg > 2)
                canContinue = 0;
              else if (mydeg < 2) // Do not prune, we have come to an end
                {
                  canContinue = 0;
                  canPrune = 0;
                }
              else // actual tracking, with mydeg == 2
                {
                  // Save position from which we are
                  GreyLevelImage::value_type* savNewP = currNewP;
                  // Apparently not used!
                  //GreyLevelImage::value_type* savOldP = currOldP;

                  // Find where we are going, without going back
                  // to where we came from
                  nextPosition = 0;
                  loadBinaryNeighbors(currNewP, n);
                  for (int k = 1 ; (k < 9) && (nextPosition == 0) ; ++k)
                    {
                    if ((n[k] != 0) && (k != iCameFrom))
                      {
                        nextPosition = k;
                      }
                    }
                  
                  // Go to next point
                  switch (nextPosition)
                    {
                    case 1:
                      --currOldP;
                      --currNewP;
                      d34 += 3;
                      iCameFrom = 5;
                      break;
                    case 2:
                      currOldP -= _width + 1;
                      currNewP -= _width + 1;
                      d34 += 4;
                      iCameFrom = 6;
                      break;
                    case 3:
                      currOldP -= _width;
                      currNewP -= _width;
                      d34 += 3;
                      iCameFrom = 7;
                      break;
                    case 4:
                      currOldP -= _width - 1;
                      currNewP -= _width - 1;
                      d34 += 4;
                      iCameFrom = 8;
                      break;
                    case 5:
                      ++currOldP;
                      ++currNewP;
                      d34 += 3;
                      iCameFrom = 1;
                      break;
                    case 6:
                      currOldP += _width + 1;
                      currNewP += _width + 1;
                      iCameFrom = 2;
                      d34 += 4;
                      break;
                    case 7:
                      currOldP += _width;
                      currNewP += _width;
                      iCameFrom = 3;
                      d34 += 3;
                      break;
                    case 8:
                      currOldP += _width - 1;
                      currNewP += _width - 1;
                      iCameFrom = 4;
                      d34 += 4;
                      break;
                    default:
                      break; 
                    }
                  // Does the pruning condition hold?
                  // ** Should we rather check from savOldP ???
                  if (((*currOldP > *startOldP) || 
                       (((*currOldP + 2)/3) == ((*startOldP +2)/3))) && 
                      ((*startOldP - *currOldP + d34) <= 3 * maxPruning))
                    // The previous point can definitively be pruned
                    {
                      possiblePointsToBePruned.push_back(savNewP);
                      mydeg = degree(currNewP);
                    }
                  else
                    {
                      canContinue = 0;
                      canPrune = 0;  // ** or should we prune what we can ?
                    }
      
                }
            }

          if (canPrune != 0)
            {
              // we actually can prune these points
              pointsToBePruned.insert(pointsToBePruned.end(), 
                                      possiblePointsToBePruned.begin(),
                                      possiblePointsToBePruned.end());
            }
          
          possiblePointsToBePruned.erase(possiblePointsToBePruned.begin(),
                                         possiblePointsToBePruned.end());
          
        }
      else
        {
          // if pruning condition does not hold, forget it
        }

      // Advance
      ++myIter;
    }

  // Now we have found all that can be pruned
  // Remove the corresponding points from

  vector<GreyLevelImage::value_type*>::iterator iterPoints =
    pointsToBePruned.begin();

  while (iterPoints != pointsToBePruned.end())
    {
      *(*iterPoints) = 255;
      ++iterPoints;
    }
}


// REMOVE JAGGEDNESS
// =================


// As noted by Sanniti di Baja in her article, we may want to make
// this condition more strict, by allowing the shift only if the
// resulting marked pixel belongs to the set of skeletal pixels
// identified during the process. This would need a special marking
// for those pixels.
// I think this is not necessary, but let us keep it in mind. (KT)

void
LabeledSkeletonImage::removeJaggedness()
{
  GreyLevelImage::value_type n[10];

  for (int i = 1 ; i < (_height-1) ; ++i)
    {
      GreyLevelImage::value_type* newp = _pixMap2 + (i * _width) + 1;
      GreyLevelImage::value_type* oldp = _pPixMap + (i * _width) + 1;

      for(int j = 1 ; j < (_width - 1) ; ++j, ++newp, ++oldp)
        {
          if ((*newp == 0) && (degree(newp) == 2))
            {
              loadBinaryNeighbors(newp, n);

              if ((n[2]==1) && (n[4]==1) && (*(oldp-_width)!=0))
                {
                  *newp          = 255;
                  *(newp-_width) = 0;
                }
              else if ((n[4]==1) && (n[6]==1) && (*(oldp+1)!=0))
                {
                  *newp     = 255;
                  *(newp+1) = 0;
                }
              else if ((n[6]==1) && (n[8]==1) && (*(oldp+_width)!=0))
                {
                  *newp          = 255;
                  *(newp+_width) = 0;
                }
              
              else if ((n[8]==1) && (n[2]==1) && (*(oldp-1)!=0))
                {
                  *newp     = 255;
                  *(newp-1) = 0;
                }
            } // END if
        } // END for j
    } // END for i
}



// LOAD THE RIGHT LABELS FOR A NEIGHBORHOOD
// ========================================


void
LabeledSkeletonImage::loadLabeledNeighbors(unsigned char* p, unsigned char* n)
{
  *n = *p;
  ++n;
  *n = *(p - 1);
  ++n;
  *n = *(p - _width - 1);
  ++n;
  *n = *(p -_width);
  ++n;
  *n = *(p -_width + 1);
  ++n;
  *n = *(p + 1);
  ++n;
  *n = *(p + _width + 1);
  ++n;
  *n = *(p + _width);
  ++n;
  *n = *(p + _width - 1);
  ++n;
  *n = *(p - 1);  // n[9] == n[1]
}


// LOAD THE BINARY IMAGE OF A NEIGHBOURHOOD
// ========================================


void
LabeledSkeletonImage::loadBinaryNeighbors(unsigned char* p, unsigned char* n)
{
  *n = *p ? 0 : 1;
  ++n;
  *n = *(p-1) ? 0 : 1;
  ++n;
  *n = *(p - _width - 1) ? 0 : 1;
  ++n;
  *n = *(p - _width) ? 0 : 1;
  ++n;
  *n = *(p - _width + 1) ? 0 : 1;
  ++n;
  *n = *(p + 1) ? 0 : 1;
  ++n;
  *n = *(p + _width + 1) ? 0 : 1;
  ++n;
  *n = *(p + _width) ? 0 : 1;
  ++n;
  *n = *(p + _width - 1) ? 0 : 1;
  ++n;
  *n = *(p - 1) ? 0 : 1;  // n[9] == n[1]
}


// LOAD SUITABLE X BINARY NEIGHBORS
// ================================


void
LabeledSkeletonImage::loadSuitXBinary(unsigned char suitbin[], unsigned char n[])
{
  for (int i = 1 ; i < 10 ; ++i)
    {
      if ((n[i] < n[0]) && (n[i] != 255))
        {
          suitbin[i] = 1;
        }
      else
        {
          suitbin[i] = 0;
        }
    }
}


// LOAD SUITABLE C BINARY NEIGHBORS
// ================================


void
LabeledSkeletonImage::loadSuitCBinary(unsigned char suitbin[], unsigned char n[])
{
  for (int i = 1 ; i < 10 ; ++i)
    {
      if((n[i] > n[0]) && (n[i] != 255))
        {
          suitbin[i] = 1;
        }
      else
        {
          suitbin[i] = 0;
        }
    }
}


// NUMBER OF 8-CONNECTED COMPONENTS IN A NEIGHBORHOOD
// ==================================================


int
LabeledSkeletonImage::connexityCnt(unsigned char n[])
{
  int c8 = 0;

  for (int k = 1 ; k <= 4 ; ++k)
    {
      int k2 = 2 * k;
      int tmp = 1 - n[k2 - 1];
      c8 += tmp - (tmp * (1 - n[k2]) * (1 - n[k2 + 1]));
    }

  return c8;
}


// NUMBER OF 4-CONNECTED COMPONENTS IN A NEIGHBORHOOD
// ==================================================


int
LabeledSkeletonImage::crossingCnt(unsigned char n[])
{
  int x4 = 0;
  
  for(int k = 1 ; k <= 8 ; ++k)
    {
      x4 += std::abs(n[k + 1] - n[k]);
    }
  
  x4 /= 2;
  return x4;
}


// DEGREE OF A POINT: ITS NUMBER OF MARKED NEIGHBORS
// =================================================


int
LabeledSkeletonImage::degree(unsigned char* p)
{
  int d = 0;

  if(*(p-1) != 255)
    {
      ++d;
    }
  if(*(p - _width - 1) != 255)
    {
      ++d;
    }
  if(*(p - _width) != 255)
    {
      ++d;
    }
  if(*(p - _width + 1) != 255)
    { 
      ++d;
    }
  if(*(p + 1) != 255)
    {
      ++d;
    }
  if(*(p + _width + 1) != 255)
    {
      ++d;
    }
  if(*(p + _width) != 255)
    {
      ++d;
    }
  if(*(p + _width - 1) != 255)  
    {
      ++d;
    }

  return d;
}


// TRIPLE CONSECUTIVE NEIGHBORS ALL LABELED 3
// ==========================================


int
LabeledSkeletonImage::oddEvenOddLab3(unsigned char n[], int debut)
{
  return n[debut] + n[debut + 1] + n[debut + 2];
}


// GRADIENT IN A 3*3 NEIGHBORHOOD
// ==============================


void
LabeledSkeletonImage::gradient(unsigned char n[], float grad[])
{
  float* g = grad;
  *g++ = 0.0;

  unsigned char* pn = n + 1;

  for (int k = 1 ; k < 9 ; ++k)
    {
      *g++ = (float) (*pn++ - *n) / (((k % 2) == 0) ? 4.0 : 3.0);
    }
}


// POSITION OF MAXIMUM GRADIENT
// ============================


int
LabeledSkeletonImage::maxGradientPos1(unsigned char* p, float grad[])
{
  int ind = 0;
  float* g = grad + 1;

  for(int k = 1 ; k < 9 ; ++k, ++g)
    {
      if (*g > grad[0])
        {
          ind = k;
          grad[0] = *g;
        }
    }

  return ind;
}


// FIND, IF ANY, SECOND POSITION WITH SAME MAXIMUM GRADIENT
// ========================================================


int
LabeledSkeletonImage::maxGradientPos2(float grad[], unsigned char* p, int indmax)
{
  int ind=0;

  if ((*p == 3) || (*p == 4))
    {
      float* g = grad + indmax + 1;

      for(int k = indmax + 1 ; k < 9 ; ++k, ++g)
        {
          if(*g == grad[0])
            {
              ind = k;
            }
        }
    }

  return ind;
}


// REDUCE TO UNIT WIDTH
// ====================


// Erase marker from any marked pixel satisfying both:
// - at least one odd-neighbor is not marked,
// - at least a triplet of neighbors (nk, nk+2, nk+5)
//   (k odd, addition modulo 8) exists such that
//   nk and nk+2 are marked, with nk+5 not marked

void
LabeledSkeletonImage::unitWidth(unsigned char* p, unsigned char lab[])
{
  if (((lab[1] == 255) || (lab[3] == 255) || 
       (lab[5] == 255) || (lab[7] == 255))
      &&
      (((lab[1] == 0) && (lab[3] == 0) && (lab[6] == 255)) ||
       ((lab[3] == 0) && (lab[5] == 0) && (lab[8] == 255)) ||
       ((lab[5] == 0) && (lab[7] == 0) && (lab[2] == 255)) ||
       ((lab[7] == 0) && (lab[1] == 0) && (lab[4] == 255))))
    {
      *p = 255;
    }
}


// -------------------------------------------------------------------

} // namespace qgar