moments.cpp

 

//////////////////////////////////////////////////////////////////////////
//
//
//                            moments.cpp
//
//
//                           Michael Boland
//                             21 May 1996 
//
//
//////////////////////////////////////////////////////////////////////////

#include 
#include 

#include "moments.h"

Moments* moments(Pixel pixel_array[], long num_pixels, int color, 
                  double total_color_intensity)
{
	int x, y, z ;
	double intensity ;

	Moments* moments_struct = new Moments ;
	if (!moments_struct) return (0) ;

    moments_struct->fxyz    = 0.0 ;
    moments_struct->xfxyz   = 0.0 ;
    moments_struct->yfxyz   = 0.0 ;
    moments_struct->x2fxyz  = 0.0 ;
    moments_struct->y2fxyz  = 0.0 ;
    moments_struct->x3fxyz  = 0.0 ;
    moments_struct->y3fxyz  = 0.0 ;
    moments_struct->xyfxyz  = 0.0 ;
    moments_struct->xy2fxyz = 0.0 ;
    moments_struct->x2yfxyz = 0.0 ;

    //
    // Note:  All formulas are from "Digital Image Processing"
    //          by Gonzalez and Woods (c) 1992. pp. 514-516.
    //
    // 2D Moments
    //

    for (int i = 0; i < num_pixels; i++)
    {
        x = pixel_array[i].x ;
        y = pixel_array[i].y ;
        z = pixel_array[i].z ;
        //
        // Use total cell intensity to normalize each pixel
        //
        if(total_color_intensity > 0.0)
            intensity = (double)(pixel_array[i].intensity[color]) /
                        (double)(total_color_intensity) ;
        else
            intensity = 0.0 ;

        moments_struct->fxyz += intensity ;

        moments_struct->xfxyz += x * intensity ;
        moments_struct->yfxyz += y * intensity ;
        moments_struct->x2fxyz += pow(x, 2.0) * intensity ;
        moments_struct->y2fxyz += pow(y, 2.0) * intensity ;

        moments_struct->x3fxyz += pow(x, 3.0) * intensity ;
        moments_struct->y3fxyz += pow(y, 3.0) * intensity ;

        moments_struct->xyfxyz += x * y * intensity ;
        moments_struct->xy2fxyz += x * pow(y, 2.0) * intensity ;
        moments_struct->x2yfxyz += pow(x, 2.0) * y * intensity ;

    } // end for

    //
    // Central Moments
    //

    if(moments_struct->fxyz > 0.0)
    {
        moments_struct->xbar = moments_struct->xfxyz / 
	  moments_struct->fxyz ;
        moments_struct->ybar = moments_struct->yfxyz / 
	  moments_struct->fxyz ;
    }
    else
    {
        moments_struct->xbar = 0 ;
        moments_struct->ybar = 0 ;
    }

    moments_struct->mu00 = moments_struct->fxyz ;
    moments_struct->mu20 = moments_struct->x2fxyz 
      - (moments_struct->mu00 * pow(moments_struct->xbar,2.0)) ;
    moments_struct->mu02 = moments_struct->y2fxyz 
      - (moments_struct->mu00 * pow(moments_struct->ybar,2.0)) ;
    moments_struct->mu11 = moments_struct->xyfxyz 
      - (moments_struct->mu00 * moments_struct->xbar * 
	 moments_struct->ybar) ;
    moments_struct->mu30 = moments_struct->x3fxyz 
      - (3 * moments_struct->x2fxyz * moments_struct->xbar)
      + (2 * moments_struct->mu00 * pow(moments_struct->xbar, 3.0)) ;
    moments_struct->mu21 = moments_struct->x2yfxyz 
      - (moments_struct->x2fxyz * moments_struct->ybar) 
      - (2 * moments_struct->xyfxyz * moments_struct->xbar)
      + (2 * moments_struct->mu00 * pow(moments_struct->xbar, 2.0) 
	 * moments_struct->ybar) ;
    moments_struct->mu12 = moments_struct->xy2fxyz - 
      (moments_struct->y2fxyz * moments_struct->xbar) 
      - (2 * moments_struct->xyfxyz * moments_struct->ybar) + 
      (2 * moments_struct->mu00 * moments_struct->xbar 
       * pow(moments_struct->ybar, 2.0)) ;
    moments_struct->mu03 = moments_struct->y3fxyz - 
      (3 * moments_struct->y2fxyz * moments_struct->ybar) 
      + (2 * moments_struct->mu00 * pow(moments_struct->ybar, 3.0)) ;

	return(moments_struct) ;
}