1 #include "mrilib.h"
2 #include "thd.h"
3
4 /*******************************************************************/
5 /********** 21 Mar 2003: Read a 3D file as an AFNI dataset *********/
6 /*******************************************************************/
7
8 /*-----------------------------------------------------------------*/
9 /*! Open a NIML 3D file as an AFNI dataset. Each column is a
10 separate sub-brick. Dataset returned is empty (no data).
11
12 Broken into modules. 31 May 2006 [rickr]
13 -------------------------------------------------------------------*/
14
THD_open_3D(char * pathname)15 THD_3dim_dataset * THD_open_3D( char *pathname )
16 {
17 THD_3dim_dataset *dset=NULL ; /* use locals only for tracing */
18 NI_element *nel ;
19
20 ENTRY("THD_open_3D") ;
21
22 nel = read_niml_file(pathname, 0) ; /* do not get data */
23 if( nel ) dset = THD_niml_3D_to_dataset(nel, pathname) ;
24
25 RETURN(dset) ;
26 }
27
28 /*-----------------------------------------------------------------*/
29 /*! Convert a NIML 3D element (no data) to an AFNI dataset.
30 This was moved from THD_open_3D(). 31 May 2006 [rickr]
31 -------------------------------------------------------------------*/
THD_niml_3D_to_dataset(NI_element * nel,char * pathname)32 THD_3dim_dataset * THD_niml_3D_to_dataset( NI_element * nel, char * pathname )
33 {
34 THD_3dim_dataset *dset=NULL ;
35 char prefix[1024] , *ppp ;
36 THD_ivec3 nxyz , orixyz ;
37 THD_fvec3 dxyz , orgxyz ;
38 int nvals , ii ;
39
40 ENTRY("THD_3D_niml_to_dataset") ;
41
42 /*-- check data element for reasonability --*/
43
44 if( nel == NULL || /* bad read */
45 NI_element_type(nel) != NI_ELEMENT_TYPE || /* bad element */
46 nel->vec_num <= 0 || /* no data */
47 nel->vec_len <= 0 || /* no data */
48 strcmp(nel->name,"AFNI_3D_dataset") != 0 || /* incorrect data */
49 nel->vec_rank > 3 ){ /* weird header */
50
51 fprintf(stderr,"** Can't read 3D head from %s\n",pathname) ;
52 NI_free_element(nel) ; RETURN(NULL) ;
53 }
54
55 STATUS("checking header") ;
56
57 /*-- check column types to make sure they are all numeric --*/
58 /* [AFNI doesn't like String or compound type datasets] */
59
60 for( ii=0 ; ii < nel->vec_num ; ii++ ){
61
62 if( !NI_IS_NUMERIC_TYPE(nel->vec_typ[ii]) ){
63 fprintf(stderr,"** 3D file %s isn't numeric!\n",pathname) ;
64 NI_free_element(nel) ; RETURN(NULL) ;
65 }
66 }
67
68 /*-- now have good data element ==> make a dataset --*/
69
70 STATUS("making dataset") ;
71
72 dset = EDIT_empty_copy(NULL) ; /* default dataset */
73
74 /* set prefix from input filename */
75
76 STATUS("setting prefix") ;
77
78 ppp = THD_trailname(pathname,0) ; /* strip directory */
79 NI_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
80
81 /* set grid sizes from element header */
82
83 STATUS("setting grid sizes") ;
84
85 nxyz.ijk[0] = nel->vec_len ; nxyz.ijk[1] = nxyz.ijk[2] = 1 ;
86 if( nel->vec_axis_len != NULL ){
87 if( nel->vec_rank >= 1 ) nxyz.ijk[0] = nel->vec_axis_len[0] ;
88 if( nel->vec_rank >= 2 ) nxyz.ijk[1] = nel->vec_axis_len[1] ;
89 if( nel->vec_rank >= 3 ) nxyz.ijk[2] = nel->vec_axis_len[2] ;
90 }
91
92 /* set grid spacings */
93
94 STATUS("setting grid spacings") ;
95
96 dxyz.xyz[0] = dxyz.xyz[1] = dxyz.xyz[2] = 1.0 ;
97 if( nel->vec_axis_delta != NULL ){
98 if( nel->vec_rank >= 1) dxyz.xyz[0] = nel->vec_axis_delta[0] ;
99 if( nel->vec_rank >= 2) dxyz.xyz[1] = nel->vec_axis_delta[1] ;
100 if( nel->vec_rank >= 3) dxyz.xyz[2] = nel->vec_axis_delta[2] ;
101 }
102
103 /* set grid origins */
104
105 STATUS("setting grid origins") ;
106
107 orgxyz.xyz[0] = orgxyz.xyz[1] = orgxyz.xyz[2] = 0.0 ;
108 if( nel->vec_axis_origin != NULL ){
109 if( nel->vec_rank >= 1) orgxyz.xyz[0] = nel->vec_axis_origin[0] ;
110 if( nel->vec_rank >= 2) orgxyz.xyz[1] = nel->vec_axis_origin[1] ;
111 if( nel->vec_rank >= 3) orgxyz.xyz[2] = nel->vec_axis_origin[2] ;
112 }
113
114 /* set grid orientations (default is RAI) */
115
116 STATUS("setting grid orientation") ;
117
118 { char orcx='R', orcy='A', orcz='I' ;
119 int oxx,oyy,ozz ;
120 if( nel->vec_rank == 3 && nel->vec_axis_label != NULL ){
121 orcx = toupper(nel->vec_axis_label[0][0]) ;
122 orcy = toupper(nel->vec_axis_label[1][0]) ;
123 orcz = toupper(nel->vec_axis_label[2][0]) ;
124 }
125 oxx = ORCODE(orcx); oyy = ORCODE(orcy); ozz = ORCODE(orcz);
126 if( !OR3OK(oxx,oyy,ozz) ){
127 oxx = ORI_R2L_TYPE; oyy = ORI_A2P_TYPE; ozz = ORI_I2S_TYPE;
128 }
129
130 orixyz.ijk[0] = oxx ; orixyz.ijk[1] = oyy ; orixyz.ijk[2] = ozz ;
131 }
132
133 /* number of sub-bricks (one per vector/column) */
134
135 nvals = nel->vec_num ;
136
137 /* set idcode from element, or take random default one */
138
139 STATUS("setting idcode") ;
140
141 ppp = NI_get_attribute( nel , "self_idcode" ) ;
142 if( ppp == NULL )
143 ppp = NI_get_attribute( nel , "ni_idcode" ) ;
144 if( ppp != NULL && *ppp != '\0' ){
145 NI_strncpy( dset->idcode.str , ppp , MCW_IDSIZE ) ;
146 } else {
147 MCW_hash_idcode( pathname , dset ) ; /* 06 May 2005 */
148 dset->idcode.str[0] = 'A' ; /* overwrite 1st 3 bytes of idcode */
149 dset->idcode.str[1] = '3' ;
150 dset->idcode.str[2] = 'D' ;
151 }
152
153 /*-- now modify the default dataset --*/
154
155 STATUS("Editing dataset") ;
156
157 EDIT_dset_items( dset ,
158 ADN_prefix , prefix ,
159 ADN_datum_array , nel->vec_typ ,
160 ADN_nxyz , nxyz ,
161 ADN_xyzdel , dxyz ,
162 ADN_xyzorg , orgxyz ,
163 ADN_xyzorient , orixyz ,
164 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
165 ADN_nvals , nvals ,
166 ADN_type , HEAD_ANAT_TYPE ,
167 ADN_func_type , ANAT_BUCK_TYPE ,
168 ADN_none ) ;
169
170 dset->dblk->diskptr->storage_mode = STORAGE_BY_3D ;
171 NI_strncpy( dset->dblk->diskptr->brick_name , pathname , THD_MAX_NAME ) ;
172
173 /*-- time axis? [03 Jun 2005] --*/
174
175 ppp = NI_get_attribute( nel , "ni_timestep" ) ;
176 if( ppp != NULL && nvals > 1 ){
177 float dt = strtod(ppp,NULL) ; if( dt <= 0.0 ) dt = 1.0 ;
178 if (dt > 360) {
179 dt *= 0.001;
180 fprintf(stderr,"Warning: ni_timestep appears incorrecly set in msec.\n"
181 "Reducing dt by a factor of 1000 to %fsec.\n", dt);
182 }
183 EDIT_dset_items( dset ,
184 ADN_func_type , ANAT_EPI_TYPE ,
185 ADN_ntt , nvals ,
186 ADN_ttdel , dt ,
187 ADN_tunits , UNITS_SEC_TYPE ,
188 ADN_none ) ;
189 }
190
191 STATUS("checking for statistics") ;
192
193 /*-- see if we have any statistics bricks --*/
194
195 ppp = NI_get_attribute( nel , "ni_stat" ) ;
196 if( ppp != NULL ){
197 NI_str_array *sar = NI_decode_string_list( ppp , ";" ) ;
198 if( sar != NULL ){
199 int itop=MIN(sar->num,nvals) , jj,ll ;
200 char *dnam , qnam[64] ;
201 for( ii=0 ; ii < itop ; ii++ ){
202 if( strcmp(sar->str[ii],"none") == 0 ) continue ;
203 for( jj=NI_STAT_FIRSTCODE ; jj <= NI_STAT_LASTCODE ; jj++ ){
204 dnam = NI_stat_distname(jj) ;
205 strcpy(qnam,dnam); strcat(qnam,"("); ll = strlen(qnam);
206 if( strncmp(sar->str[ii],qnam,ll) == 0 ) break ;
207 }
208 if( jj >= AFNI_FIRST_STATCODE && jj <= AFNI_LAST_STATCODE ){
209 float parm[4]={1,1,1,1} ; int np,kk,mm , sp ;
210 np = NI_stat_numparam(jj) ; sp = ll ;
211 for( kk=0 ; kk < np ; kk++ ){
212 mm = 0 ; sscanf(sar->str[ii]+sp,"%f%n",parm+kk,&mm) ; sp += mm+1 ;
213 }
214 EDIT_STATAUX4( dset , ii , jj , parm[0],parm[1],parm[2],parm[3] ) ;
215 }
216 }
217 NI_delete_str_array(sar) ;
218 }
219 }
220
221 /*-- purge the NIML data element and return the new dataset --*/
222
223 STATUS("freeing element") ;
224
225 NI_free_element( nel ) ;
226 RETURN(dset) ;
227 }
228
229 /*-----------------------------------------------------------------*/
230 /*! Mark a dataset's storage mode.
231 -------------------------------------------------------------------*/
232
THD_set_storage_mode(THD_3dim_dataset * dset,int mm)233 void THD_set_storage_mode( THD_3dim_dataset *dset , int mm )
234 {
235 if( !ISVALID_DSET(dset) ||
236 mm < 0 ||
237 mm > LAST_STORAGE_MODE ||
238 dset->dblk == NULL ||
239 dset->dblk->diskptr == NULL ) return ;
240
241 dset->dblk->diskptr->storage_mode = mm ;
242 }
243
244 /*-----------------------------------------------------------------*/
245 /*! Load a 3D dataset's data into memory.
246 Called from THD_load_datablock() in thd_loaddblk.c.
247 -------------------------------------------------------------------*/
248
THD_load_3D(THD_datablock * dblk)249 void THD_load_3D( THD_datablock *dblk )
250 {
251 THD_diskptr *dkptr ;
252 NI_element *nel ;
253 NI_stream ns ;
254 int nxyz , iv,nv ;
255 char *ppp ;
256
257 ENTRY("THD_load_3D") ;
258
259 /*-- open input [these errors should never occur] --*/
260
261 if( !ISVALID_DATABLOCK(dblk) ||
262 dblk->diskptr->storage_mode != STORAGE_BY_3D ||
263 dblk->brick == NULL ) EXRETURN ;
264
265 dkptr = dblk->diskptr ;
266 nxyz = dkptr->dimsizes[0] * dkptr->dimsizes[1] * dkptr->dimsizes[2] ;
267 nv = dkptr->nvals ;
268
269 if( nxyz*nv > 1000000 ) fprintf(stderr,"++ Reading %s\n",dkptr->brick_name) ;
270
271 ppp = (char*)calloc( sizeof(char) , strlen(dkptr->brick_name)+16 ) ;
272
273 strcpy(ppp,"file:") ; strcat(ppp,dkptr->brick_name) ;
274 ns = NI_stream_open( ppp , "r" ) ; free(ppp) ;
275 if( ns == NULL ) EXRETURN ;
276
277 NI_skip_procins(1) ;
278 nel = NI_read_element(ns,333); NI_stream_close(ns);
279 NI_skip_procins(0) ;
280 if( nel == NULL ) EXRETURN ;
281
282 /*-- allocate space for data --*/
283
284 if( nxyz != nel->vec_len || nv != nel->vec_num ){
285 fprintf(stderr,"THD_load_3D(%s): nxyz or nv mismatch!\n",dkptr->brick_name) ;
286 fprintf(stderr," expect nxyz=%d; got %d\n",nxyz, nel->vec_len);
287 fprintf(stderr," expect nv =%d; got %d\n",nv , nel->vec_num);
288 NI_free_element(nel) ; EXRETURN ;
289 }
290
291 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
292
293 /** malloc space for each brick separately,
294 and copy data from element into place **/
295
296 for( iv=0 ; iv < nv ; iv++ ){
297 if( DBLK_ARRAY(dblk,iv) == NULL ){ /* needs data */
298 ppp = AFMALL(char, DBLK_BRICK_BYTES(dblk,iv) ); /* make space */
299 if( ppp == NULL ) break ; /* bad bad bad */
300 mri_fix_data_pointer( ppp, DBLK_BRICK(dblk,iv) ) ;
301 memcpy( ppp, nel->vec[iv], DBLK_BRICK_BYTES(dblk,iv) ) ;
302 NI_free(nel->vec[iv]) ; nel->vec[iv] = NULL ;
303 }
304 }
305
306 NI_free_element(nel) ;
307
308 /* if malloc failed, then delete all bricks */
309
310 if( iv < nv ){
311 fprintf(stderr,"\n** malloc failed for 3D dataset input!\n");
312 for( iv=0 ; iv < nv ; iv++ ){
313 if( DBLK_ARRAY(dblk,iv) != NULL ){
314 free( DBLK_ARRAY(dblk,iv) ) ;
315 mri_fix_data_pointer( NULL, DBLK_BRICK(dblk,iv) ) ;
316 }
317 }
318 }
319
320 EXRETURN ;
321 }
322
323 /*------------------------------------------------------------------*/
324 /*! Write a dataset to disk as a 3D file.
325 Called from THD_write_3dim_dataset().
326 This is kind of cheating, since we just call THD_write_1D()
327 instead.
328 --------------------------------------------------------------------*/
329
THD_write_3D(char * sname,char * pname,THD_3dim_dataset * dset)330 void THD_write_3D( char *sname, char *pname , THD_3dim_dataset *dset )
331 {
332 THD_write_1D( sname,pname,dset ) ;
333 }
334