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