1 /*---------------------------------------------------------------------
2 01 Feb 2001: compute the center of a dataset (DICOM coord order)
3 Adapted from part of thd_info.c -- RWCox
4 -----------------------------------------------------------------------*/
5
6 #include "mrilib.h"
7
THD_dataset_center(THD_3dim_dataset * dset)8 THD_fvec3 THD_dataset_center( THD_3dim_dataset * dset )
9 {
10 THD_dataxes * daxes ;
11 THD_fvec3 fv1 , fv2 ;
12
13 ENTRY("THD_dataset_center") ;
14
15 if( !ISVALID_DSET(dset) ){ LOAD_FVEC3(fv1,0,0,0); RETURN(fv1); }
16
17 daxes = dset->daxes ;
18
19 LOAD_FVEC3(fv1 , daxes->xxorg , daxes->yyorg , daxes->zzorg) ;
20 fv1 = THD_3dmm_to_dicomm( dset , fv1 ) ;
21
22 LOAD_FVEC3(fv2 , daxes->xxorg + (daxes->nxx-1)*daxes->xxdel ,
23 daxes->yyorg + (daxes->nyy-1)*daxes->yydel ,
24 daxes->zzorg + (daxes->nzz-1)*daxes->zzdel ) ;
25 fv2 = THD_3dmm_to_dicomm( dset , fv2 ) ;
26
27 fv1.xyz[0] = 0.5f * (fv1.xyz[0]+fv2.xyz[0]) ;
28 fv1.xyz[1] = 0.5f * (fv1.xyz[1]+fv2.xyz[1]) ;
29 fv1.xyz[2] = 0.5f * (fv1.xyz[2]+fv2.xyz[2]) ;
30
31 RETURN(fv1) ;
32 }
33
34 /*-------------------------------------------------------------------------*/
35 /*! Get the center of mass of this volume, in DICOM coords, taken from 3dCM.
36
37 // [PT, 20 Dec, 2016] Updated to output either: ijk, or the (default)
38 // xyz in DICOM coords. New input variable 'cmode' controls this
39 // (default, cmode=0).
40
41 ---------------------------------------------------------------------------*/
42
THD_cmass(THD_3dim_dataset * xset,int iv,byte * mmm,int cmode)43 THD_fvec3 THD_cmass( THD_3dim_dataset *xset , int iv , byte *mmm,
44 int cmode)
45 {
46 THD_fvec3 cmv ;
47 THD_ivec3 cmvi ;
48 MRI_IMAGE *im ;
49 float *far , icm,jcm,kcm ;
50 int ii , nvox ;
51
52 LOAD_FVEC3(cmv,0,0,0) ;
53
54 nvox = DSET_NVOX(xset) ;
55 im = mri_to_float( DSET_BRICK(xset,iv) ) ;
56 if( im == NULL ) return cmv ;
57 far = MRI_FLOAT_PTR(im) ; if( far == NULL ) return cmv ;
58
59 if( mmm != NULL ){
60 for( ii=0 ; ii < nvox ; ii++ )
61 /* if mask is NOT set, clear it 2 Jan 2013 [rickr/dglen] */
62 if( ! mmm[ii] ) far[ii] = 0.0 ;
63 }
64
65 mri_get_cmass_3D( im , &icm,&jcm,&kcm ) ; mri_free(im) ;
66
67 LOAD_FVEC3(cmv,icm,jcm,kcm) ;
68 cmv = THD_3dfind_to_3dmm( xset , cmv ) ;
69
70 // [PT, Dec, 2016]: a bit convoluted, but return ijk (as floats) if
71 // the mode suits it.
72 if( cmode == 1 ) { // return int ijk in dset-orientation
73 cmvi = THD_3dmm_to_3dind( xset , cmv );
74 LOAD_FVEC3(cmv,cmvi.ijk[0], cmvi.ijk[1], cmvi.ijk[2]);
75 return cmv;
76 }
77
78 cmv = THD_3dmm_to_dicomm( xset , cmv ) ;
79 return cmv ;
80 }
81
82 /*-------------------------------------------------------------------------*/
83 /*! Get the center of mass of integer labeled ROIs in a sub-brick. Returns
84 a float vector N_rois XYZ triplets.
85 ---------------------------------------------------------------------------*/
86
THD_roi_cmass(THD_3dim_dataset * xset,int iv,int * rois,int N_rois,int cmode)87 float *THD_roi_cmass(THD_3dim_dataset *xset , int iv , int *rois,
88 int N_rois, int cmode)
89 {
90 float *xyz=NULL, roi;
91 THD_fvec3 cmr ;
92 int ir;
93 byte *mmm;
94
95 ENTRY("THD_roi_cmass");
96
97 if (!xset || !rois || N_rois < 1) RETURN(NULL);
98
99 xyz = (float *)calloc(N_rois*3, sizeof(float));
100 for (ir = 0; ir < N_rois; ++ir) {
101 roi = rois[ir];
102 mmm = THD_makemask( xset, iv , roi , roi );
103 cmr = THD_cmass( xset, iv , mmm, cmode);
104 free(mmm); mmm = NULL;
105 xyz[3*ir] = cmr.xyz[0];
106 xyz[3*ir+1] = cmr.xyz[1];
107 xyz[3*ir+2] = cmr.xyz[2];
108 }
109
110 RETURN(xyz);
111 }
112
113 /* find internal center "Icent" in volume
114 * for some shapes, cmass can be outside the region
115 * this method forces the result to be in a voxel
116 * Parameters: xset=dataset, iv = subbrick, mmm=mask array
117 * cmode = 0/1 return xyz position or ijk index
118 * cmxyz = cm center of mass in xyz */
THD_Icent(THD_3dim_dataset * xset,int iv,byte * mmm,int cmode,THD_fvec3 cmxyz)119 THD_fvec3 THD_Icent( THD_3dim_dataset *xset , int iv , byte *mmm,
120 int cmode, THD_fvec3 cmxyz)
121 {
122 THD_fvec3 cmv ;
123 THD_ivec3 cmvi ;
124 MRI_IMAGE *im ;
125 float *far, xcm, ycm, zcm ;
126 int ii, jj, kk, nvox ;
127 THD_dataxes *daxes ;
128 int koff, joff, nx,ny,nz,nxy ;
129 double xx,yy,zz , vsum , val ;
130 double wbest;
131 float xmi,ymi,zmi;
132
133 xcm = cmxyz.xyz[0];
134 ycm = cmxyz.xyz[1];
135 zcm = cmxyz.xyz[2];
136
137 LOAD_FVEC3(cmv,0,0,0) ;
138
139 nvox = DSET_NVOX(xset) ;
140 im = mri_to_float( DSET_BRICK(xset,iv) ) ;
141 if( im == NULL ) return cmv ;
142 far = MRI_FLOAT_PTR(im) ; if( far == NULL ) return cmv ;
143
144 if( mmm != NULL ){
145 for( ii=0 ; ii < nvox ; ii++ )
146 /* if mask is NOT set, clear it 2 Jan 2013 [rickr/dglen] */
147 if( ! mmm[ii] ) far[ii] = 0.0 ;
148 }
149
150
151 daxes = CURRENT_DAXES(xset) ;
152
153 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
154 /* some default distance to beat and a default center coordinate */
155 wbest = 1.e+37f ; xmi=ymi=zmi = 0.0f ;
156
157 for( kk=0 ; kk < nz ; kk++ ){
158 koff = kk * nxy ;
159 zz = daxes->zzorg + kk * daxes->zzdel;
160 for( jj=0 ; jj < ny ; jj++ ){
161 joff = koff + jj * nx ;
162 yy = daxes->yyorg + jj * daxes->yydel;
163 for( ii=0 ; ii < nx ; ii++ ){
164 if(far[ii+joff]!=0.0) {
165 xx = daxes->xxorg + ii * daxes->xxdel;
166 vsum = SQR(xx-xcm) + SQR(yy-ycm) + SQR(zz-zcm); /* use dist^2 here for a bit of speed*/
167 if( vsum < wbest ){ xmi = xx; ymi = yy; zmi = zz; wbest=vsum; }
168 }
169 }
170 }
171 }
172
173 mri_free(im) ;
174
175 LOAD_FVEC3(cmv,xmi,ymi,zmi) ;
176
177 // [PT, Dec, 2016]: a bit convoluted, but return ijk (as floats) if
178 // the mode suits it.
179 if( cmode == 1 ) { // return int ijk in dset-orientation
180 cmvi = THD_3dmm_to_3dind( xset , cmv );
181 LOAD_FVEC3(cmv,cmvi.ijk[0], cmvi.ijk[1], cmvi.ijk[2]);
182 return cmv;
183 }
184
185 cmv = THD_3dmm_to_dicomm( xset , cmv ) ;
186 return cmv ;
187 }
188
189 /* find distance center "Dcent" in volume
190 * for some shapes, cmass can be outside the region
191 * like Icent this method forces the result to be in a voxel
192 ** but computes the center as the voxel location with the smallest average
193 * distance to all the others
194 * Parameters: xset=dataset, iv = subbrick, mmm=mask array
195 * cmode = 0/1 return xyz position or ijk index
196 * cmxyz = cm center of mass in xyz */
THD_Dcent(THD_3dim_dataset * xset,int iv,byte * mmm,int cmode,THD_fvec3 cmxyz)197 THD_fvec3 THD_Dcent( THD_3dim_dataset *xset , int iv , byte *mmm,
198 int cmode, THD_fvec3 cmxyz)
199 {
200 THD_fvec3 cmv ;
201 THD_ivec3 cmvi ;
202 MRI_IMAGE *im ;
203 float *far ;
204 int ii, jj, kk, nvox ;
205 THD_dataxes *daxes ;
206 int koff, joff, nx,ny,nz,nxy ;
207 double xx,yy,zz , vsum , val ;
208 double wbest;
209 float xmi,ymi,zmi, xorg,yorg, zorg, xdel, ydel, zdel;
210
211 LOAD_FVEC3(cmv,0,0,0) ;
212
213 nvox = DSET_NVOX(xset) ;
214 im = mri_to_float( DSET_BRICK(xset,iv) ) ;
215 if( im == NULL ) return cmv ;
216 far = MRI_FLOAT_PTR(im) ; if( far == NULL ) return cmv ;
217
218 if( mmm != NULL ){
219 for( ii=0 ; ii < nvox ; ii++ )
220 /* if mask is NOT set, clear it 2 Jan 2013 [rickr/dglen] */
221 if( ! mmm[ii] ) far[ii] = 0.0 ;
222 }
223
224
225 daxes = CURRENT_DAXES(xset) ;
226
227 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
228 /* some default distance to beat and a default center coordinate */
229 wbest = 1.e+37f ; xmi=ymi=zmi = 0.0f ;
230 xorg = daxes->xxorg;
231 yorg = daxes->yyorg;
232 zorg = daxes->zzorg;
233 xdel = daxes->xxdel;
234 ydel = daxes->yydel;
235 zdel = daxes->zzdel;
236
237 for( kk=0 ; kk < nz ; kk++ ){
238 koff = kk * nxy ;
239 for( jj=0 ; jj < ny ; jj++ ){
240 joff = koff + jj * nx ;
241 for( ii=0 ; ii < nx ; ii++ ){
242 if(far[ii+joff]!=0.0) {
243 zz = zorg + kk * zdel;
244 yy = yorg + jj * ydel;
245 xx = xorg + ii * xdel;
246 /* find average distance of all voxels in the mask to this voxel */
247 vsum = THD_xyz_distance(xset, im, xx,yy,zz);
248 if( vsum < wbest ){ xmi = xx; ymi = yy; zmi = zz; wbest=vsum; }
249 }
250 }
251 }
252 }
253
254 mri_free(im) ;
255
256 LOAD_FVEC3(cmv,xmi,ymi,zmi) ;
257
258 // [PT, Dec, 2016]: a bit convoluted, but return ijk (as floats) if
259 // the mode suits it.
260 if( cmode == 1 ) { // return int ijk in dset-orientation
261 cmvi = THD_3dmm_to_3dind( xset , cmv );
262 LOAD_FVEC3(cmv,cmvi.ijk[0], cmvi.ijk[1], cmvi.ijk[2]);
263 return cmv;
264 }
265
266 cmv = THD_3dmm_to_dicomm( xset , cmv ) ;
267 return cmv ;
268 }
269
270 /* compute average distance in mask to xyz location */
THD_xyz_distance(THD_3dim_dataset * xset,MRI_IMAGE * im,double xcm,double ycm,double zcm)271 double THD_xyz_distance( THD_3dim_dataset *xset , MRI_IMAGE *im ,
272 double xcm, double ycm, double zcm)
273 {
274 float *far;
275 int ii, jj, kk, nvox, nmask ;
276 THD_dataxes *daxes ;
277 int koff, joff, nx,ny,nz,nxy ;
278 double xx,yy,zz , vsum , bad_vsum, val ;
279
280 nmask = 0; bad_vsum = 1.e+37f;
281 vsum = 0.0;
282
283 nvox = DSET_NVOX(xset) ;
284 far = MRI_FLOAT_PTR(im) ; if( far == NULL ) return bad_vsum ;
285
286 daxes = CURRENT_DAXES(xset) ;
287
288 nx = im->nx ; ny = im->ny ; nz = im->nz ; nxy = nx*ny ;
289
290 for( kk=0 ; kk < nz ; kk++ ){
291 koff = kk * nxy ;
292 for( jj=0 ; jj < ny ; jj++ ){
293 joff = koff + jj * nx ;
294 for( ii=0 ; ii < nx ; ii++ ){
295 if(far[ii+joff]!=0.0) {
296 zz = daxes->zzorg + kk * daxes->zzdel;
297 yy = daxes->yyorg + jj * daxes->yydel;
298 xx = daxes->xxorg + ii * daxes->xxdel;
299 vsum += SQR(xx-xcm) + SQR(yy-ycm) + SQR(zz-zcm) ; /* use dist^2 here for a bit of speed*/
300 nmask++;
301 }
302 }
303 }
304 }
305
306 if(nmask != 0){
307 vsum = vsum / nmask;
308 return(vsum);
309 }
310 else return(bad_vsum);
311 }
312