 /*
  * Khoros: $Id: lvshape.c,v 1.1 1991/05/10 15:41:54 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvshape.c,v 1.1 1991/05/10 15:41:54 khoros Exp $";
#endif

 /*
  * $Log: lvshape.c,v $
 * Revision 1.1  1991/05/10  15:41:54  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as too the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvshape.c
 >>>>
 >>>>      Program Name: vshape
 >>>>
 >>>> Date Last Updated: Tue Mar 26 15:36:45 1991 
 >>>>
 >>>>          Routines: lvshape - the library call for vshape
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#include "vipl/lvshape.h"
#define         ALL_OPTS        61
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvshape - library call for vshape
*
* Purpose:
*    
*    Performs shape analysis on images.
*    
*    
* Input:
*    
*         image -- xvimage structure
*    
*         image_type -- if 0, image containing  clusters  or  labelled
*         objs. if 1, grey level image.
*    
*         axis_flg -- if 1, generates the image with axis  representa-
*         tion.
*    
*         outline_flg -- if  1,  generates  the  image  with  outlines
*         representation.
*    
*         file_flg --  if 1, generates the ASCII  file  with  objs  or
*         image information.
*    
*         input_file -- name of input file being analyzed
*    
*         printdev --  output device for printing
*    
*         print_flg -- print to printdev if set.
*    
*         standard: if 1, computes and prints standard moments.
*    
*         central: if 1, computes and prints central moments.
*    
*         invariants: if 1, computes and prints invariant moments.
*    
*    
* Output:
*    
*         image -- xvimage structure, holds the resulting image.
*    
*         axis -- xvimage structure, axis image.
*    
*         outline -- xvimage structure, outline image.
*    
*    
*
* Written By: Pascal ADAM, Shirley Lee
*    
*    Ghosted by  Pascal  ADAM   Fri  Jun  15  13:38:24  MDT  1990
*    Modified  vshape  to  output VIFF file format by Shirley Lee
*    Feb/26/1991.
*    
*    
****************************************************************/


/* -library_def */
int lvshape(image_type,
             image,
             typ_axe,
             typ_color,
             axis_flg,
             axis,
             outline_flg,
             outline,
             viff_flg,
             img3,
             file_flg,
             input_file,
             printdev,
             standard,
             central,
             invariant)

struct xvimage *image, **outline, **axis, **img3; 
int image_type, file_flg, axis_flg, outline_flg, viff_flg;
int typ_axe, typ_color, standard, central, invariant;
char *input_file;
FILE *printdev;

/* -library_def_end */

