1 #include "mrilib.h"
2
3 /** not 7D SAFE **/
4
5 /*--------------------------------------------------------------------
6 12 Nov 2001:
7 Get the center of mass of a 2D image;
8 store it in the user-supplied locations *xcm, *ycm.
9 Note that *xcm and *ycm are the INDEX center of mass coordinates.
10 ----------------------------------------------------------------------*/
11
mri_get_cmass_2D(MRI_IMAGE * im,float * xcm,float * ycm)12 void mri_get_cmass_2D( MRI_IMAGE *im , float *xcm , float *ycm )
13 {
14 int ii,jj,joff , nx,ny ;
15 double xx , yy , sum , val ; float *far ;
16 MRI_IMAGE *flim ;
17
18 ENTRY("mri_get_cmass_2D") ;
19
20 if( im == NULL || xcm == NULL || ycm == NULL ) EXRETURN ;
21
22 if( im->kind != MRI_float ) flim = mri_to_float( im ) ;
23 else flim = im ;
24
25 far = MRI_FLOAT_PTR(flim) ;
26 nx = im->nx ; ny = im->ny ;
27
28 sum = xx = yy = 0.0 ;
29 for( jj=0 ; jj < ny ; jj++ ){
30 joff = jj * nx ;
31 for( ii=0 ; ii < nx ; ii++ ){
32 val = fabs(far[ii+joff]) ;
33 sum += val ;
34 xx += val * ii ;
35 yy += val * jj ;
36 }
37 }
38 if( flim != im ) mri_free(flim) ;
39
40 if( sum > 0.0 ){ xx /= sum ; yy /= sum ; }
41 else { xx = 0.5*(nx-1); yy=0.5*(ny-1); }
42
43 *xcm = (float)xx ; *ycm = (float)yy ; EXRETURN ;
44 }
45
46 /*--------------------------------------------------------------------*/
47 /*! Get the center of mass of a 3D image (in index coordinates). */
48
mri_get_cmass_3D(MRI_IMAGE * im,float * xcm,float * ycm,float * zcm)49 void mri_get_cmass_3D( MRI_IMAGE *im, float *xcm, float *ycm, float *zcm )
50 {
51 int ii,jj,kk,koff,joff , nx,ny,nz,nxy ;
52 double xx,yy,zz , sum , val ; float *far ;
53 MRI_IMAGE *flim ;
54
55 ENTRY("mri_get_cmass_3D") ;
56
57 if( im == NULL || xcm == NULL || ycm == NULL || zcm == NULL ) EXRETURN ;
58
59 if( im->kind != MRI_float ) flim = mri_to_float( im ) ;
60 else flim = im ;
61
62 far = MRI_FLOAT_PTR(flim) ;
63 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
64
65 sum = xx = yy = zz = 0.0 ;
66 for( kk=0 ; kk < nz ; kk++ ){
67 koff = kk * nxy ;
68 for( jj=0 ; jj < ny ; jj++ ){
69 joff = koff + jj * nx ;
70 for( ii=0 ; ii < nx ; ii++ ){
71 val = fabs(far[ii+joff]) ;
72 sum += val ;
73 xx += val * ii ;
74 yy += val * jj ;
75 zz += val * kk ;
76 }
77 }
78 }
79 if( flim != im ) mri_free(flim) ;
80
81 if( sum > 0.0 ){ xx /= sum ; yy /= sum ; zz /= sum ; }
82 else { xx = 0.5*(nx-1); yy=0.5*(ny-1); zz=0.5*(nz-1); }
83
84 *xcm = (float)xx ; *ycm = (float)yy ; *zcm = (float)zz ; EXRETURN ;
85 }
86