1 #include "mrilib.h"
2
3 /****************************************************************************
4 Functions for dealing with data extracted from a neighborhood of a voxel
5 in a dataset.
6 -- 19 Aug 2005 - RWCox
7 *****************************************************************************/
8
9 #undef INMASK
10 #define INMASK(i) ( mask == NULL || mask[i] != 0 )
11
12 /*---------------------------------------------------------------------------*/
13 /*! Extract a nbhd from a dataset sub-brick.
14 The 3-index of a voxel from the row is in (xx,yy,zz). The result
15 is a 1D MRI_IMAGE *im, with im->nx = number of points extracted.
16 If mask!=NULL, then it is a mask of allowed points.
17 -----------------------------------------------------------------------------*/
18
THD_get_dset_nbhd(THD_3dim_dataset * dset,int ival,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)19 MRI_IMAGE * THD_get_dset_nbhd( THD_3dim_dataset *dset, int ival, byte *mask,
20 int xx, int yy, int zz, MCW_cluster *nbhd )
21 {
22 int kind , nx,ny,nz,nxy , npt , nout , aa,bb,cc,kk,ii ;
23 void *brick ;
24 MRI_IMAGE *im ;
25 float fac ;
26
27 ENTRY("THD_get_dset_nbhd") ;
28
29 if( dset == NULL || nbhd == NULL || nbhd->num_pt <= 0 ) RETURN(NULL) ;
30 if( ival < 0 || ival >= DSET_NVALS(dset) ) RETURN(NULL) ;
31
32 im = mri_get_nbhd( DSET_BRICK(dset,ival) , mask , xx,yy,zz , nbhd ) ;
33 if( im == NULL ) RETURN(NULL) ;
34
35 fac = DSET_BRICK_FACTOR(dset,ival) ;
36 if( fac != 0.0f && fac != 1.0f ){
37 MRI_IMAGE *qim = mri_scale_to_float( fac , im ) ;
38 mri_free(im) ; im = qim ;
39 }
40
41 RETURN(im) ;
42 }
43
44 /*-------------------------------------------------------------------------*/
45
THD_get_dset_indexed_nbhd(THD_3dim_dataset * dset,int ival,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)46 MRI_IMARR * THD_get_dset_indexed_nbhd( THD_3dim_dataset *dset,
47 int ival, byte *mask,
48 int xx, int yy, int zz,
49 MCW_cluster *nbhd )
50 {
51 int kind , nx,ny,nz,nxy , npt , nout , aa,bb,cc,kk,ii ;
52 void *brick ;
53 MRI_IMAGE *im , *qm ;
54 MRI_IMARR *imar ;
55 float fac ;
56
57 ENTRY("THD_get_dset_indexed_nbhd") ;
58
59 if( dset == NULL || nbhd == NULL || nbhd->num_pt <= 0 ) RETURN(NULL) ;
60 if( ival < 0 || ival >= DSET_NVALS(dset) ) RETURN(NULL) ;
61
62 imar = mri_get_indexed_nbhd( DSET_BRICK(dset,ival), mask, xx,yy,zz, nbhd );
63 if( imar == NULL ) RETURN(NULL) ;
64
65 fac = DSET_BRICK_FACTOR(dset,ival) ;
66 if( fac == 0.0f || fac == 1.0f ) RETURN(imar) ;
67
68 im = IMARR_SUBIM(imar,0) ;
69 qm = mri_scale_to_float( fac , im ) ; mri_free(im) ;
70 IMARR_SUBIM(imar,0) = qm ;
71 RETURN(imar) ;
72 }
73
74 /*-------------------------------------------------------------------------*/
75
76 static byte SearchAboutMaskedVoxel = 0;
77
SetSearchAboutMaskedVoxel(int v)78 void SetSearchAboutMaskedVoxel(int v)
79 {
80 SearchAboutMaskedVoxel = v;
81 }
82
83 /*---------------------------------------------------------------------------*/
84 /*!
85 This function loads the 1D indices of voxels that are in the nbhd
86 Indices are loaded into array nind. The function returns the number
87 of neighboring voxels in the mask. -1 in case of error */
mri_load_nbhd_indices(int nx,int ny,int nz,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,int * nind)88 int mri_load_nbhd_indices ( int nx, int ny, int nz , byte *mask ,
89 int xx, int yy, int zz, MCW_cluster *nbhd,
90 int *nind)
91 {
92 int nxy, nxyz , npt , nout , aa,bb,cc,kk,ii ;
93
94 ENTRY("mri_load_nbhd_indices") ;
95
96 if( nbhd == NULL || nind == NULL ) RETURN(-1) ;
97
98 nxy = nx*ny ;
99 nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
100
101 kk = xx + yy*nx + zz*nxy ;
102 if (SearchAboutMaskedVoxel) {
103 if( npt == 0 || kk < 0 || kk >= nxyz ) RETURN(0) ;
104 } else {
105 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) RETURN(0) ;
106 }
107
108 for( ii=0 ; ii < npt ; ii++ ){
109 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
110 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
111 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
112 kk = aa + bb*nx + cc*nxy ;
113 if( INMASK(kk) ) nind[nout++] = kk ;
114 }
115
116 RETURN(nout);
117 }
118
119 /*---------------------------------------------------------------------------*/
120
mri_get_nbhd(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)121 MRI_IMAGE * mri_get_nbhd( MRI_IMAGE *inim , byte *mask ,
122 int xx, int yy, int zz, MCW_cluster *nbhd )
123 {
124 int kind , nx,ny,nz,nxy,nxyz , npt , nout , aa,bb,cc,kk,ii ;
125 MRI_IMAGE *im ;
126 void *brick ;
127
128 ENTRY("mri_get_nbhd") ;
129
130 if( inim == NULL || nbhd == NULL ) RETURN(NULL) ;
131
132 nx = inim->nx ;
133 ny = inim->ny ; nxy = nx*ny ;
134 nz = inim->nz ; nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
135
136 kk = xx + yy*nx + zz*nxy ;
137 if (SearchAboutMaskedVoxel) {
138 if( npt == 0 || kk < 0 || kk >= nxyz ) RETURN(NULL) ;
139 } else {
140 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) RETURN(NULL) ;
141 }
142 kind = inim->kind ;
143 brick = mri_data_pointer(inim) ; if( brick == NULL ) RETURN(NULL) ;
144 im = mri_new( npt , 1 , kind ) ;
145
146 /*-- extract data, based on kind of data in sub-brick --*/
147
148 switch( kind ){
149
150 default: break ; /* bad bad bad */
151
152 case MRI_short:{
153 short *rr = MRI_SHORT_PTR(im) , *vv = (short *)brick ;
154 for( ii=0 ; ii < npt ; ii++ ){
155 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
156 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
157 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
158 kk = aa + bb*nx + cc*nxy ;
159 if( INMASK(kk) ) rr[nout++] = vv[kk] ;
160 }
161 }
162 break ;
163
164 case MRI_byte:{
165 byte *rr = MRI_BYTE_PTR(im) , *vv = (byte *)brick ;
166 for( ii=0 ; ii < npt ; ii++ ){
167 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
168 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
169 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
170 kk = aa + bb*nx + cc*nxy ;
171 if( INMASK(kk) ) rr[nout++] = vv[kk] ;
172 }
173 }
174 break ;
175
176 case MRI_float:{
177 float *rr = MRI_FLOAT_PTR(im) , *vv = (float *)brick ;
178 for( ii=0 ; ii < npt ; ii++ ){
179 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
180 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
181 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
182 kk = aa + bb*nx + cc*nxy ;
183 if( INMASK(kk) ) rr[nout++] = vv[kk] ;
184 }
185 thd_floatscan( nout , rr ) ; /* 10 Jun 2021 */
186 }
187 break ;
188
189 case MRI_complex:{
190 complex *rr = MRI_COMPLEX_PTR(im) , *vv = (complex *)brick ;
191 for( ii=0 ; ii < npt ; ii++ ){
192 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
193 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
194 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
195 kk = aa + bb*nx + cc*nxy ;
196 if( INMASK(kk) ) rr[nout++] = vv[kk] ;
197 }
198 thd_complexscan( nout , rr ) ; /* 10 Jun 2021 */
199 }
200 break ;
201
202 case MRI_rgb:{
203 byte *rr = MRI_BYTE_PTR(im) , *vv = (byte *)brick ;
204 for( ii=0 ; ii < npt ; ii++ ){
205 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
206 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
207 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
208 kk = aa + bb*nx + cc*nxy ;
209 if( INMASK(kk) ){
210 rr[3*nout ] = vv[3*kk ] ;
211 rr[3*nout+1] = vv[3*kk+1] ;
212 rr[3*nout+2] = vv[3*kk+2] ; nout++ ;
213 }
214 }
215 }
216 break ;
217 }
218
219 if( nout == 0 ){ mri_free(im) ; RETURN(NULL) ; }
220
221 im->nx = im->nvox = nout ; RETURN(im) ;
222 }
223
224 /*---------------------------------------------------------------------------*/
225
mri_get_nbhd_array_float(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,float * nar)226 static int mri_get_nbhd_array_float( MRI_IMAGE *inim , byte *mask ,
227 int xx, int yy, int zz, MCW_cluster *nbhd , float *nar )
228 {
229 int kind , nx,ny,nz,nxy,nxyz , npt , nout , aa,bb,cc,kk,ii ;
230 MRI_IMAGE *im ;
231 float *brick ;
232
233 if( inim == NULL || nbhd == NULL || nar == NULL ) return(0) ;
234
235 nx = inim->nx ;
236 ny = inim->ny ; nxy = nx*ny ;
237 nz = inim->nz ; nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
238
239 kk = xx + yy*nx + zz*nxy ; /* voxel index in inim array */
240 if (SearchAboutMaskedVoxel) {
241 if( npt == 0 || kk < 0 || kk >= nxyz ) return(0) ;
242 } else {
243 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) return(0) ;
244 }
245 kind = inim->kind ; if( kind != MRI_float ) return(0) ;
246 brick = MRI_FLOAT_PTR(inim) ; if( brick == NULL ) return(0) ;
247
248 /*-- extract data, based on kind of data in sub-brick --*/
249
250 for( ii=0 ; ii < npt ; ii++ ){
251 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
252 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
253 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
254 kk = aa + bb*nx + cc*nxy ;
255 if( INMASK(kk) ) nar[nout++] = brick[kk] ;
256 }
257
258 thd_floatscan( nout , nar ) ; /* 10 Jun 2021 */
259
260 return(nout) ;
261 }
262
263 /*---------------------------------------------------------------------------*/
264
mri_get_nbhd_array_short(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,short * nar)265 static int mri_get_nbhd_array_short( MRI_IMAGE *inim , byte *mask ,
266 int xx, int yy, int zz, MCW_cluster *nbhd , short *nar )
267 {
268 int kind , nx,ny,nz,nxy,nxyz , npt , nout , aa,bb,cc,kk,ii ;
269 MRI_IMAGE *im ;
270 short *brick ;
271
272 if( inim == NULL || nbhd == NULL || nar == NULL ) return(0) ;
273
274 nx = inim->nx ;
275 ny = inim->ny ; nxy = nx*ny ;
276 nz = inim->nz ; nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
277
278 kk = xx + yy*nx + zz*nxy ; /* voxel index in inim array */
279 if (SearchAboutMaskedVoxel) {
280 if( npt == 0 || kk < 0 || kk >= nxyz ) return(0) ;
281 } else {
282 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) return(0) ;
283 }
284 kind = inim->kind ; if( kind != MRI_short ) return(0) ;
285 brick = MRI_SHORT_PTR(inim) ; if( brick == NULL ) return(0) ;
286
287 /*-- extract data, based on kind of data in sub-brick --*/
288
289 for( ii=0 ; ii < npt ; ii++ ){
290 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
291 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
292 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
293 kk = aa + bb*nx + cc*nxy ;
294 if( INMASK(kk) ) nar[nout++] = brick[kk] ;
295 }
296
297 return(nout) ;
298 }
299
300 /*---------------------------------------------------------------------------*/
301
mri_get_nbhd_array_byte(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,byte * nar)302 static int mri_get_nbhd_array_byte( MRI_IMAGE *inim , byte *mask ,
303 int xx, int yy, int zz, MCW_cluster *nbhd , byte *nar )
304 {
305 int kind , nx,ny,nz,nxy,nxyz , npt , nout , aa,bb,cc,kk,ii ;
306 MRI_IMAGE *im ;
307 byte *brick ;
308
309 if( inim == NULL || nbhd == NULL || nar == NULL ) return(0) ;
310
311 nx = inim->nx ;
312 ny = inim->ny ; nxy = nx*ny ;
313 nz = inim->nz ; nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
314
315 kk = xx + yy*nx + zz*nxy ; /* voxel index in inim array */
316 if (SearchAboutMaskedVoxel) {
317 if( npt == 0 || kk < 0 || kk >= nxyz ) return(0) ;
318 } else {
319 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) return(0) ;
320 }
321 kind = inim->kind ; if( kind != MRI_byte ) return(0) ;
322 brick = MRI_BYTE_PTR(inim) ; if( brick == NULL ) return(0) ;
323
324 /*-- extract data, based on kind of data in sub-brick --*/
325
326 for( ii=0 ; ii < npt ; ii++ ){
327 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
328 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
329 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
330 kk = aa + bb*nx + cc*nxy ;
331 if( INMASK(kk) ) nar[nout++] = brick[kk] ;
332 }
333
334 return(nout) ;
335 }
336
337 /*---------------------------------------------------------------------------*/
338
mri_get_nbhd_array(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd,void * nar)339 int mri_get_nbhd_array( MRI_IMAGE *inim , byte *mask ,
340 int xx, int yy, int zz, MCW_cluster *nbhd , void *nar )
341 {
342 if( inim == NULL || nbhd == NULL || nar == NULL ) return 0 ;
343
344 switch( inim->kind ){
345 default: break ;
346 case MRI_float:
347 return mri_get_nbhd_array_float(inim,mask,xx,yy,zz,nbhd,(float *)nar) ;
348 case MRI_short:
349 return mri_get_nbhd_array_short(inim,mask,xx,yy,zz,nbhd,(short *)nar) ;
350 case MRI_byte:
351 return mri_get_nbhd_array_byte (inim,mask,xx,yy,zz,nbhd,(byte *) nar) ;
352 }
353 return 0 ;
354 }
355
356 /*---------------------------------------------------------------------------*/
357
mri_get_indexed_nbhd(MRI_IMAGE * inim,byte * mask,int xx,int yy,int zz,MCW_cluster * nbhd)358 MRI_IMARR * mri_get_indexed_nbhd( MRI_IMAGE *inim , byte *mask ,
359 int xx, int yy, int zz, MCW_cluster *nbhd )
360 {
361 int kind , nx,ny,nz,nxy,nxyz , npt , nout , aa,bb,cc,kk,ii ;
362 MRI_IMAGE *im , *jm ; int *jar ;
363 MRI_IMARR *imar ;
364 void *brick ;
365
366 ENTRY("mri_get_indexed_nbhd") ;
367
368 if( inim == NULL || nbhd == NULL ) RETURN(NULL) ;
369
370 nx = inim->nx ;
371 ny = inim->ny ; nxy = nx*ny ;
372 nz = inim->nz ; nxyz = nxy*nz ; npt = nbhd->num_pt ; nout = 0 ;
373
374 kk = xx + yy*nx + zz*nxy ;
375 if (SearchAboutMaskedVoxel) {
376 if( npt == 0 || kk < 0 || kk >= nxyz ) RETURN(NULL);
377 } else {
378 if( npt == 0 || kk < 0 || kk >= nxyz || !INMASK(kk) ) RETURN(NULL);
379 }
380 kind = inim->kind ;
381 brick = mri_data_pointer(inim) ; if( brick == NULL ) RETURN(NULL);
382 im = mri_new( npt , 1 , kind ) ;
383 jm = mri_new( npt , 1 , MRI_int ) ; jar = MRI_INT_PTR(jm) ;
384
385 /*-- extract data, based on kind of data in sub-brick --*/
386
387 switch( kind ){
388
389 default: break ; /* bad bad bad user */
390
391 case MRI_short:{
392 short *rr = MRI_SHORT_PTR(im) , *vv = (short *)brick ;
393 for( ii=0 ; ii < npt ; ii++ ){
394 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
395 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
396 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
397 kk = aa + bb*nx + cc*nxy ;
398 if( INMASK(kk) ){ jar[nout] = kk; rr[nout++] = vv[kk]; }
399 }
400 }
401 break ;
402
403 case MRI_byte:{
404 byte *rr = MRI_BYTE_PTR(im) , *vv = (byte *)brick ;
405 for( ii=0 ; ii < npt ; ii++ ){
406 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
407 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
408 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
409 kk = aa + bb*nx + cc*nxy ;
410 if( INMASK(kk) ){ jar[nout] = kk; rr[nout++] = vv[kk]; }
411 }
412 }
413 break ;
414
415 case MRI_float:{
416 float *rr = MRI_FLOAT_PTR(im) , *vv = (float *)brick ;
417 for( ii=0 ; ii < npt ; ii++ ){
418 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
419 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
420 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
421 kk = aa + bb*nx + cc*nxy ;
422 if( INMASK(kk) ){ jar[nout] = kk; rr[nout++] = vv[kk]; }
423 }
424 thd_floatscan(nout,rr) ; /* 10 Jun 2021 */
425 }
426 break ;
427
428 case MRI_complex:{
429 complex *rr = MRI_COMPLEX_PTR(im) , *vv = (complex *)brick ;
430 for( ii=0 ; ii < npt ; ii++ ){
431 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
432 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
433 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
434 kk = aa + bb*nx + cc*nxy ;
435 if( INMASK(kk) ){ jar[nout] = kk; rr[nout++] = vv[kk]; }
436 }
437 thd_complexscan(nout,rr) ; /* 10 Jun 2021 */
438 }
439 break ;
440
441 case MRI_rgb:{
442 byte *rr = MRI_BYTE_PTR(im) , *vv = (byte *)brick ;
443 for( ii=0 ; ii < npt ; ii++ ){
444 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
445 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
446 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
447 kk = aa + bb*nx + cc*nxy ;
448 if( INMASK(kk) ){
449 jar[nout] = kk ;
450 rr[3*nout ] = vv[3*kk ] ;
451 rr[3*nout+1] = vv[3*kk+1] ;
452 rr[3*nout+2] = vv[3*kk+2] ; nout++ ;
453 }
454 }
455 }
456 break ;
457 }
458
459 if( nout == 0 ){ mri_free(im); mri_free(jm); RETURN(NULL); }
460
461 jm->nx = jm->nvox = im->nx = im->nvox = nout ;
462 INIT_IMARR(imar) ;
463 ADDTO_IMARR(imar,im) ; ADDTO_IMARR(imar,jm) ; RETURN(imar) ;
464 }
465