Properties of chain code i.e. perimeter, area and center of mass

The chain code is compact and efficient to store and transfer. Though it is not that straight forward to calculate some properties of shape .i.e area, perimeter and center of mass.

Perimeter is relatively easy and from Bernd Jähne book ( Digital Image Processing 5th: 508-511 ), it can be calculated as below code.

  double CalculatePerimeter( const ChainCode& chains )
  {
    double p = 0;
    double sqrt2 = sqrt( 2.0 );

    for( size_t i=0; i < chains.size(); ++i )
    {
      if( chains[i]%2 == 0 ) 
      {
        p += 1.0;
      }
      else
      {
        p += sqrt2;
      }
    }

    return p;
  }

Bernd Jähne also showed how to get the area using a chain code but it doesn't work well when the chain code is 1 pixel width or height e.g. ( 0 0 4 4 ) will return 2 instead of 3. Refer below image how the area end up 1 pixel short.

Cris Luengo ( Cris's Image Analysis Blog : More chain code measures ) shows a matlab code that works well with 1 pixel wide chain code. In above case of 3 pixels, it adds two more addition when the direction changes which add up 1 more pixel. Refer below image.

Here is the C/C++ code snippet that is converted from the above Cris Luengo's matlab code.

// making a circular chain code by inserting the last code to the front 
ChainCode cc;
cc.push_back( chains[ chains.size()-1 ] );
for( size_t i=0; i < chains.size(); ++i )
{
  cc.push_back( chains[i] );
}
m_Area = CalculateArea( cc );


int CalculateArea( const ChainCode& cc ) // a circular chain
{
  int (&M)[8][8]( TableAreaIncrement );
  int B = 10, A = 0;

  for( size_t i=1; i < cc.size(); ++i )
  {
    uint8_t c_i = cc[i];
    uint8_t c_im1 = cc[i-1];

    switch( M[c_im1][c_i] )
    {
    case 'X' :
      throw Exception( StringFormat::As( _T("Invalid chain code found as %d->%d"), 
        c_im1, c_i ));
    case 'A': A += -B+1;  break;
    case 'B': A += -B;    break;
    case 'C': A += B;    break;
    case 'D': A += B-1;    break;
    }

    switch( c_i )
    {
    case 0: A+=B;      break;
    case 1: A+=(++B);    break;
    case 2: B++;      break;
    case 3: A+=-(++B)+1;  break;
    case 4: A+=-B+1;    break;
    case 5: A+=-(--B)+1;  break;
    case 6: B--;      break;
    case 7: A+=(--B);    break;
    }
  }

  return A;
}
  
int ChainCodeProperties::TableAreaIncrement[8][8] = 
{
  {'0','1','B','X','A','A','6','7'},
  {'0','1','B','B','X','A','6','7'},
  {'C','C','2','3','4','X','C','C'},
  {'C','C','2','3','4','5','X','C'},
  {'C','C','2','3','4','5','D','X'},
  {'X','C','2','3','4','5','D','D'},
  {'0','X','A','A','A','A','6','7'},
  {'0','1','X','A','A','A','6','7'},
};

The center of mass can be calculated by summing up x and y coordinates. But just simply adding up does not cover 1 pixel wide chain code as below image. The (1,1) pixel added two times instead one time. The compensation can be done by adding more pixels when a direction changes in opposite way as the idea of Cris Luengo's area calculation.

Below is the C/C++ code snippet that calculate the centroid.

void CalculateCentroid( const ChainCode& cc, double* cx, double *cy )
{
  int (&O)[8][2]( TableOffsetIncrement );
  int vpc = 0;      // visitied pixel count
  int sx = 0, sy = 0;    // sum of x and y
  int px = 0, py = 0;
  
  for( size_t i=1; i < cc.size(); ++i )
  {
    uint8_t c_i = cc[i];
    uint8_t c_im1 = cc[i-1];

    if( c_i != c_im1 && (c_i & 0x3 ) == ( c_im1 & 0x3 ) )
    {
      // when direction changes in opposite direction 
      // e.g. 0 -> 4, 4 -> 0, 1 -> 5, ...
      vpc++; sx += px;  sy+=py;
    }

    px += O[ c_i ][0];
    py += O[ c_i ][1];

    vpc++; sx += px;  sy+=py;
  }

  *cx = (double)sx / (double)vpc;
  *cy = (double)sy / (double)vpc;
}

int ChainCodeProperties::TableOffsetIncrement[8][2] =
{
  { 1,  0  },
  { 1,  -1  },
  { 0,  -1  },
  { -1,   -1  },
  { -1,  0  },
  { -1,  1 },
  { 0,  1 },
  { 1,  1 },
};  

It works well most of the case but it fails if there is a pixel which is visited more than 2 times as below image. It has (1,1) pixel counted 3 times but it should be counted 2 times.

Contour tracing by Moore Neighbor Algorithm

Having a binary image with defects, it is time to find out how the defects looks like. The location and size is what we want to figure out. Though we need to extract the defect from the image first. This is called segmentation and there are two methods. One is contour detection and the other is region coloring.

Region coloring is to put a unique label to a separated component. At first it labels all the pixels in a way and at a second phase, those labels are merged when connected. The result of the region coloring is still a area image with each component assigned with a unique number. And if you want to count how many defects are there, the area has to be scanned and count unique number. Or if you want to transmit the defect information, the whole labeled image has to be sent out - though the size can be much smaller than original image if run length encoding is used.

Contour detection is to trace the defect pixels and extract contour only. This start on a pixel and keep following on nearby defect pixels until it comes back to the start pixel. The result of this segmentation is only contour pixels and can apply chain coding which takes 3 bits per pixel. And it makes quite small size and compact to transmit to other. This contour information can be used to calculate area or moment for analysis.

There are number of algorithms for contour tracing. Refer [CTA]. The square tracing doesn't work with 8 neighbor hood. The Moore-Neighbor tracing and radial sweep algorithms are almost identical. Theo Pavlidis's algorithm is not that easy to implement.

In here, I will show implementation of Moore-Neighbor Tracing based on the above web site. Though the stop condition at the web site is not so simple and will use approach from other book - P192 of [IPA] as below exerpt.

  If the current boundary element Pn is equal to the second border element P1 , and
  if the previous border element Pn-1 is equal to Po , stop. Otherwise repeat step 2.

First the image has to be traversed to find the start point. And when found a pixel, do the contour tracing and clear out the traced object so that it does not get picked up again. Here is the code snippet.

void Build( BinaryImageBase& image, const Region& roi, PointsContainer* contours ) 
{
  for( int sy = roi.Top; sy<= roi.Bottom; ++sy )
  {
    for( int sx = roi.Left; sx<=roi.Right; ++sx )
    {
      if( image.GetPixel( sx, sy ) == false )
        continue;

      Points contour;
      m_Tracer->Trace( image, roi, Point((uint16_t)sx,(uint16_t)sy), &contour ); // MooreNeighbor Tracing
      m_Filler->Fill( image, contour[0], false );  // FloodFill
    
      contours->push_back( contour );
    }
  }
}

The Moore-Neighbor tracing can be implemented as below.

struct WalkTableItem
{
  Point Step;
  int Next;
} WalkTable[8] =
{
  { Point(1,0),   6, },  // from North to North-East, when found pixel, start from index 6 of this table
  { Point(0,1),  0, },   // from North-East to East
  { Point(0,1),  0, },   // from East to South-East
  { Point(-1,0),  2, },  // from South-East to South
  { Point(-1,0),  2, },  // from South to South-West
  { Point(0,-1),  4, },  // from South-West to West
  { Point(0,-1),  4, },  // from West to North-West
  { Point(1,0),  6, },   // from North-West to North
};

void Trace( const BinaryImageBase& image, const Region& roi, const Point& start, Points* trace )
{
  const int SizeOfT = sizeof(WalkTable)/sizeof(WalkTable[0]);
  WalkTableItem (&T)[SizeOfT]( WalkTable );   // access WalkTable with variable name T for short

  Points B; B.push_back( start );             // B is the container of points
  int sk = 0, k = 0;                          // sk is the start index, k is current index
  Point c = start - T[k].Step;                // c is current point
  k = T[k].Next;
  sk = ( SizeOfT + k - 1 ) % SizeOfT;

  for( ; ; ) 
  {
    Point n = c + T[k].Step;                  // n is next point
    if( image.GetPixel( n.As.XY.X, n.As.XY.Y ) == true ) 
    {
      size_t bc = B.size();

      if( bc > 2 && n == B[1] && B[bc-1] == B[0] )  
      {
        break;                                 // end condition by [IPA]
      }
      else
      {
        B.push_back( n );                      // found contour pixel
        n = c;                                 // back-track
        k = T[k].Next;                         // next index to search 
        sk = ( SizeOfT + k - 1 ) % SizeOfT;    // stop index if no pixel at moore-neighbor
      }
    }
    else
    {
      k = (k+1) % SizeOfT;                      // move to next index 
    }

    c = n;
    if( sk == k )                               // stop condition with 1 pixel defect
    {
      break;
    }
  }

  if( B.size() == 1 )                           
  {
    trace->swap( B );                           // 1 pixel defect
  }
  else
  {
    trace->swap( Points( &B[0], &B[ B.size() - 1 ] ) );  // The last pixel is the start pixel so exclude it.
  }
}

Here is code of stack-based FloodFill based on [WFF].

void Fill( BinaryImageBase& image, const Point& start, bool fillValue )
{
  std::stack< Point > Q;
  Q.push( start );

  while( Q.size() != 0 )
  {
    Point p = Q.top();
    Q.pop();

    if( p.As.XY.X < 0 || p.As.XY.X >= image.GetWidth() ||
      p.As.XY.Y < 0 || p.As.XY.Y >= image.GetHeight() )
      continue;

    if( image.GetPixel( p.As.XY.X, p.As.XY.Y ) == fillValue ) 
      continue;

    image.SetPixel( p.As.XY.X, p.As.XY.Y, fillValue );

    Q.push( Point( p.As.XY.X + 1, p.As.XY.Y ) );
    Q.push( Point( p.As.XY.X - 1, p.As.XY.Y ) );

    Q.push( Point( p.As.XY.X + 1, p.As.XY.Y - 1 ) );
    Q.push( Point( p.As.XY.X,     p.As.XY.Y - 1 ) );
    Q.push( Point( p.As.XY.X - 1, p.As.XY.Y - 1 ) );

    Q.push( Point( p.As.XY.X + 1, p.As.XY.Y + 1 ) );
    Q.push( Point( p.As.XY.X,     p.As.XY.Y + 1 ) );
    Q.push( Point( p.As.XY.X - 1, p.As.XY.Y + 1 ) );
  }
}

The result of contour tracing is boundary pixels and one or more pixels can be repeated as it needs to return to the original image. Refer below image and it's result.

Contour pixels : (1,1),(2,2),(3,1),(2,2),(3,3),(2,2)

Reference :
[CTA] Contour Tracing WEB Site by Abeer George Ghuneim
[IPA] Image Processing, Analysis and Machine Vision by Milan Sonka, Vaclav Hlavac, Roger Boyle.
[WFF] FloodFill from Wikipedia