/* -library_code */
{
struct xvimage *createimage(),
               *o1_img,
               *o2_img,
               *o3_img;
int  nc, 
     nr, 
     color,
     row, 
     col, 
     row_deb, 
     row_end, 
     col_deb, 
     col_end, 
     ind, 
     numreg, 
     num;

char *dummy = '\0';

int *ptr, 
    *ptr_line, 
    *ptr_axe;

double pixel, 
       a1, 
       a2, 
       meters_surf;

char *malloc();
void draw_line();
void comp_shape();
char *unload_vector();
int printout();

int i;
float **tmp_obj;

XVSHAPE  **array_obj;



   nc = image->row_size;
   nr = image->col_size;

   /* Allocate space for output image (axis image) */

   o2_img=NULL;
   if (axis_flg == TRUE) 
   {
      o2_img  = createimage(nr,
                          nc, 
                          VFF_TYP_4_BYTE,
                          1,
                          1,
                          dummy, 
                          0, 
                          0, 
                          VFF_MS_NONE,
                          VFF_MAPTYP_NONE, 
                          VFF_LOC_IMPLICIT, 
                          0);

      if (o2_img  == NULL) 
      {
         (void) fprintf(stderr,"lvshape: Unable to allocate new image!\n\n");
         exit(1);
      }
   }

   /* Allocate space for output image (Outlines image) */

   o1_img = NULL;
   if (outline_flg == TRUE)
   {
      o1_img = createimage(nr,
                            nc, 
                            VFF_TYP_4_BYTE, 
                            1,
                            1,
                            dummy, 
                            0, 
                            0, 
                            VFF_MS_NONE,
                            VFF_MAPTYP_NONE, 
                            VFF_LOC_IMPLICIT, 
                            0);

      if (o1_img == NULL) 
      {
         (void) fprintf(stderr,"lvshape: Unable to allocate new image!\n\n");
         exit(1);
      }
   }


   /* assign image data address to ptr.................................. */
   ptr = (int  *) image->imagedata;

   /* computes the square meters of a pixel with info: pixsizx, pixsizy  */
   meters_surf = (double) (image->pixsizx * image->pixsizy); 

   /* 2 differents cases can occur...................................... */
   /*  -  1.  the user had a grey level as input .................... */
   /*         In this case the moment mpq are computed with respect to   */
   /*         the grey level of the pixel at this point. cf vshape.h     */
   /*  -  2.  The user had an image in which the value of the pixel      */
   /*         represents a class (cluster image or labelled image).      */
   /*         This time the algorithm will compute the moments... for    */
   /*         each class.  ............................................. */
   /*                                                                    */
   /* ____ These different cases are processed beyond .................. */
   
   /* Computation of mpq for a grey level image ........................ */
   if (image_type == 0) { /* grey level image  ......................... */

      numreg = 1;  /* in this case there is only one region ............ */

      /* allocate memory for numreg XVSHAPE pointed by array_obj ... */
      array_obj= (XVSHAPE  **) malloc(sizeof(XVSHAPE *)*numreg);
      if (array_obj == NULL) {
         (void)fprintf(stderr,
           "lvshape:  No space for array_obj _ malloc failed! \n");
         exit(1);
          }
      for (ind = 0; ind < numreg; ind++) {
         array_obj[ind] = (XVSHAPE *)malloc(sizeof(XVSHAPE));
         if (array_obj[ind] == NULL) {
            (void)fprintf(stderr,
                "lvshape:  No space for array_obj[%d] _ malloc failed! \n",
                ind);
            exit(1);
            }
         }

      /* initialisation ................................................ */
      array_obj[0]->right = array_obj[0]->bottom = 0;
      array_obj[0]->left = nc;
      array_obj[0]->top  = nr;

      for(row = 0; row < nr; row++){
         for(col = 0; col < nc; col++){
            pixel = (double)ptr[row * nc + col];
            a1    = pixel * (double)col;
            a2    = pixel * (double)row;
            array_obj[0]->m00 += pixel; 
            array_obj[0]->m01 += a1; 
            array_obj[0]->m10 += a2;
            array_obj[0]->m11 += a1 * (double)row; 
            a1   *= (double)col;
            a2   *= (double)row; 
            array_obj[0]->m02 += a1; 
            array_obj[0]->m20 += a2; 
            array_obj[0]->m12 += a1 * (double)row; 
            array_obj[0]->m21 += a2 * (double)col; 
            array_obj[0]->m03 += a1 * (double)col; 
            array_obj[0]->m30 += a2 * (double)row; 
            }
         }  
      }

   /* Computation of mpq for a cluster or labelled image ............... */
   else {  
     
      /* scan image to find number of different regions or clusters .... */
      numreg = 0;
      for(ind = 0; ind < nc*nr; ind++) {
         if (ptr[ind] > numreg) 
           numreg = ptr[ind];
      }

      numreg++;  /* increment numreg to get the total number of reg ... */
 
      /* allocate memory for numreg XVSHAPE pointed by array_obj ... */
      array_obj=(XVSHAPE  **)malloc(sizeof(XVSHAPE *)*numreg); 
      if (array_obj == NULL) {
         (void)fprintf(stderr,
           "lvshape:  No space for array_obj _ malloc failed! \n");
         exit(1);
          }

      for(ind = 0; ind < numreg; ind++) {
         array_obj[ind]=(XVSHAPE *)malloc(sizeof(XVSHAPE));

         if (array_obj[ind] == NULL) {
            (void)fprintf(stderr,
                "lvshape:  No space for array_obj[%d] _ malloc failed! \n",
                 ind);
            exit(1);
            }
        
         /* initialisation ............................................. */
         array_obj[ind]->right = array_obj[ind]->bottom = 0;
         array_obj[ind]->left = nc;
         array_obj[ind]->top  = nr;
         }

      /* scan the image and compute mpq for each regionodr cluster  */
      for(row = 0; row < nr; row++){
         for(col = 0; col < nc; col++){

            ind = ptr[row * nc + col]; /* ind contains the class number .. */

            /* gets the top, bottom, right and left border coord of obj. */
            if (array_obj[ind]->right<col) array_obj[ind]->right=col;
            if (array_obj[ind]->left>col) array_obj[ind]->left=col;
            if (array_obj[ind]->bottom<row) array_obj[ind]->bottom=row;
            if (array_obj[ind]->top>row) array_obj[ind]->top=row;

        
            array_obj[ind]->m00 += 1.0; 

            a1 = (double)col; a2 = (double)row;
            array_obj[ind]->m01 += a1;
            array_obj[ind]->m10 += a2; 
            array_obj[ind]->m11 += a1 * a2; 
            a1 *= a1; a2 *= a2;
            array_obj[ind]->m02 += a1; 
            array_obj[ind]->m20 += a2; 
            array_obj[ind]->m12 += (double)row * a1; 
            array_obj[ind]->m21 += (double)col * a2;
            array_obj[ind]->m03 += (double)col * a1;
            array_obj[ind]->m30 += (double)row * a2; 
            }
         }
      }
      
   /* In both case we know the information on mpq .................. */
   /* Let's compute all the other fields of the XVSHAPE ............ */

   (void)comp_shape( meters_surf,array_obj,numreg,central,invariant,nc, nr); 

   /* print out the results ........................................ */
   if (file_flg) printout(image_type, array_obj,numreg, standard, central, 
                          invariant, input_file, printdev);

   /* create output image files ONLY if image_type==1 (cannot represent */
   /* axis or outline when a grey level image is the input image....... */
   if ((image_type == 1) && (axis_flg == 1)) {
      ptr_axe = (int *)o2_img->imagedata;
      /* give to the background image its values ...................... */
      for (row = 0; row < nr ; row++) {
         for (col = 0; col < nc ; col++) {
            ptr_axe[row * nc + col] = 0;
            }
         }

      /* write  y&x axis of the found objs over the background image */ 
      for (ind = 1; ind < numreg; ind++) {
         if (typ_color == 0) color = 255; 
         else color = ind; 

         if (typ_axe == 1) {    
            draw_line(array_obj[ind]->axis_x1deb, 
                      array_obj[ind]->axis_y1deb,
                      array_obj[ind]->axis_x1end,
                      array_obj[ind]->axis_y1end,
                      ptr_axe,
                      nc,
                      nr,
                      color);

            draw_line(array_obj[ind]->axis_x2deb,
                      array_obj[ind]->axis_y2deb,
                      array_obj[ind]->axis_x2end,
                      array_obj[ind]->axis_y2end,
                      ptr_axe,
                      nc,
                      nr,
                      color);
            }
        
         else {
            row_deb = (int)(array_obj[ind]->top);
            row_end = (int)(array_obj[ind]->bottom);
            col_deb = (int)(array_obj[ind]->left);
            col_end = (int)(array_obj[ind]->right);
            for ( row = row_deb; row < row_end; row++) {
               for ( col = col_deb; col < col_end; col++) {
                  if ( (col == (int)(array_obj[ind]->ycent)) ||
                       (row == (int)(array_obj[ind]->xcent))  )
                     ptr_axe[row * nc + col] = color; 
                  }
               }
            }
         }
         *axis=o2_img;
      }
     
   if ((image_type == 1) && (outline_flg == 1)) {
      ptr_line = (int *)o1_img->imagedata;
      /* give to the background image its values .................. */
      for (row = 0; row < nr ; row++) {
         for (col = 0; col < nc ; col++) {
            ptr_line[row * nc + col] = 0;
            }
         }


      for ( row = 1; row < nr-1; row++){
         for ( col = 1; col < nc-1; col++){

            num = ptr[ row * nc + col];
            if (typ_color == 0) color = 255;
            else color = num;

            if ((ptr[ row * nc + col -1] != num) || 
                (ptr[ row * nc + col +1] != num) ||
                (ptr[(row -1)* nc + col] != num) ||
                (ptr[(row +1)* nc + col] != num)  )  
               ptr_line[row * nc + col] = color;
               
            }
         }
         *outline=o1_img;
      }

      if (viff_flg)  /* output in VIFF file format */ {

        o3_img = createimage((unsigned long) 1,      
                         (unsigned long) (numreg-1),        
                         (unsigned long) VFF_TYP_FLOAT,
                         (unsigned long) 1,      
                         (unsigned long) ALL_OPTS,      
                         dummy,                 
                         (unsigned long) 0,           
                         (unsigned long) 0,            
                         (unsigned long) VFF_MS_NONE,    
                         (unsigned long) VFF_MAPTYP_NONE, 
                         (unsigned long) VFF_LOC_IMPLICIT, 
                         (unsigned long) 0                
                        );

        if (o3_img == NULL)
        {
          (void) fprintf(stderr,"lvshape: Unable to allocate new image!\n\n");
          exit(1);
        }

        tmp_obj = (float **)calloc(numreg,sizeof(float*));

        if (tmp_obj == NULL) {
           (void)fprintf(stderr,
                "lvshape:  No space for tmp_obj - malloc failed! \n");
            exit(1);
        }

        for (i=0; i<numreg; i++) {
          tmp_obj[i] = (float *)calloc(ALL_OPTS,sizeof(float));

          if (tmp_obj[i] == NULL) {
            (void)fprintf(stderr,
                 "lvshape:  No space for tmp_obj[%d] - malloc failed! \n",
                 i);
             exit(1);
          }
        }

        /* change the array_obj[] struct to 2D array */

        for (i=0; i<(numreg-1); i++){
                tmp_obj[i][0] = (float)array_obj[i+1]->surf_image;
                tmp_obj[i][1] = (float)array_obj[i+1]->area;
                tmp_obj[i][2] = (float)array_obj[i+1]->percent_area;
                tmp_obj[i][3] = (float)array_obj[i+1]->real_area;
                tmp_obj[i][4] = (float)array_obj[i+1]->xcent;
                tmp_obj[i][5] = (float)array_obj[i+1]->ycent;
                tmp_obj[i][6] = (float)array_obj[i+1]->sigma_x;
                tmp_obj[i][7] = (float)array_obj[i+1]->sigma_y;
                tmp_obj[i][8] = (float)array_obj[i+1]->sigma_MX;
                tmp_obj[i][9] = (float)array_obj[i+1]->sigma_MY;
                tmp_obj[i][10] = (float)array_obj[i+1]->theta;
                tmp_obj[i][11] = (float)array_obj[i+1]->eccentricity;
                tmp_obj[i][12] = (float)array_obj[i+1]->m00;
                tmp_obj[i][13] = (float)array_obj[i+1]->m01;
                tmp_obj[i][14] = (float)array_obj[i+1]->m10;
                tmp_obj[i][15] = (float)array_obj[i+1]->m11;
                tmp_obj[i][16] = (float)array_obj[i+1]->m20;
                tmp_obj[i][17] = (float)array_obj[i+1]->m02;
                tmp_obj[i][18] = (float)array_obj[i+1]->m21;
                tmp_obj[i][19] = (float)array_obj[i+1]->m12;
                tmp_obj[i][20] = (float)array_obj[i+1]->m03;
                tmp_obj[i][21] = (float)array_obj[i+1]->m30;
                tmp_obj[i][22] = (float)array_obj[i+1]->mu00;
                tmp_obj[i][23] = (float)array_obj[i+1]->mu10;
                tmp_obj[i][24] = (float)array_obj[i+1]->mu01;
                tmp_obj[i][25] = (float)array_obj[i+1]->mu20;
                tmp_obj[i][26] = (float)array_obj[i+1]->mu02;
                tmp_obj[i][27] = (float)array_obj[i+1]->mu11;
                tmp_obj[i][28] = (float)array_obj[i+1]->mu30;
                tmp_obj[i][29] = (float)array_obj[i+1]->mu03;
                tmp_obj[i][30] = (float)array_obj[i+1]->mu21;
                tmp_obj[i][31] = (float)array_obj[i+1]->mu12;
                tmp_obj[i][32] = (float)array_obj[i+1]->nu00;
                tmp_obj[i][33] = (float)array_obj[i+1]->nu10;
                tmp_obj[i][34] = (float)array_obj[i+1]->nu01;
                tmp_obj[i][35] = (float)array_obj[i+1]->nu20;
                tmp_obj[i][36] = (float)array_obj[i+1]->nu02;
                tmp_obj[i][37] = (float)array_obj[i+1]->nu11;
                tmp_obj[i][38] = (float)array_obj[i+1]->nu30;
                tmp_obj[i][39] = (float)array_obj[i+1]->nu03;
                tmp_obj[i][40] = (float)array_obj[i+1]->nu21;
                tmp_obj[i][41] = (float)array_obj[i+1]->nu12;
                tmp_obj[i][42] = (float)array_obj[i+1]->M1;
                tmp_obj[i][43] = (float)array_obj[i+1]->M2;
                tmp_obj[i][44] = (float)array_obj[i+1]->M3;
                tmp_obj[i][45] = (float)array_obj[i+1]->M4;
                tmp_obj[i][46] = (float)array_obj[i+1]->M5;
                tmp_obj[i][47] = (float)array_obj[i+1]->M6;
                tmp_obj[i][48] = (float)array_obj[i+1]->M7;
                tmp_obj[i][49] = (float)array_obj[i+1]->right;
                tmp_obj[i][50] = (float)array_obj[i+1]->left;
                tmp_obj[i][51] = (float)array_obj[i+1]->top;
                tmp_obj[i][52] = (float)array_obj[i+1]->bottom;
                tmp_obj[i][53] = (float)array_obj[i+1]->axis_x1deb;
                tmp_obj[i][54] = (float)array_obj[i+1]->axis_x1end;
                tmp_obj[i][55] = (float)array_obj[i+1]->axis_y1deb;
                tmp_obj[i][56] = (float)array_obj[i+1]->axis_y1end;
                tmp_obj[i][57] = (float)array_obj[i+1]->axis_x2deb;
                tmp_obj[i][58] = (float)array_obj[i+1]->axis_x2end;
                tmp_obj[i][59] = (float)array_obj[i+1]->axis_y2deb;
                tmp_obj[i][60] = (float)array_obj[i+1]->axis_y2end;
          }

        freeimage(o3_img->imagedata);

        o3_img->imagedata = (char *)unload_vector(tmp_obj,0,VFF_TYP_FLOAT,(numreg-1),ALL_OPTS,1,(numreg-1)); 

        *img3 = o3_img;
        } /* viff output */

        return(1);
} 



