1 #include "mrilib.h"
2 #include "pbardefs.h"
3
4 /********************************************************************
5 ****** Functions to create and deal with SUMA_surface structs ******
6 ********************************************************************/
7
8 #define SUMA_EXTEND_NUM 64
9 #define SUMA_EXTEND_FAC 1.05
10
11 /*------------------------------------------------------------------*/
12 /*! Create an empty surface mask description.
13 --------------------------------------------------------------------*/
14
SUMA_create_empty_mask(void)15 SUMA_mask * SUMA_create_empty_mask(void)
16 {
17 SUMA_mask *msk ;
18
19 ENTRY("SUMA_create_empty_mask") ;
20
21 msk = (SUMA_mask *)calloc(1,sizeof(SUMA_mask)) ;
22 msk->type = SUMA_MASK_TYPE ;
23
24 msk->idcode[0] = '\0' ;
25 msk->num_surf = 0 ;
26 msk->idcode_surf = NULL ;
27 msk->surf = NULL ;
28 LOAD_FVEC3( msk->init_cen , 0,0,0 ) ;
29 LOAD_FVEC3( msk->show_cen , 0,0,0 ) ;
30
31 RETURN(msk) ;
32 }
33
34 /*------------------------------------------------------------------*/
35
SUMA_destroy_mask(SUMA_mask * msk,int kill_surfaces_too)36 void SUMA_destroy_mask( SUMA_mask *msk , int kill_surfaces_too )
37 {
38 int ss ;
39
40 ENTRY("SUMA_destroy_mask") ;
41
42 if( msk == NULL ) EXRETURN ;
43
44 if( msk->idcode_surf != NULL ){
45 for( ss=0 ; ss < msk->num_surf ; ss++ )
46 if( msk->idcode_surf[ss] != NULL ) free(msk->idcode_surf[ss]) ;
47 free(msk->idcode_surf) ;
48 }
49
50 if( msk->surf != NULL ){
51 if( kill_surfaces_too ){
52 for( ss=0 ; ss < msk->num_surf ; ss++ )
53 SUMA_destroy_surface(msk->surf[ss]) ;
54 }
55 free(msk->surf) ;
56 }
57
58 free(msk) ; EXRETURN ;
59 }
60
61 /*------------------------------------------------------------------*/
62 /*! Create an empty surface description.
63 --------------------------------------------------------------------*/
64
SUMA_create_empty_surface(void)65 SUMA_surface * SUMA_create_empty_surface(void)
66 {
67 SUMA_surface *ag ;
68
69 ENTRY("SUMA_create_empty_surface") ;
70
71 ag = (SUMA_surface *) calloc(1,sizeof(SUMA_surface)) ;
72 ag->type = SUMA_SURFACE_TYPE ;
73
74 ag->num_ixyz = ag->num_ijk = 0 ;
75 ag->nall_ixyz = ag->nall_ijk = 1 ;
76 ag->ixyz = (SUMA_ixyz *) malloc(sizeof(SUMA_ixyz)) ; /* space for */
77 ag->ijk = (SUMA_ijk *) malloc(sizeof(SUMA_ijk) ) ; /* 1 of each */
78 ag->norm = NULL ; ; /* none of this */
79
80 if( ag->ixyz == NULL || ag->ijk == NULL )
81 ERROR_exit("SUMA_create_empty_surface: can't malloc") ;
82
83 ag->idcode[0] = ag->idcode_dset[0] = ag->idcode_ldp[0] =
84 ag->label[0] = ag->label_ldp[0] = '\0' ;
85
86 ag->xbot = ag->ybot = ag->zbot = WAY_BIG ;
87 ag->xtop = ag->ytop = ag->ztop = -WAY_BIG ;
88 ag->xcen = ag->ycen = ag->zcen = 0.0 ;
89
90 ag->seq = ag->seqbase = ag->sorted = 0 ; /* not sequential; not sorted */
91
92 ag->vv = NULL ; /* 16 Jun 2003 */
93 ag->vn = NULL ; /* 22 Jan 2004 */
94
95 ag->mask_code = 0 ; /* mask stuff [04 Apr 2014] */
96 ag->mask_parent_idcode[0] = '\0' ;
97 ag->line_color[0] = '\0' ;
98 ag->box_color[0] = '\0' ;
99 ag->line_width = 0 ;
100 ag->mask = NULL ;
101 ag->parent = NULL ;
102
103 RETURN( ag ) ;
104 }
105
106 /*------------------------------------------------------------------*/
107 /*! Throw out some trash (i.e., free the contents of a surface).
108 --------------------------------------------------------------------*/
109
SUMA_destroy_surface(SUMA_surface * ag)110 void SUMA_destroy_surface( SUMA_surface *ag )
111 {
112 ENTRY("SUMA_destroy_surface") ;
113
114 if( ag == NULL ) EXRETURN ;
115 if( ag->ixyz != NULL ) free((void *)ag->ixyz) ;
116 if( ag->ijk != NULL ) free((void *)ag->ijk) ;
117 if( ag->norm != NULL ) free((void *)ag->norm) ;
118
119 if( ag->vv != NULL ) DESTROY_VVLIST(ag->vv) ;
120 if( ag->vn != NULL ) SUMA_destroy_vnlist(ag->vn) ;
121
122 free((void *)ag) ; EXRETURN ;
123 }
124
125 /*------------------------------------------------------------------*/
126 /* Add a bunch of nodes to a surface.
127 --------------------------------------------------------------------*/
128
SUMA_add_nodes_ixyz(SUMA_surface * ag,int nadd,int * iadd,float * xadd,float * yadd,float * zadd)129 void SUMA_add_nodes_ixyz( SUMA_surface *ag, int nadd,
130 int *iadd, float *xadd, float *yadd, float *zadd )
131 {
132 int ii , nup ;
133
134 ENTRY("SUMA_add_nodes_ixyz") ;
135
136 if( ag == NULL || nadd < 1 ) EXRETURN ;
137 if( xadd == NULL || yadd == NULL || zadd == NULL || iadd == NULL ) EXRETURN ;
138
139 nup = ag->num_ixyz + nadd ;
140
141 if( nup >= SUMA_MAX_NODES ){ /* 07 Sep 2001 */
142 fprintf(stderr,
143 "** SUMA surface can't have more than %d nodes!\n",
144 SUMA_MAX_NODES-1 ) ;
145 EXRETURN ;
146 }
147
148 if( nup > ag->nall_ixyz ){ /* extend length of array */
149 ag->nall_ixyz = nup = nup*SUMA_EXTEND_FAC + SUMA_EXTEND_NUM ;
150 ag->ixyz = (SUMA_ixyz *) realloc( (void *)ag->ixyz, sizeof(SUMA_ixyz)*nup );
151 if( ag->ixyz == NULL ){
152 fprintf(stderr,"SUMA_add_nodes_ixyz: can't malloc!\n"); EXIT(1);
153 }
154 }
155
156 nup = ag->num_ixyz ;
157
158 for( ii=0 ; ii < nadd ; ii++ ){
159 ag->ixyz[ii+nup].x = xadd[ii] ;
160 ag->ixyz[ii+nup].y = yadd[ii] ;
161 ag->ixyz[ii+nup].z = zadd[ii] ;
162 ag->ixyz[ii+nup].id = iadd[ii] ;
163 }
164
165 ag->num_ixyz += nadd ;
166
167 ag->seq = ag->sorted = 0 ; EXRETURN ;
168 }
169
170 /*--------------------------------------------------------------------
171 * Add/replace normals on the given surface. 05 Oct 2004 [rickr]
172 *
173 * This function requires one normal per node, and that the
174 * indices match.
175 *--------------------------------------------------------------------
176 */
SUMA_add_norms_xyz(SUMA_surface * ag,int nadd,float * xadd,float * yadd,float * zadd)177 int SUMA_add_norms_xyz( SUMA_surface *ag, int nadd,
178 float *xadd, float *yadd, float *zadd )
179 {
180 int ii ;
181
182 ENTRY("SUMA_add_norms_xyz") ;
183
184 if( ag == NULL || nadd < 1 ) RETURN(-1) ;
185 if( xadd == NULL || yadd == NULL || zadd == NULL ) RETURN(-1) ;
186
187 if( nadd != ag->num_ixyz ){
188 fprintf(stderr, "** SUMA surface has %d nodes but %d normals!\n",
189 ag->num_ixyz, nadd ) ;
190 RETURN(-1) ;
191 }
192
193 /* if norm is NULL, memory is needed */
194 if( ag->norm == NULL ){
195 ag->norm = (THD_fvec3 *)calloc(nadd, sizeof(THD_fvec3));
196 if( ag->norm == NULL ){
197 fprintf(stderr,"SUMA_add_norms_xyz: can't malloc!\n"); EXIT(1);
198 }
199 }
200
201 for( ii=0 ; ii < nadd ; ii++ ){
202 ag->norm[ii].xyz[0] = xadd[ii] ;
203 ag->norm[ii].xyz[1] = yadd[ii] ;
204 ag->norm[ii].xyz[2] = zadd[ii] ;
205 }
206
207 RETURN(0) ;
208 }
209
210 /*------------------------------------------------------------------*/
211 /*! Add 1 pitiful node to a surface.
212 --------------------------------------------------------------------*/
213
SUMA_add_node_ixyz(SUMA_surface * ag,int i,float x,float y,float z)214 void SUMA_add_node_ixyz( SUMA_surface *ag , int i,float x,float y,float z )
215 {
216 SUMA_add_nodes_ixyz( ag , 1 , &i,&x,&y,&z ) ;
217 }
218
219 /*------------------------------------------------------------------*/
220 /*! Erase all triangles from a surface.
221 --------------------------------------------------------------------*/
222
SUMA_clear_triangles(SUMA_surface * ag)223 void SUMA_clear_triangles( SUMA_surface *ag )
224 {
225
226 ENTRY("SUMA_clear_triangles") ;
227
228 if( ag == NULL || ag->num_ijk <= 0 ) EXRETURN ;
229
230 free(ag->ijk) ; ag->ijk = NULL ; ag->num_ijk = 0 ;
231 EXRETURN ;
232 }
233
234 /*------------------------------------------------------------------*/
235
SUMA_clear_normals(SUMA_surface * ag)236 void SUMA_clear_normals( SUMA_surface *ag )
237 {
238
239 ENTRY("SUMA_clear_normals") ;
240
241 if( ag == NULL || ag->norm == NULL ) EXRETURN ;
242
243 free(ag->norm) ; ag->norm = NULL ;
244 EXRETURN ;
245 }
246
247 /*------------------------------------------------------------------*/
248 /*! Add a bunch of triangles (node id triples) to a surface.
249 --------------------------------------------------------------------*/
250
SUMA_add_triangles(SUMA_surface * ag,int nadd,int * it,int * jt,int * kt)251 void SUMA_add_triangles( SUMA_surface *ag, int nadd, int *it, int *jt, int *kt )
252 {
253 int ii , nup ;
254
255 ENTRY("SUMA_add_triangles") ;
256
257 if( ag == NULL || nadd < 1 ) EXRETURN ;
258 if( it == NULL || jt == NULL || kt == NULL ) EXRETURN ;
259
260 nup = ag->num_ijk + nadd ;
261 if( nup > ag->nall_ijk ){ /* extend length of array */
262 ag->nall_ijk = nup = nup*SUMA_EXTEND_FAC + SUMA_EXTEND_NUM ;
263 ag->ijk = (SUMA_ijk *) realloc( (void *)ag->ijk , sizeof(SUMA_ijk)*nup ) ;
264 if( ag->ijk == NULL ){
265 fprintf(stderr,"SUMA_add_triangles: can't malloc!\n"); EXIT(1);
266 }
267 }
268
269 nup = ag->num_ijk ;
270 for( ii=0 ; ii < nadd ; ii++ ){
271 ag->ijk[ii+nup].id = it[ii] ;
272 ag->ijk[ii+nup].jd = jt[ii] ;
273 ag->ijk[ii+nup].kd = kt[ii] ;
274 }
275
276 ag->num_ijk += nadd ; EXRETURN ;
277 }
278
279 /*------------------------------------------------------------------*/
280 /*! Add 1 pitiful triangle to a surface.
281 --------------------------------------------------------------------*/
282
SUMA_add_triangle(SUMA_surface * ag,int it,int jt,int kt)283 void SUMA_add_triangle( SUMA_surface *ag, int it, int jt, int kt )
284 {
285 SUMA_add_triangles( ag , 1 , &it,&jt,&kt ) ;
286 }
287
288 /*------------------------------------------------------------------*/
289 /*! Truncate the memory used by the node and triangle arrays back
290 to the minimum they need.
291 --------------------------------------------------------------------*/
292
SUMA_truncate_memory(SUMA_surface * ag)293 void SUMA_truncate_memory( SUMA_surface *ag )
294 {
295 int nn ;
296
297 ENTRY("SUMA_truncate_memory") ;
298
299 if( ag == NULL ) EXRETURN ;
300
301 if( ag->num_ixyz < ag->nall_ixyz && ag->num_ixyz > 0 ){
302 ag->nall_ixyz = nn = ag->num_ixyz ;
303 ag->ixyz = (SUMA_ixyz *) realloc( (void *)ag->ixyz, sizeof(SUMA_ixyz)*nn );
304 }
305
306 if( ag->num_ijk < ag->nall_ijk && ag->num_ijk > 0 ){
307 ag->nall_ijk = nn = ag->num_ijk ;
308 ag->ijk = (SUMA_ijk *) realloc( (void *)ag->ijk , sizeof(SUMA_ijk)*nn ) ;
309 }
310
311 EXRETURN ;
312 }
313
314 /*------------------------------------------------------------------
315 Generate a function to sort array of SUMA_ixyz-s by their id-s
316 --------------------------------------------------------------------*/
317
318 #undef STYPE
319 #define STYPE SUMA_ixyz /* name of type to sort */
320 #define SLT(a,b) ((a).id < (b).id) /* macro to decide on order */
321 #define SNAME SUMA_ixyz /* function qsort_SUMA_ixyz */
322 #include "cs_sort_template.h" /* generate the function */
323 #undef STYPE
324
325 /*** Can now use qsort_SUMA_ixyz(int,SUMA_ixyz *) ***/
326
327 /*------------------------------------------------------------------*/
328 /*! Sort the nodes by id-s, and mark if the id-s are sequential.
329 --------------------------------------------------------------------*/
330
SUMA_ixyzsort_surface(SUMA_surface * ag)331 void SUMA_ixyzsort_surface( SUMA_surface *ag )
332 {
333 int nn , ii , ndup ;
334 float xb,yb,zb , xt,yt,zt , xc,yc,zc ;
335
336 ENTRY("SUMA_ixyzsort_surface") ;
337
338 if( ag == NULL || ag->num_ixyz < 1 ) EXRETURN ;
339
340 SUMA_truncate_memory( ag ) ;
341
342 nn = ag->num_ixyz ;
343
344 /* check if nodes are already sorted [26 Oct 2001] */
345
346 for( ii=1 ; ii < nn ; ii++ )
347 if( ag->ixyz[ii].id <= ag->ixyz[ii-1].id ) break ;
348
349 /* if not in increasing order,
350 sort them using the function generated above */
351
352 if( ii < nn ){
353 qsort_SUMA_ixyz( nn , ag->ixyz ) ;
354 }
355
356 ag->sorted = 1 ; /* mark as sorted */
357
358 /* check if node id-s are sequential */
359
360 for( ii=1 ; ii < nn ; ii++ )
361 if( ag->ixyz[ii].id != ag->ixyz[ii-1].id+1 ) break ;
362
363 /* if we finished that loop all the way,
364 mark the nodes as being sequential, and
365 store the base of the sequence (id of node #0) */
366
367 if( ii == nn ){
368 ag->seq = 1 ; ag->seqbase = ag->ixyz[0].id ;
369 }
370
371 /* 07 Sep 2001: check for duplicate node id-s */
372
373 for( ndup=0,ii=1 ; ii < nn ; ii++ )
374 if( ag->ixyz[ii].id == ag->ixyz[ii-1].id ) ndup++ ;
375
376 if( ndup > 0 )
377 fprintf(stderr,"** SUMA WARNING: %d duplicate surface node id's found!\n",ndup);
378
379 /* find bounding box of all nodes (it's useful on occasion) */
380
381 xb = xt = ag->ixyz[0].x ;
382 yb = yt = ag->ixyz[0].y ;
383 zb = zt = ag->ixyz[0].z ;
384 xc = yc = zc = 0.0 ;
385 for( ii=1 ; ii < nn ; ii++ ){
386 xc += ag->ixyz[ii].x ;
387 yc += ag->ixyz[ii].y ;
388 zc += ag->ixyz[ii].z ;
389
390 if( ag->ixyz[ii].x < xb ) xb = ag->ixyz[ii].x ;
391 else if( ag->ixyz[ii].x > xt ) xt = ag->ixyz[ii].x ;
392
393 if( ag->ixyz[ii].y < yb ) yb = ag->ixyz[ii].y ;
394 else if( ag->ixyz[ii].y > yt ) yt = ag->ixyz[ii].y ;
395
396 if( ag->ixyz[ii].z < zb ) zb = ag->ixyz[ii].z ;
397 else if( ag->ixyz[ii].z > zt ) zt = ag->ixyz[ii].z ;
398 }
399
400 ag->xbot = xb ; ag->xtop = xt ;
401 ag->ybot = yb ; ag->ytop = yt ;
402 ag->zbot = zb ; ag->ztop = zt ;
403
404 ag->xcen = xc/nn ; ag->ycen = yc/nn ; ag->zcen = zc/nn ;
405
406 EXRETURN ;
407 }
408
409 /*--------------------------------------------------------------------*/
410 /*! Find a node id in a surface, and return its index into the node
411 array; return -1 if not found.
412 ----------------------------------------------------------------------*/
413
SUMA_find_node_id(SUMA_surface * ag,int target)414 int SUMA_find_node_id( SUMA_surface *ag , int target )
415 {
416 int nn , ii,jj,kk ;
417
418 if( ag == NULL || ag->num_ixyz < 1 || target < 0 ) return( -1 );
419
420 if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
421
422 if( ag->seq ){ /* node id-s are sequential (the easy case) */
423 kk = target - ag->seqbase ;
424 if( kk >= 0 && kk < ag->num_ixyz ) return( kk );
425 return( -1 );
426 }
427
428 /* node id-s are in increasing order, but not sequential;
429 so, use binary search to find the node id (if present) */
430
431 ii = 0 ; jj = ag->num_ixyz - 1 ; /* search bounds */
432
433 if( target < ag->ixyz[0].id ) return( -1 ); /* not present */
434 else if( target == ag->ixyz[0].id ) return( ii ); /* at start! */
435
436 if( target > ag->ixyz[jj].id ) return( -1 ); /* not present */
437 else if( target == ag->ixyz[jj].id ) return( jj ); /* at end! */
438
439 while( jj - ii > 1 ){ /* while search bounds not too close */
440
441 kk = (ii+jj) / 2 ; /* midway between search bounds */
442
443 nn = ag->ixyz[kk].id - target ;
444 if( nn == 0 ) return( kk ); /* AHA! */
445
446 if( nn < 0 ) ii = kk ; /* kk before target => bottom = kk */
447 else jj = kk ; /* kk after target => top = kk */
448 }
449
450 return( -1 );
451 }
452
453 /*-------------------------------------------------------------------------*/
454 /*! Create the voxel-to-node list for this surface/dataset combo.
455 ---------------------------------------------------------------------------*/
456
SUMA_make_vnlist(SUMA_surface * ag,THD_3dim_dataset * dset)457 SUMA_vnlist * SUMA_make_vnlist( SUMA_surface *ag , THD_3dim_dataset *dset )
458 {
459 int ii,jj,kk , nx,ny,nz , nxy,nxyz , nnode , pp,qq,nn,nvox ;
460 THD_fvec3 fv ;
461 THD_ivec3 iv ;
462 int *vlist , *nlist , wodsave ;
463 SUMA_vnlist *vnlist ;
464 float xbot,xtop , ybot,ytop , zbot,ztop ;
465
466 ENTRY("SUMA_make_vnlist") ;
467
468 if( ag == NULL || ag->num_ixyz < 1 || !ISVALID_DSET(dset) ) RETURN(NULL) ;
469
470 if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
471
472 /* setup: create arrays for voxel list and node list */
473
474 nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
475 nxy = nx*ny ; nxyz = nxy*nz ; nnode = ag->num_ixyz ;
476 vlist = (int *) malloc(sizeof(int)*nnode) ;
477 nlist = (int *) malloc(sizeof(int)*nnode) ;
478 if( vlist == NULL || nlist == NULL ){
479 fprintf(stderr,"SUMA_make_vnlist: can't malloc!\n"); EXIT(1);
480 }
481
482 /* for each node, find which voxel it is in */
483
484 wodsave = dset->wod_flag ; dset->wod_flag = 0 ;
485
486 xbot = DSET_XXMIN(dset) ; xtop = DSET_XXMAX(dset) ;
487 ybot = DSET_YYMIN(dset) ; ytop = DSET_YYMAX(dset) ;
488 zbot = DSET_ZZMIN(dset) ; ztop = DSET_ZZMAX(dset) ;
489
490 for( nn=pp=0 ; pp < nnode ; pp++ ){
491 LOAD_FVEC3( fv , ag->ixyz[pp].x, ag->ixyz[pp].y, ag->ixyz[pp].z ) ;
492 fv = THD_dicomm_to_3dmm( dset , fv ) ; /* convert Dicom coords */
493
494 if( fv.xyz[0] < xbot || fv.xyz[0] > xtop ) continue ;
495 if( fv.xyz[1] < ybot || fv.xyz[1] > ytop ) continue ;
496 if( fv.xyz[2] < zbot || fv.xyz[2] > ztop ) continue ;
497
498 iv = THD_3dmm_to_3dind( dset , fv ) ; /* in surface to */
499 UNLOAD_IVEC3( iv , ii,jj,kk ) ; /* dataset indexes */
500
501 nlist[nn] = pp ; /* list of nodes */
502 vlist[nn] = ii + jj*nx + kk*nxy ; /* list of voxels */
503 nn++ ;
504 }
505
506 nnode = nn ; /* number of nodes inside dataset volume */
507 if( nnode == 0 ){ free(nlist); free(vlist); RETURN(NULL); }
508
509 dset->wod_flag = wodsave ;
510
511 /* now sort the 2 lists so that vlist is increasing
512 (and each nlist still corresponds to its original vlist) */
513
514 qsort_intint( nnode , vlist , nlist ) ;
515
516 /* count how many distinct voxels we found */
517
518 nvox = 1 ; ii = vlist[0] ;
519 for( pp=1 ; pp < nnode ; pp++ ){
520 if( vlist[pp] != ii ){ nvox++; ii = vlist[pp]; }
521 }
522
523 /* now create the output vnlist */
524
525 /* from malloc (no effect now) 12 Feb 2009 [lesstif patrol] */
526 vnlist = (SUMA_vnlist *) calloc( 1, sizeof(SUMA_vnlist) ) ;
527 vnlist->nvox = nvox ;
528 vnlist->voxijk = (int *) malloc(sizeof(int) *nvox) ;
529 vnlist->numnod = (int *) calloc(sizeof(int) ,nvox) ;
530 vnlist->nlist = (int **)malloc(sizeof(int*)*nvox);
531 vnlist->dset = dset ;
532
533 if( vnlist->voxijk==NULL || vnlist->numnod==NULL || vnlist->nlist==NULL ){
534 fprintf(stderr,"SUMA_make_vnlist: can't malloc!\n"); EXIT(1);
535 }
536
537 /* now count how many nodes are at each voxel in the list */
538
539 ii = vlist[0] ; qq = nn = 0 ;
540 for( pp=1 ; pp < nnode ; pp++ ){
541 if( vlist[pp] != ii ){ /* qq..pp-1 are the same */
542 vnlist->voxijk[nn] = ii ;
543 vnlist->numnod[nn] = jj = pp-qq ;
544 vnlist->nlist[nn] = (int *) malloc(sizeof(int)*jj) ;
545 memcpy( vnlist->nlist[nn] , nlist+qq , sizeof(int)*jj ) ;
546 ii = vlist[pp] ; nn++ ; qq = pp ;
547 }
548 }
549 vnlist->voxijk[nn] = ii ;
550 vnlist->numnod[nn] = jj = pp-qq ;
551 vnlist->nlist[nn] = (int *) malloc(sizeof(int)*jj) ;
552 memcpy( vnlist->nlist[nn] , nlist+qq , sizeof(int)*jj ) ;
553
554 /* and we're done! */
555
556 free(nlist) ; free(vlist) ; RETURN( vnlist ) ;
557 }
558
559 /*-------------------------------------------------------------------------*/
560 /*! Destroy a SUMA_vnlist struct.
561 ---------------------------------------------------------------------------*/
562
SUMA_destroy_vnlist(SUMA_vnlist * vnlist)563 void SUMA_destroy_vnlist( SUMA_vnlist *vnlist )
564 {
565 int ii ;
566 if( vnlist == NULL ) return ;
567 if( vnlist->voxijk != NULL ) free( vnlist->voxijk ) ;
568 if( vnlist->numnod != NULL ) free( vnlist->numnod ) ;
569 if( vnlist->nlist != NULL ){
570 for( ii=0 ; ii < vnlist->nvox ; ii++ )
571 if( vnlist->nlist[ii] != NULL ) free( vnlist->nlist[ii] ) ;
572 free( vnlist->nlist ) ;
573 }
574 free( vnlist ) ;
575 }
576
577 /*--------------------------------------------------------------------------
578 The following routines are used to convert DICOM order coordinates
579 (used in AFNI) to SureFit order coordinates -- 25 Oct 2001 - RWCox
580 ----------------------------------------------------------------------------*/
581
THD_dicomm_to_surefit(THD_3dim_dataset * dset,THD_fvec3 fv)582 THD_fvec3 THD_dicomm_to_surefit( THD_3dim_dataset *dset , THD_fvec3 fv )
583 {
584 float xx,yy,zz , xbase,ybase,zbase ;
585 THD_fvec3 vout ;
586
587 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ; /* xyz now LPI */
588
589 if( dset != NULL ){
590 THD_fvec3 v1 , v2 ;
591 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
592 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
593 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
594 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
595 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
596 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
597 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ; /* Left-most */
598 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ; /* Posterior */
599 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ; /* Inferior */
600 } else {
601 xbase = ybase = zbase = 0.0 ;
602 }
603
604 vout.xyz[0] = xx - xbase ;
605 vout.xyz[1] = yy - ybase ;
606 vout.xyz[2] = zz - zbase ; return vout ;
607 }
608
609 /* --------------------------------------------------------------------------*/
610
THD_surefit_to_dicomm(THD_3dim_dataset * dset,THD_fvec3 fv)611 THD_fvec3 THD_surefit_to_dicomm( THD_3dim_dataset *dset , THD_fvec3 fv )
612 {
613 float xx,yy,zz , xbase,ybase,zbase ;
614 THD_fvec3 vout ;
615
616 xx = -fv.xyz[0] ; yy = -fv.xyz[1] ; zz = fv.xyz[2] ; /* xyz now RAI */
617
618 if( dset != NULL ){
619 THD_fvec3 v1 , v2 ;
620 LOAD_FVEC3(v1, DSET_XORG(dset),DSET_YORG(dset),DSET_ZORG(dset)) ;
621 v1 = THD_3dmm_to_dicomm( dset , v1 ) ;
622 LOAD_FVEC3(v2, DSET_XORG(dset)+(DSET_NX(dset)-1)*DSET_DX(dset) ,
623 DSET_YORG(dset)+(DSET_NY(dset)-1)*DSET_DY(dset) ,
624 DSET_ZORG(dset)+(DSET_NZ(dset)-1)*DSET_DZ(dset) ) ;
625 v2 = THD_3dmm_to_dicomm( dset , v2 ) ;
626 xbase = MAX( v1.xyz[0] , v2.xyz[0] ) ; xbase = -xbase ;
627 ybase = MAX( v1.xyz[1] , v2.xyz[1] ) ; ybase = -ybase ;
628 zbase = MIN( v1.xyz[2] , v2.xyz[2] ) ;
629 } else {
630 xbase = ybase = zbase = 0.0 ;
631 }
632
633 vout.xyz[0] = xx - xbase ;
634 vout.xyz[1] = yy - ybase ;
635 vout.xyz[2] = zz + zbase ; return vout ;
636 }
637
638 /****************************************************************************
639 ********** AFNI no longer reads surface from files [22 Jan 2004] ***********
640 ****************************************************************************/
641
642 #undef ALLOW_SURFACE_FILES
643 #ifdef ALLOW_SURFACE_FILES
644 /*----------------------------------------------------------*/
645 /*! Will load this many items (nodes, triangles) at a time. */
646
647 #undef NBUF
648 #define NBUF 64
649
650 /*------------------------------------------------------------------------*/
651 /*! Read a surface description file that is associated with an AFNI
652 dataset. Return a surface ready to rock-n-roll.
653 NULL is returned if something bad happens.
654 --------------------------------------------------------------------------*/
655
SUMA_read_surface(char * fname,THD_3dim_dataset * dset)656 SUMA_surface * SUMA_read_surface( char *fname , THD_3dim_dataset *dset )
657 {
658 SUMA_surface *ag ;
659 FILE *fp ;
660 char lbuf[1024] , *cpt ;
661 int do_nod=1 , ii ;
662 float xx[NBUF],yy[NBUF],zz[NBUF] ;
663 int pp[NBUF],qq[NBUF],rr[NBUF] , nn ;
664
665 THD_vecmat mv ;
666 int have_mv=0 ;
667
668 ENTRY("SUMA_read_surface") ;
669
670 if( fname == NULL || fname[0] == '\0' ) RETURN( NULL );
671
672 /*-- open input --*/
673
674 if( strcmp(fname,"-") == 0 ){ /* special case */
675 fp = stdin ;
676 } else {
677 fp = fopen( fname , "r" ) ;
678 if( fp == NULL ) RETURN( NULL );
679 }
680
681 /*-- initialize surface that we will fill up here */
682
683 ag = SUMA_create_empty_surface() ;
684
685 /*-- set IDCODEs of surface (from filename) and of its dataset --*/
686
687 cpt = UNIQ_hashcode(fname); strcpy(ag->idcode,cpt); free(cpt);
688
689 strcpy( ag->idcode_dset , dset->idcode.str ) ;
690
691 MCW_strncpy( ag->label, THD_trailname(fname,0), 32 ) ; /* 19 Aug 2002 */
692
693 /*-- read data --*/
694
695 nn = 0 ;
696 while(1){
697 cpt = afni_fgets(lbuf,1024,fp) ; /* read a line */
698 if( cpt == NULL ) break ; /* end of file */
699
700 /*-- read a transformation matrix-vector? --*/
701
702 if( strncmp(lbuf,"<MATVEC>",8) == 0 ){ /* 07 Sep 2001 */
703 float a11,a12,a13 , v1 ,
704 a21,a22,a23 , v2 ,
705 a31,a32,a33 , v3 ;
706
707 ii = sscanf( lbuf+8 , "%f%f%f%f%f%f%f%f%f%f%f%f" ,
708 &a11,&a12,&a13 , &v1 ,
709 &a21,&a22,&a23 , &v2 ,
710 &a31,&a32,&a33 , &v3 ) ;
711
712 if( ii < 12 ){
713 fprintf(stderr,"** SUMA: Illegal MATVEC in %s\n",fname) ;
714 have_mv = 0 ;
715 } else {
716 LOAD_FVEC3(mv.vv , v1,v2,v3 ) ;
717 LOAD_MAT (mv.mm , a11,a12,a13,a21,a22,a23,a31,a32,a33) ;
718 have_mv = 1 ;
719 }
720 continue ; /* skip to next line */
721 }
722
723 /*-- read data from SureFit? --*/
724
725 if( strncmp(lbuf,"<SureFit",8) == 0 ){ /* 26 Oct 2001 */
726
727 if( nn > 0 ){ /* process existing inputs, if any */
728 if( do_nod )
729 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
730 else
731 SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
732 nn = 0 ;
733 }
734
735 SUMA_import_surefit( ag , lbuf , dset ) ;
736 continue ; /* skip to next input line */
737
738 } /* end of SureFit input */
739
740 /*-- end of node input? --*/
741
742 if( strstr(lbuf,"</NODES>") != NULL ){
743 if( do_nod && nn > 0 ){
744 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
745 nn = 0 ;
746 }
747 #if 1
748 do_nod = 0 ; /* from now on, triangles */
749 continue ; /* skip to next line */
750 #else
751 break ; /* don't create triangles */
752 #endif
753 }
754
755 /*-- process line as a node? --*/
756
757 if( do_nod ){ /* nodes */
758
759 ii = sscanf(lbuf,"%d%f%f%f",pp+nn,xx+nn,yy+nn,zz+nn) ;
760 if( ii < 4 ) continue ;
761 if( have_mv ){ /* transform coords */
762 THD_fvec3 fv , qv ;
763 LOAD_FVEC3(fv , xx[nn],yy[nn],zz[nn] ) ;
764 qv = VECSUB_MAT( mv.mm , fv , mv.vv ) ;
765 UNLOAD_FVEC3( qv , xx[nn],yy[nn],zz[nn] ) ;
766 }
767 nn++ ;
768 if( nn == NBUF ){ /* add nodes in chunks */
769 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
770 nn = 0 ;
771 }
772
773 /*-- process line as a triangle --*/
774
775 } else { /* triangles */
776
777 ii = sscanf(lbuf,"%d%d%d",pp+nn,qq+nn,rr+nn) ;
778 if( ii < 3 ) continue ;
779 nn++ ;
780 if( nn == NBUF ){ /* add triangles in chunks */
781 SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
782 nn = 0 ;
783 }
784 }
785 } /* end of loop over input lines */
786
787 /*-- finish up, eh? --*/
788
789 if( fp != stdin ) fclose(fp) ;
790 if( nn > 0 ){ /* process leftover data lines */
791 if( do_nod )
792 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
793 else
794 SUMA_add_triangles( ag,nn , pp,qq,rr ) ;
795 }
796
797 /*-- hope we got something --*/
798
799 if( ag->num_ixyz < 1 ){
800 SUMA_destroy_surface(ag) ; RETURN(NULL) ;
801 }
802
803 /*-- sort the nodes (if needed), et cetera --*/
804
805 SUMA_ixyzsort_surface(ag) ;
806
807 /*-- done --*/
808
809 RETURN( ag );
810 }
811
812 /*-----------------------------------------------------------------------------*/
813
814 /*! 26 point 3x3x3 nbhd of {0,0,0}
815 (not including the central point itself) */
816
817 static int ip[26][3] = { {-1,-1,-1},{-1,-1, 0},{-1,-1, 1},
818 {-1, 0,-1},{-1, 0, 0},{-1, 0, 1},
819 {-1, 1,-1},{-1, 1, 0},{-1, 1, 1},
820 { 0,-1,-1},{ 0,-1, 0},{ 0,-1, 1},
821 { 0, 0,-1}, { 0, 0, 1},
822 { 0, 1,-1},{ 0, 1, 0},{ 0, 1, 1},
823 { 1,-1,-1},{ 1,-1, 0},{ 1,-1, 1},
824 { 1, 0,-1},{ 1, 0, 0},{ 1, 0, 1},
825 { 1, 1,-1},{ 1, 1, 0},{ 1, 1, 1} } ;
826
827 /*------------------------------------------------------------------------*/
828 /*! Given a surface and a dataset, compute a voxel-to-node map.
829 This is an int array (the return value) the size of a dataset
830 brick - one value per voxel. Each entry "v" is
831 - v < 0 means that voxel is not close to a surface node
832 - otherwise, the closest node (index in ag->ixyz, NOT id) to the
833 voxel is SUMA_VMAP_UNMASK(v)
834 - the "level" of that node is SUMA_VMAP_LEVEL(v),
835 where level is the number of nbhd expansions outward from the
836 voxel that were needed to hit a node (0<=level<=7).
837 - The node index (26 bit integer) is stored in bits 0..25 of v
838 - The level (3 bit integer) is stored in bits 26..28 of v
839 - Bits 29..30 are reserved for future purposes
840 - Bit 31 is the sign bit, and so signifies "no node nearby"
841
842 Other notes:
843 - The input surface should have been sorted by SUMA_ixyzsort_surface().
844 If it wasn't, it will be now.
845 - Although this function was crafted to be efficient, it still takes
846 a few seconds to run.
847 - Don't add nodes to the surface after calling this, since they
848 won't be indexed here properly. Or you could free() the output
849 map and re-invoke this function.
850 ---------------------------------------------------------------------------*/
851
SUMA_map_dset_to_surf(SUMA_surface * ag,THD_3dim_dataset * dset)852 int * SUMA_map_dset_to_surf( SUMA_surface *ag , THD_3dim_dataset *dset )
853 {
854 int *vmap , ii,jj,kk , nx,ny,nz , nxy,nxyz , pp,qq,pbest ;
855 THD_fvec3 fv ;
856 THD_ivec3 iv ;
857 int ibot,jbot,kbot , itop,jtop,ktop , lev , ijk ;
858 float xv,yv,zv , dd,dbest=0 , xp,yp,zp ;
859 char *elev ; int ltop , ntop , lmask ;
860
861 ENTRY("SUMA_map_dset_to_surf") ;
862
863 if( ag == NULL || ag->num_ixyz < 1 || !ISVALID_DSET(dset) ) RETURN( NULL );
864
865 if( !ag->sorted ) SUMA_ixyzsort_surface( ag ) ;
866
867 /* setup: create the output vmap and fill it with -1 */
868
869 nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
870 nxy = nx*ny ; nxyz = nxy*nz ;
871 vmap = (int *) malloc(sizeof(int)*nxyz) ;
872 if( vmap == NULL ){
873 fprintf(stderr,"SUMA_map_dset_to_surf: can't malloc!\n"); EXIT(1);
874 }
875 for( ii=0 ; ii < nxyz ; ii++ ) vmap[ii] = -1 ; /* not mapped yet */
876
877 /* put nodes directly into voxels (level 0) */
878
879 STATUS("putting nodes into voxels") ;
880
881 for( pp=0 ; pp < ag->num_ixyz ; pp++ ){
882 LOAD_FVEC3( fv , ag->ixyz[pp].x, ag->ixyz[pp].y, ag->ixyz[pp].z ) ;
883 fv = THD_dicomm_to_3dmm( dset , fv ) ; /* convert Dicom coords */
884 iv = THD_3dmm_to_3dind( dset , fv ) ; /* in surface to */
885 UNLOAD_IVEC3( iv , ii,jj,kk ) ; /* dataset indexes */
886 qq = vmap[ii+jj*nx+kk*nxy] ; /* previously mapped index? */
887 if( qq < 0 ){ /* not mapped before */
888 vmap[ii+jj*nx+kk*nxy] = pp ; /* index, not id */
889 } else {
890 LOAD_IVEC3(iv,ii,jj,kk) ; /* get Dicom coords of voxel */
891 fv = THD_3dind_to_3dmm( dset , iv ) ;
892 fv = THD_3dmm_to_dicomm( dset , fv ) ;
893 UNLOAD_FVEC3(fv,xv,yv,zv) ; /* voxel is at (xv,yv,zv) */
894
895 /* find if old node in this voxel (#qq)
896 is closer to voxel center than current node (#pp) */
897
898 xp=xv-ag->ixyz[qq].x;
899 yp=yv-ag->ixyz[qq].y;
900 zp=zv-ag->ixyz[qq].z; dd=xp*xp+yp*yp+zp*zp; /* dist^2 to old node */
901
902 xp=xv-ag->ixyz[pp].x;
903 yp=yv-ag->ixyz[pp].y;
904 zp=zv-ag->ixyz[pp].z; dbest=xp*xp+yp*yp+zp*zp; /* dist^2 to new node */
905
906 if( dbest < dd ) vmap[ii+jj*nx+kk*nxy] = pp ; /* new is better */
907 }
908 }
909
910 LOAD_FVEC3( fv , ag->xbot,ag->ybot,ag->zbot ) ; /* find dataset */
911 fv = THD_dicomm_to_3dmm( dset , fv ) ; /* indexes for */
912 iv = THD_3dmm_to_3dind( dset , fv ) ; /* bot point */
913 UNLOAD_IVEC3( iv , ibot,jbot,kbot ) ; /* of surface */
914
915 LOAD_FVEC3( fv , ag->xtop,ag->ytop,ag->ztop ) ; /* and for top */
916 fv = THD_dicomm_to_3dmm( dset , fv ) ;
917 iv = THD_3dmm_to_3dind( dset , fv ) ;
918 UNLOAD_IVEC3( iv , itop,jtop,ktop ) ;
919
920 /* we will only try to fill voxels inside
921 the rectangular subvolume (ibot..itop,jbot..jtop,kbot..ktop) */
922
923 if( ibot > itop ){ ii=ibot ; ibot=itop; itop=ii; }
924 if( jbot > jtop ){ jj=jbot ; jbot=jtop; jtop=jj; }
925 if( kbot > ktop ){ kk=kbot ; kbot=ktop; ktop=kk; }
926 if( ibot < 1 ) ibot = 1 ; if( itop >= nx ) itop = nx-1 ;
927 if( jbot < 1 ) jbot = 1 ; if( jtop >= ny ) jtop = ny-1 ;
928 if( kbot < 1 ) kbot = 1 ; if( ktop >= nz ) ktop = nz-1 ;
929
930 /* larger than the largest node id */
931
932 ntop = ag->ixyz[ag->num_ixyz-1].id + 100 ;
933
934 /* scan for voxels that are next to those already mapped,
935 out to a given level (default max level is 4) */
936
937 elev = getenv("SUMA_NBHD_LEVEL") ; /* find level for expanding out */
938 if( elev != NULL ){
939 char *cpt ;
940 ltop = strtol( elev , &cpt , 10 ) ;
941 if( ltop < 0 || ltop > 7 || (ltop == 0 && *cpt != '\0') ) ltop = 4 ;
942 } else {
943 ltop = 4 ;
944 }
945
946 /* loop over levels */
947
948 for( lev=1 ; lev <= ltop ; lev++ ){ /* if ltop = 0, won't be executed */
949
950 if(PRINT_TRACING){
951 char str[256]; sprintf(str,"expansion level %d",lev); STATUS(str);
952 }
953
954 /* loop over voxel 3D indexes */
955
956 for( kk=kbot ; kk < ktop ; kk++ ){
957 for( jj=jbot ; jj < jtop ; jj++ ){
958 for( ii=ibot ; ii < itop ; ii++ ){
959
960 ijk = ii+jj*nx+kk*nxy ; /* index into brick array */
961 if( vmap[ijk] >= 0 ) continue ; /* already mapped */
962
963 LOAD_IVEC3(iv,ii,jj,kk) ; /* get Dicom coords */
964 fv = THD_3dind_to_3dmm( dset , iv ) ;
965 fv = THD_3dmm_to_dicomm( dset , fv ) ;
966 UNLOAD_FVEC3(fv,xv,yv,zv) ; /* coords = (xv,yv,zv) */
967
968 for( pbest=-1,qq=0 ; qq < 26 ; qq++ ){ /* loop over neighbors and */
969 /* find closest mapped pt */
970 /* (at a lower level) */
971
972 /* pp is the volume map at the qq'th neighboring voxel */
973
974 pp = vmap[(ii+ip[qq][0]) + (jj+ip[qq][1])*nx + (kk+ip[qq][2])*nxy];
975
976 if( pp >= 0 ){ /* if qq'th nbhr is mapped */
977 pp = SUMA_VMAP_UNMASK(pp) ; /* index of mapped node in qq */
978 xp=xv-ag->ixyz[pp].x ; /* coords of mapped node */
979 yp=yv-ag->ixyz[pp].y ;
980 zp=zv-ag->ixyz[pp].z ;
981 dd=xp*xp+yp*yp+zp*zp ; /* dist^2 to mapped pt */
982
983 /* save pbest = index of mapped node closest to (ii,jj,kk)
984 dbest = dist^2 of pbest to (ii,jj,kk) voxel center */
985
986 if( pbest >= 0 ){
987 if( dd < dbest ){ pbest = pp ; dbest = dd ; }
988 } else {
989 pbest = pp ; dbest = dd ;
990 }
991 }
992 } /* end of loop over 26 neighbors of (ii,jj,kk) voxel */
993
994 /* save closest of the neighbors;
995 temporarily as a large negative number,
996 so we won't hit it again in this level of expansion */
997
998 if( pbest >= 0 ) vmap[ijk] = -(pbest+ntop) ; /* closest of the neighbors */
999
1000 }}} /* end of loop over voxels (ii,jj,kk) */
1001
1002 STATUS(".. masking") ;
1003
1004 /* lmask = 3 bit int for level, shifted to bits 26..28 */
1005
1006 lmask = SUMA_VMAP_LEVMASK(lev) ; /* 07 Sep 2001: put on a mask */
1007 /* to indicate which level of */
1008 /* indirection this voxel was */
1009
1010 for( kk=kbot ; kk < ktop ; kk++ ){ /* change all the nodes we found */
1011 for( jj=jbot ; jj < jtop ; jj++ ){ /* at this level to non-negative */
1012 for( ii=ibot ; ii < itop ; ii++ ){
1013 ijk = ii+jj*nx+kk*nxy ;
1014 if( vmap[ijk] < -1 ) vmap[ijk] = (-vmap[ijk] - ntop) | lmask ;
1015 }}}
1016
1017 } /* end of loop over lev */
1018
1019 RETURN( vmap );
1020 }
1021
1022 /*-------------------------------------------------------------------------*/
1023 /*! Load the surface for this dataset from its file (if any).
1024 ---------------------------------------------------------------------------*/
1025
SUMA_load(THD_3dim_dataset * dset)1026 void SUMA_load( THD_3dim_dataset *dset )
1027 {
1028 int ks ; /* 14 Aug 2002: surface index */
1029
1030 ENTRY("SUMA_load") ;
1031
1032 if( !ISVALID_DSET(dset) ||
1033 dset->su_num == 0 ||
1034 dset->su_sname == NULL ) EXRETURN ;
1035
1036 /* 1st time in: allocate arrays to hold surface data */
1037
1038 if( dset->su_surf == NULL ){
1039 dset->su_surf = (SUMA_surface **) calloc(dset->su_num,sizeof(SUMA_surface *));
1040 dset->su_vmap = (int **) calloc(dset->su_num,sizeof(int *) );
1041 dset->su_vnlist = (SUMA_vnlist **) calloc(dset->su_num,sizeof(SUMA_vnlist *) );
1042 }
1043
1044 for( ks=0 ; ks < dset->su_num ; ks++ ){ /* loop over surfaces */
1045
1046 if( dset->su_surf[ks] != NULL ) continue ; /* already loaded */
1047
1048 dset->su_surf[ks] = SUMA_read_surface( dset->su_sname[ks] , dset ) ;
1049
1050 if( dset->su_surf[ks] == NULL ) continue ;
1051
1052 if( dset->su_vmap[ks] != NULL ) free(dset->su_vmap[ks]) ;
1053 #if 0
1054 dset->su_vmap[ks] = SUMA_map_dset_to_surf( dset->su_surf , dset ) ;
1055 #else
1056 dset->su_vmap[ks] = NULL ;
1057 #endif
1058
1059 if( dset->su_vnlist[ks] != NULL ){
1060 SUMA_destroy_vnlist( dset->su_vnlist[ks] ) ;
1061 dset->su_vnlist[ks] = NULL ;
1062 }
1063 #if 0
1064 dset->su_vnlist[ks] = SUMA_make_vnlist( dset->su_surf[ks] , dset ) ;
1065 #endif
1066
1067 }
1068
1069 EXRETURN ;
1070 }
1071
1072 /*--------------------------------------------------------------------------*/
1073 /*! Free the surface structs for this dataset (if any).
1074 ----------------------------------------------------------------------------*/
1075
SUMA_unload(THD_3dim_dataset * dset)1076 void SUMA_unload( THD_3dim_dataset *dset )
1077 {
1078 int ks ; /* 14 Aug 2002: surface index */
1079
1080 ENTRY("SUMA_unload") ;
1081
1082 if( !ISVALID_DSET(dset) ||
1083 dset->su_num == 0 ||
1084 dset->su_sname == NULL ||
1085 dset->su_surf == NULL ) EXRETURN ;
1086
1087 for( ks=0 ; ks < dset->su_num ; ks++ ){
1088
1089 if( dset->su_sname[ks] == NULL ||
1090 strstr(dset->su_sname[ks],"++LOCK++") != NULL ) continue ;
1091
1092 if( dset->su_surf[ks] != NULL ){
1093 SUMA_destroy_surface( dset->su_surf[ks] ) ; dset->su_surf[ks] = NULL ;
1094 }
1095
1096 if( dset->su_vmap[ks] != NULL ){
1097 free( dset->su_vmap[ks] ) ; dset->su_vmap[ks] = NULL ;
1098 }
1099
1100 if( dset->su_vnlist[ks] != NULL ){
1101 SUMA_destroy_vnlist( dset->su_vnlist[ks] ) ; dset->su_vnlist[ks] = NULL ;
1102 }
1103
1104 }
1105
1106 EXRETURN ;
1107 }
1108
1109 /*---------------------------------------------------------------------
1110 Read a <SureFit .../> line into a surface
1111 ag = already existing surface (nodes loaded into this)
1112 lbuf = line containing "<SureFit"
1113 dset = dataset for SureFit-to-DICOM coordinate conversion
1114 -----------------------------------------------------------------------*/
1115
SUMA_import_surefit(SUMA_surface * ag,char * lbuf,THD_3dim_dataset * dset)1116 void SUMA_import_surefit( SUMA_surface *ag, char *lbuf, THD_3dim_dataset *dset )
1117 {
1118 float xx[NBUF],yy[NBUF],zz[NBUF] ;
1119 int pp[NBUF] ;
1120 int nn , ii , idadd=0 ;
1121 FILE *sfp ;
1122 char sname[1024] , *cpt ;
1123 THD_fvec3 fv ;
1124
1125 ENTRY("SUMA_import_surefit") ;
1126
1127 /* scan input line for coord=sname, and extract into sname */
1128
1129 cpt = strstr(lbuf,"coord=") ;
1130 if( cpt == NULL ){
1131 fprintf(stderr,"** SUMA: Illegal SureFit: no coord=\n** %s\n",lbuf) ;
1132 EXRETURN ;
1133 }
1134 cpt += 6 ; /* skip coord= */
1135 if( *cpt == '\"' || *cpt == '\'' ) cpt++ ; /* skip quote */
1136 ii = sscanf(cpt,"%s",sname) ; /* get sname */
1137 if( ii == 0 ){
1138 fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
1139 EXRETURN ;
1140 }
1141 ii = strlen(sname) ;
1142 if( ii == 0 ){
1143 fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
1144 EXRETURN ;
1145 }
1146 if( sname[ii-1] == '\'' || sname[ii-1] == '\"' ) sname[ii-1] = '\0' ;
1147 if( strlen(sname) == 0 ){
1148 fprintf(stderr,"** SUMA: Illegal SureFit: bad coord=\n** %s\n",lbuf) ;
1149 EXRETURN ;
1150 }
1151
1152 /* add dataset directory name to start of sname? */
1153
1154 if( sname[0] != '/' ){
1155 char buf[1024] ;
1156 sprintf(buf,"%s%s",DSET_DIRNAME(dset),sname) ;
1157 strcpy(sname,buf) ;
1158 }
1159
1160 /* scan line for IDadd=value, and extract into idadd */
1161
1162 cpt = strstr(lbuf,"IDadd=") ;
1163 if( cpt != NULL ){
1164 cpt += 6 ;
1165 if( *cpt == '\"' || *cpt == '\'' ) cpt++ ;
1166 ii = sscanf(cpt,"%d",&idadd) ;
1167 if( ii == 0 || idadd < 0 ){
1168 fprintf(stderr,"** SUMA: Illegal SureFit: bad IDadd=\n** %s\n",lbuf) ;
1169 EXRETURN ;
1170 }
1171 }
1172
1173 /* open sname */
1174
1175 sfp = fopen( sname , "r" ) ;
1176 if( sfp == NULL ){
1177 fprintf(stderr,"** SUMA: Illegal SureFit: can't open file %s\n** %s\n",sname,lbuf) ;
1178 EXRETURN ;
1179 }
1180
1181 nn = 0 ;
1182
1183 while(1){
1184 cpt = afni_fgets(sname,1024,sfp) ; /* read a line */
1185 if( cpt == NULL ) break ; /* end of file */
1186
1187 if( strstr(sname,"BeginHeader") != NULL ){ /* skip SureFit header */
1188 do{
1189 cpt = afni_fgets(sname,1024,sfp) ; /* get next line */
1190 if( cpt == NULL ){ fclose(sfp); EXRETURN; } /* bad */
1191 } while( strstr(sname,"EndHeader") == NULL ) ;
1192 cpt = afni_fgets(sname,1024,sfp) ; /* 1 more line */
1193 if( cpt == NULL ){ fclose(sfp); EXRETURN; } /* bad */
1194 continue ; /* start over */
1195 }
1196
1197 ii = sscanf(sname,"%d%f%f%f",pp+nn,xx+nn,yy+nn,zz+nn) ;
1198 if( ii < 4 ) continue ; /* skip this line; it's bad */
1199
1200 /* process value to AFNI-ize it from SureFit */
1201
1202 pp[nn] += idadd ;
1203 LOAD_FVEC3(fv,xx[nn],yy[nn],zz[nn]) ;
1204 fv = THD_surefit_to_dicomm( dset , fv ) ;
1205 UNLOAD_FVEC3(fv,xx[nn],yy[nn],zz[nn]) ;
1206
1207 nn++ ;
1208 if( nn == NBUF ){
1209 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
1210 nn = 0 ;
1211 }
1212 } /* end of loop over input lines */
1213
1214 fclose(sfp) ;
1215
1216 if( nn > 0 ){
1217 SUMA_add_nodes_ixyz( ag,nn , pp,xx,yy,zz ) ;
1218 }
1219
1220 EXRETURN ;
1221 }
1222
1223 /*-------------------------------------------------------------------------*/
1224 /*! Load the surface name into the dataset struct (derived by replacing
1225 .HEAD with .SURF).
1226 ---------------------------------------------------------------------------*/
1227
SUMA_get_surfname(THD_3dim_dataset * dset)1228 void SUMA_get_surfname( THD_3dim_dataset *dset )
1229 {
1230 char *snam ;
1231 int ii , ks ;
1232
1233 ENTRY("THD_get_surfname") ;
1234
1235 if( !ISVALID_DSET(dset) || dset->su_num > 0 ) EXRETURN ;
1236
1237 snam = strdup( DSET_HEADNAME(dset) ) ;
1238 ii = strlen(snam) ;
1239 if( ii > 5 ){
1240 strcpy(snam+ii-4,"SURF") ;
1241 if( THD_is_file(snam) ){
1242 dset->su_num = 1 ;
1243 dset->su_sname = (char **) malloc(sizeof(char *)) ;
1244 dset->su_sname[0] = snam;
1245 EXRETURN;
1246 }
1247 }
1248 free(snam) ; EXRETURN ; /* .SURF file does not exist */
1249 }
1250
1251
1252 #endif /* ALLOW_SURFACE_FILES */
1253