zernike.cpp

 
//////////////////////////////////////////////////////////////////////////
//
//
//                              zernike.cpp
//
//
//                           Michael Boland
//                          21 February 1996 
//
//
// Revisions:
//    13Jan96:  Added test code for use in determining z=1 for Zernike
//               calculations. -M. Boland
//    19Jan96:  Moved cell_zernike() to cell_z.cpp -M. Boland
//    21Feb96:  Moved cell_z.cpp to zernike.cpp. -M. Boland
//    28Feb96:  Added functions to support the external calculation of 
//                features for the Cell or for Objects. -M. Boland
//    08Mar96:  Fixed behavior for case where total intensity = 0
//                - M. Boland
//    15Jul96:  Added ability to specify the order of Zernike moments
//                to be calculated and to specify the radius of 
//                integration. -M. Boland
//
//
//////////////////////////////////////////////////////////////////////////

#include 
#include 
#include 
#include 
#include 

#include "cell.h"
#include "object.h"
#include "zernike.h"


//
// Calculates n! (uses double arithmetic to avoid overflow)
//
double factorial(long n)
{
	if(n < 0)
		return(0.0) ;
	if(n == 0)
		return(1.0) ;
	else
		return(n * factorial(n-1)) ;
}	


//
// Equations for Zernike moments from Prokop and Reeves, CVGIP: Graphical
//   Models and Image Processing, v 54 n 5, Sep 92, pp 438-460.
//

//
// Function to calculate the Zernike polynomial Rnl(r).
//  Note: 0 <= r <= 1
//

double zernikeR(int n, int l, double r)
{
  int m ;
  double sum = 0.0 ;
  
  if( ((n-l) % 2) != 0 )
    {
      cerr << "zernikeR(): improper values of n,l\n" ;
      return(0) ;
    }
  
  for(m = 0; m <= (n-l)/2; m++)
    {
      sum += (pow((double)-1.0,(double)m)) * ( factorial(n-m) ) / 
	( factorial(m) * (factorial((n - 2*m + l) / 2)) *
	  (factorial((n - 2*m - l) / 2)) ) *
	( pow((double)r, (double)(n - 2*m)) );
    }
  
  return(sum) ;
}

//
// Function for the Zernike moments.
//  Notes: requires xbar and ybar to perform translation and r_max to
//    normalize rho to fall between 0 and 1.
//

complex zernikeZ(Pixel pixel_array[], long num_pixels, 
			 double total_color_intensity, int n, int l, 
			 int color, double xbar, double ybar, double r_max)
{
  double rho ;		// radius of pixel from COM
  double theta ;		// angle of pixel
  complex integral(0,0) ;
  // summation of zernike calculations over all pixels 
  double x, y ;		// pixel coordinates
  double intensity ;	// pixel intensity
  
  if (!pixel_array)
    {
      cerr << "zernikeZ(): no pixels.\n" ;
      return(-1.0) ;
    }
  
  for(int i = 0; i < num_pixels; i++)
    {
      x = (double)(pixel_array[i].x) - (double)xbar ;
      y = (double)(pixel_array[i].y) - (double)ybar ;
      
      //
      // Normalize the intensity of each pixel to the total intensity
      //
      
      if(total_color_intensity > 0.0)
	{
	  intensity = (double)(pixel_array[i].intensity[color]) / 
	    (double)(total_color_intensity) ;
	}
      else
	{
	  intensity = 0.0 ;
	}

      //
      // Normalize the radial distance to each pixel to be <= 1
      // 
      
      if(r_max > 0.0)
	rho = (double)sqrt(x*x + y*y) / (double)r_max ;
      else
	rho = 0.0 ;
      
      if(rho <= 1.0)
	{
	  if(y == 0.0 && x == 0.0)
	    theta = 0.0 ;
	  else
	    theta = atan2(y, x) ;
	  integral += conj(zernikeR(n,l,rho) * polar(1.0, l*theta)) *
	    // 
	    // Confirm this formula
	    //	    intensity ;
	    intensity * rho ;
	}
    } // end for
  return(integral * (n+1)/PI) ;
}


int zernike(Moments* moments_struct, Pixel pixel_array[], long z_order,
	    long num_pixels, long radius, int color, 
	    double total_color_intensity, ofstream& outfile)
{
  // 
  // Number of features generated by Zernike moments up to order z_order.
  //
  long   z_num_features = (long)( ((z_order + 4) * 
				   (z_order + 1) - 2 * 
				   (long)(((long)z_order + 1) / (long)2) ) / 4 ) ;
  //
  // Storage for new features
  //
  double* zfeatures = new double[z_num_features + 1] ;
  int     zindex = 0 ;
  double  r_max ;
  
  if (z_order > 20 || z_order < 0)
    {
      cerr << "zernike(): You choice of z_order is invalid "
	   << "choose a value between 0 and 20.\n" ;
      exit(1) ;
    }
  if (!pixel_array)
    {
      cerr << "zernike(): No objects.\n" ;
      return(0) ;
    }
  
  if(radius == 0)
    {
      //
      // rog = radius of gyration about the com
      //
      double rog ;
      
      //
      // Radius of Gyration
      //
      // Formula from Prokop RJ, Reeves AP, CVGIP:Graphical Models and Image
      //  Processing, v54 n5, Sep 92, pp. 438-460.
      //
      
      if(moments_struct->mu00 > 0.0)
	rog = sqrt((moments_struct->mu20 + moments_struct->mu02) 
		   / moments_struct->mu00) ;
      else
	rog = 0.0 ;
      
      r_max = 2.0 * rog ;
    }
  else
    {
      r_max = (double)radius ;
    }
  
  
  //
  // Loop to calculate Zernike features
  //
  
  for (int n = 0; n <= z_order; n++)
    {
      for (int l = 0; l <= n; l++)
	if ((n-l) % 2 == 0)
	  {
	    zindex++ ;
	    zfeatures[zindex] = abs(zernikeZ(pixel_array, num_pixels,
					     total_color_intensity, n, l, color, 
					     moments_struct->xbar, 
					     moments_struct->ybar, r_max)) ;
	  }
    }
  
  //
  // Output of zernike features
  //
  
  outfile << "zernike_color_" << color << "\t" ;
  
  for (int i = 1; i <= z_num_features; i++) 
    {
      outfile << zfeatures[i] << "\t" ;
    }
  
  delete[] zfeatures ;
  return(1) ;
}