/* draws a line between (a1,b1) and (a2,b2) in the image data pointed */
/* by the color of this line is color ............................... */
void draw_line(a1,b1,a2,b2,ptr,nc,nr,color)
int a1,b1,a2,b2,color,nr,nc;
int *ptr;
{
float slope;
float deb;

int row, col;

   if (a1 == a2) {
     for (col = (int)b1; col <= (int)b2; col++) {
         row = (int) a1;
         ptr[row * nc + col] = color;
         }
     }
   else {
      slope = (float)(b1 - b2) / (float)(a1 - a2);
      deb   = ((float) a1 * (float) b2 - (float) a2 * (float) b1) /
              (float)(a1 - a2);
  
      /* the equation of the line segment is : y = slope * x + deb . */ 

      /* to draw the line we must take care of the slope .... 2 cases . */
      if ((slope < 1.0) && (slope > -1.0)) {
         for (row = (int) a1; row <= (int) a2; row++) {
            col  = (int) ( (float) ( (float)(row) * slope  + deb ) );
            ptr[row * nc + col] = color;
            }
         }
      else {
         for (col = (int) b1; col <= (int)b2; col++ ) {
            row = (int) ( (float)( ( (float)(col) - deb) / slope) );
            ptr[row * nc + col] = color;
            }
         }
      }
   }

    
/* Computes area, x_cent, y_cent, variance, sigma_x, sigma_y, theta, eccentri- */
/* city, and if required the central and invariant moments.................... */
void comp_shape(meters_surf,obj,numreg,central,invariant,nc,nr)
XVSHAPE **obj;
int  numreg, central, invariant, nc, nr;
double meters_surf;
{
#define POWER2(x) ((double)(x) * (double)(x))
int ind;  
double m00,m10,m01,m11,m02,m20,m21,m12,m03,m30,
       mu11,mu20,mu02,mu12,mu21,mu30,mu03,
       nu11,nu20,nu02,nu12,nu21,nu30,nu03,
       xm,ym,xm2,ym2,
       thet,cos2,sin2,cos_2,sin_2,
       a1,a2,a3,a4,a5,a6,a7,a8,a9,
       puissance;

     /* computes area, x_centroid, y_centroid ........................... */
     for(ind = 0; ind < numreg; ind++) {

        obj[ind]->surf_image = nc * nr;

        m00 = obj[ind]->area = obj[ind]->m00;


        obj[ind]->percent_area = obj[ind]->area*100.0 / obj[ind]->surf_image;

        obj[ind]->real_area = obj[ind]->area * meters_surf;

        if (m00 != 0.0) {
           m01 = obj[ind]->m01;
           m10 = obj[ind]->m10;
           m02 = obj[ind]->m02;
           m20 = obj[ind]->m20;
           m11 = obj[ind]->m11;
           m21 = obj[ind]->m21;
           m12 = obj[ind]->m12;
           m30 = obj[ind]->m30;
           m03 = obj[ind]->m03;
        

           xm = m10 / m00;
           obj[ind]->xcent = xm;

           ym = m01 / m00;
           obj[ind]->ycent = ym;
  
           mu20 = m20 - ( m00 * POWER2( xm ) ) ;

           mu02 = m02 - ( m00 * POWER2( ym ) ) ;

           mu11 = m11 - ( m00 * xm * ym );

           obj[ind]->sigma_x = sqrt( mu20 / m00 ); 
           obj[ind]->sigma_y = sqrt( mu02 / m00 ); 

           if (mu20 == mu02) {
              thet = XV_PI_4; 
              obj[ind]->theta = thet;
              }
           else { 
              thet = ((double) 2.0) * mu11 / (mu20 - mu02) ;
              thet = ((double)0.5)* atan( thet ); 
              obj[ind]->theta = thet;
              }

           cos2  = POWER2((cos( thet )) );
           sin2  = POWER2((sin( thet )) );
           cos_2 = cos( (((double)2.0)*thet) );
           sin_2 = sin( (((double)2.0)*thet) );  

           /* computes the eccentricity of the shape......*/
           a1    = mu02 * cos2 + mu20 * sin2 - mu11 * sin_2;
           a2    = mu02 * sin2 + mu20 * cos2 + mu11 * cos_2; 
           if ((a1 == (double)0.0) || (a2 == (double)0.0)) {
              obj[ind]->eccentricity = (double) -1.0;
              }
           else if ((a1/a2) > 0) obj[ind]->eccentricity = sqrt(a1/a2);
                else obj[ind]->eccentricity = sqrt(-a1/a2);


           /* computes the moments variance on X axis (1st principal axis */
           if (m00 == (double)0.0) {
              obj[ind]->sigma_MX = (double) -1.0;
              obj[ind]->sigma_MY = (double) -1.0;
              }
           else {
              obj[ind]->sigma_MX = sqrt((cos2*mu20+sin2*mu02+sin_2*mu11)/m00);
              obj[ind]->sigma_MY = sqrt((sin2*mu20+cos2*mu02-sin_2*mu11)/m00);
              }

           /* coeff 1.5 to make the axis bigger in the axis image */
           a1 = (double)(1.50) * cos(thet) * obj[ind]->sigma_MX; 
           a2 = (double)(1.50) * sin(thet) * obj[ind]->sigma_MX;
           obj[ind]->axis_x1deb = obj[ind]->xcent - (int)a1;
           if (obj[ind]->axis_x1deb < (double)0.0) 
              obj[ind]->axis_x1deb = (double)2.0;
           obj[ind]->axis_x1end = obj[ind]->xcent + (int)a1;
           if (obj[ind]->axis_x1end > (double)nr) 
              obj[ind]->axis_x1end = (double)(nr-2);
           obj[ind]->axis_y1deb = obj[ind]->ycent - (int)a2;
           if (obj[ind]->axis_y1deb < (double)0.0) 
              obj[ind]->axis_y1deb = (double)2.0;
           obj[ind]->axis_y1end = obj[ind]->ycent + (int)a2;
           if (obj[ind]->axis_y1end > (double)nc) 
              obj[ind]->axis_y1end = (double)(nc -2);

           a1 = -(double)(1.50) * sin(thet)*obj[ind]->sigma_MY;
           a2 =  (double)(1.50) * cos(thet)*obj[ind]->sigma_MY;
           obj[ind]->axis_x2deb = obj[ind]->xcent - (int)a1;
           if (obj[ind]->axis_x1deb < (double)0.0)
              obj[ind]->axis_x1deb = (double)2.0;
           obj[ind]->axis_x2end = obj[ind]->xcent + (int)a1;
           if (obj[ind]->axis_x1end > (double)nr)
              obj[ind]->axis_x1end = (double)(nr-2);
           obj[ind]->axis_y2deb = obj[ind]->ycent - (int)a2;
           if (obj[ind]->axis_y1deb < (double)0.0)
              obj[ind]->axis_y1deb = (double)2.0;
           obj[ind]->axis_y2end = obj[ind]->ycent + (int)a2;
           if (obj[ind]->axis_y1end > (double)nc)
              obj[ind]->axis_y1end = (double)(nc -2);


           if (central == 1 || invariant == 1) {

              /* pre_computation */
              xm2 = POWER2( xm );
              ym2 = POWER2( ym );

              /* The central moments of order 3 are as follows */
              obj[ind]->mu00 = m00;

              obj[ind]->mu10 = obj[ind]->mu01 = 0.0;

              obj[ind]->mu11 = mu11;
              obj[ind]->mu20 = mu20;
              obj[ind]->mu02 = mu02;

              mu30 = m30 - ((double)3.0)*xm*m20 +((double)2.0)*m10*xm2;
              obj[ind]->mu30 = mu30;

              mu03 = m03 - ((double)3.0)*ym*m02 +((double)2.0)*m01*ym2;
              obj[ind]->mu03 = mu03;

              mu21 = m21-m20*ym-((double)2.0)*m11*xm+((double)2.0)*xm2*m01; 
              obj[ind]->mu21 = mu21;
 
              mu12 = m12 -m02*xm - ((double)2.0)*m11*ym + ((double)2.0)*ym2*m10; 
              obj[ind]->mu12 = mu12;
              }

           if (invariant == 1) {
              /*
              * The set of moments that are invariant to translation, 
              * rotation and scale change ( HU [1962] ) are given by 
              * the following relations
              */
              /* Pre-computations ..................................*/
              obj[ind]->nu00 = (double)1.0;
              obj[ind]->nu01 = (double)0.0;
              obj[ind]->nu10 = (double)0.0;

           
              if (m00 == (double)0.0) {

              /*
                 nu02 = obj[ind]->nu02 = (double) -1.0; 
                 nu20 = obj[ind]->nu20 = (double) -1.0; 
                 nu11 = obj[ind]->nu11 = (double) -1.0; 

                 nu12 = obj[ind]->nu12 = (double) -1.0;
                 nu21 = obj[ind]->nu21 = (double) -1.0;
                 nu03 = obj[ind]->nu03 = (double) -1.0; 
                 nu30 = obj[ind]->nu30 = (double) -1.0;

                 a1  =(double) -1.0; 
                 a2  =(double) -1.0; 
                 a3  =(double) -1.0; 
                 a4  =(double) -1.0; 
                 a5  =(double) -1.0; 
                 a6  =(double) -1.0; 
                 a7  =(double) -1.0; 
                 a8  =(double) -1.0; 
                 a9  =(double) -1.0; 
              */

                 obj[ind]->M1  = (double) -1.0; 
                 obj[ind]->M2  = (double) -1.0; 
                 obj[ind]->M3  = (double) -1.0; 
                 obj[ind]->M4  = (double) -1.0; 
                 obj[ind]->M5  = (double) -1.0; 
                 obj[ind]->M6  = (double) -1.0; 
                 obj[ind]->M7  = (double) -1.0; 
                 }

              else {

                 puissance      = m00 * m00;
                 nu02 = obj[ind]->nu02 = mu02 / puissance;
                 nu20 = obj[ind]->nu20 = mu20 / puissance;
                 nu11 = obj[ind]->nu11 = mu11 / puissance;

                 puissance      = pow(m00,(double)2.5);
                 nu12 = obj[ind]->nu12 = mu12 / puissance;
                 nu21 = obj[ind]->nu21 = mu21 / puissance;
                 nu03 = obj[ind]->nu03 = mu03 / puissance;
                 nu30 = obj[ind]->nu30 = mu30 / puissance;

                 a1  = nu20 - nu02;
                 a2  = nu30 + nu12;
                 a3  = nu21 + nu03;
                 a4  = nu30 - 3.0 * nu12;
                 a5  = ((double)3.0) * nu21 - nu03;
                 a6  = POWER2(a2);
                 a7  = POWER2(a3);
                 a8  = a2 * (a6 - ((double)3.0) * a7);
                 a9  = a3 * (((double)3.0) * a6 - a7);

                 obj[ind]->M1  = nu20 + nu02;
                 obj[ind]->M2  = POWER2 (a1) + ((double)4.0) * POWER2(nu11);
                 obj[ind]->M3  = POWER2 (a4) + POWER2(a5);
                 obj[ind]->M4  = a6 +a7;
                 obj[ind]->M5  = a4 *a8  +a5 *a9;
                 obj[ind]->M6  = a1 *(a6 -a7) +((double)4.0) *nu11 *a2 *a3;
                 obj[ind]->M7 =  ((double)3.0 *nu21 - nu30) *a8 -a4 *a9;
                 }
              }
           }
        }
     }

/* Print information in INPUT FILE *********************************/
int
printout(typ_im,obj,numreg,standard,central,invariant,input_file, printdev)
XVSHAPE **obj;
int numreg,typ_im;
char    *input_file;
FILE    *printdev;
{
    int   i,new_i,i_deb;

    (void) fprintf(printdev,"\nShape Information For File Name: %s\n",
                   input_file);
    new_i = 0;

    if (typ_im == 0 ) i_deb=0; /* if grey level input image print region 0 */
    else i_deb = 1;            /* if labelled image doesn't print region 0 */
    for (i = i_deb ; i < numreg ; i++) { 
     if (obj[i]->area != 0.0) {
      (void) fprintf(printdev,"\"Object number  %d \"\n",new_i);        
      (void) fprintf(printdev,
      "Entire Image    :  %.0f pixels \n",obj[i]->surf_image); 
      (void) fprintf(printdev,
      "Object          :  %.0f pixels \n",obj[i]->area); 
      (void) fprintf(printdev,
      "Ratio Obj/Image :  %.2f %%\n",obj[i]->percent_area); 
      (void) fprintf(printdev,
      "Area            :  %e square meters \n",obj[i]->area); 
      (void) fprintf(printdev,"Xcenter         :  %e\n",obj[i]->xcent); 
      (void) fprintf(printdev,"Ycenter         :  %e\n",obj[i]->ycent); 
      (void) fprintf(printdev,"Sigma_x         :  %e\n",obj[i]->sigma_x); 
      (void) fprintf(printdev,"Sigma_y         :  %e\n",obj[i]->sigma_y); 
      (void) fprintf(printdev,"Theta           :  %e\n",obj[i]->theta); 
      (void) fprintf(printdev,"Eccentricity    :  %e\n",obj[i]->eccentricity); 
      (void) fprintf(printdev,"Uppermost point :  %d\n",obj[i]->top); 
      (void) fprintf(printdev,"Leftmost point  :  %d\n",obj[i]->left); 
      (void) fprintf(printdev,"Lowest point    :  %d\n",obj[i]->bottom); 
      (void) fprintf(printdev,"Rightmost point :  %d\n\n",obj[i]->right); 

      if (standard == 1) {
         (void) fprintf(printdev,"Standard Moments\n");
         (void) fprintf(printdev,"m00             :  %e\n",obj[i]->m00); 
         (void) fprintf(printdev,"m10             :  %e\n",obj[i]->m10); 
         (void) fprintf(printdev,"m01             :  %e\n",obj[i]->m01); 
         (void) fprintf(printdev,"m11             :  %e\n",obj[i]->m11); 
         (void) fprintf(printdev,"m20             :  %e\n",obj[i]->m20); 
         (void) fprintf(printdev,"m02             :  %e\n",obj[i]->m02); 
         (void) fprintf(printdev,"m30             :  %e\n",obj[i]->m30); 
         (void) fprintf(printdev,"m03             :  %e\n",obj[i]->m03); 
         (void) fprintf(printdev,"m12             :  %e\n",obj[i]->m12); 
         (void) fprintf(printdev,"m21             :  %e\n\n",obj[i]->m21); 
         }
      if (central == 1) {
         (void) fprintf(printdev,"Central  Moments\n");
         (void) fprintf(printdev,"mu00            :  %e\n",obj[i]->mu00); 
         (void) fprintf(printdev,"mu11            :  %e\n",obj[i]->mu11); 
         (void) fprintf(printdev,"mu20            :  %e\n",obj[i]->mu20); 
         (void) fprintf(printdev,"mu02            :  %e\n",obj[i]->mu02); 
         (void) fprintf(printdev,"mu30            :  %e\n",obj[i]->mu30); 
         (void) fprintf(printdev,"mu03            :  %e\n",obj[i]->mu03); 
         (void) fprintf(printdev,"mu12            :  %e\n",obj[i]->mu12); 
         (void) fprintf(printdev,"mu21            :  %e\n\n",obj[i]->mu21); 
         }
      if (invariant == 1) {
         (void) fprintf(printdev,"Normalized Central  Moments\n");
         (void) fprintf(printdev,"nu00            :  %e\n",obj[i]->nu00); 
         (void) fprintf(printdev,"nu11            :  %e\n",obj[i]->nu11); 
         (void) fprintf(printdev,"nu20            :  %e\n",obj[i]->nu20); 
         (void) fprintf(printdev,"nu02            :  %e\n",obj[i]->nu02); 
         (void) fprintf(printdev,"nu30            :  %e\n",obj[i]->nu30); 
         (void) fprintf(printdev,"nu03            :  %e\n",obj[i]->nu03); 
         (void) fprintf(printdev,"nu12            :  %e\n",obj[i]->nu12); 
         (void) fprintf(printdev,"nu21            :  %e\n\n",obj[i]->nu21); 

         (void) fprintf(printdev,"Invariant Moments\n");
         (void) fprintf(printdev,"M1              :  %e\n",obj[i]->M1);
         (void) fprintf(printdev,"M2              :  %e\n",obj[i]->M2);
         (void) fprintf(printdev,"M3              :  %e\n",obj[i]->M3);
         (void) fprintf(printdev,"M4              :  %e\n",obj[i]->M4);
         (void) fprintf(printdev,"M5              :  %e\n",obj[i]->M5);
         (void) fprintf(printdev,"M6              :  %e\n",obj[i]->M6);
         (void) fprintf(printdev,"M7              :  %e\n\n",obj[i]->M7);;
         }

      (void) fprintf(printdev,".....................................\n\n");
      }
     new_i++;
     }
   }

/* -library_code_end */
