1 #include "mrilib.h"
2 #include "thd.h"
3
4 /*******************************************************************/
5 /* Stuff to read CTF MRI and SAM datasets from the NIH MEG system. */
6 /* 04 Dec 2002: first cut by RWCox. */
7 /*******************************************************************/
8
9 /****************************************************************/
10 /*********** header structs for CTF MRI version 2.0 file */
11 /****************************************************************/
12
13 enum { coronal=0, sagittal, axial };
14
15 #define LEFT_ON_LEFT 0
16 #define LEFT_ON_RIGHT 1
17
18 enum { Modality_MRI=0, Modality_CT, Modality_PET, Modality_SPECT, Modality_OTHER };
19
20 #define VERSION_1_STR "CTF_MRI_FORMAT VER 1.0"
21 #define VERSION_2_STR "CTF_MRI_FORMAT VER 2.0"
22 #define VERSION_21_STR "CTF_MRI_FORMAT VER 2.1"
23 #define VERSION_22_STR "CTF_MRI_FORMAT VER 2.2"
24
25 typedef struct HeadModel_Info {
26 short Nasion_Sag; /* fiduciary point voxel locations */
27 short Nasion_Cor; /* Sag = sagittal direction */
28 short Nasion_Axi; /* Cor = coronal direction */
29 short LeftEar_Sag; /* Axi = axial direction */
30 short LeftEar_Cor;
31 short LeftEar_Axi;
32 short RightEar_Sag;
33 short RightEar_Cor;
34 short RightEar_Axi;
35 float defaultSphereX; /* default sphere parameters in mm */
36 float defaultSphereY; /* (in head based coordinate system */
37 float defaultSphereZ;
38 float defaultSphereRadius;
39 } HeadModel_Info;
40
41 /* this struct isn't used in AFNI */
42 typedef struct Image_Info { /* scan and/or sequence parameters */
43 short modality; /* 0=MRI, 1=CT, 2=PET, 3=SPECT, 4=OTHER */
44 char manufacturerName[64];
45 char instituteName[64];
46 char patientID[32];
47 char dateAndTime[32];
48 char scanType[32];
49 char contrastAgent[32];
50 char imagedNucleus[32];
51 float Frequency;
52 float FieldStrength;
53 float EchoTime;
54 float RepetitionTime;
55 float InversionTime;
56 float FlipAngle;
57 short NoExcitations;
58 short NoAcquisitions;
59 char commentString[256];
60 char forFutureUse[64];
61 } Image_Info;
62
63 /* the header for the .mri files */
64 typedef struct Version_2_Header {
65 char identifierString[32]; /* "CTF_MRI_FORMAT VER 2.x" */
66 short imageSize; /* always = 256 */
67 short dataSize; /* 1 = 8 bit data, 2 = 16 bit data */
68 short clippingRange; /* max. integer value in data */
69 short imageOrientation; /* 0 = left on left, 1 = left on right */
70 float mmPerPixel_sagittal; /* voxel dimensions in mm */
71 float mmPerPixel_coronal; /* voxel dimensions in mm */
72 float mmPerPixel_axial; /* voxel dimensions in mm */
73 HeadModel_Info headModel; /* structure defined above (34 bytes) */
74 Image_Info imageInfo; /* structure defined above (638 bytes) */
75 float headOrigin_sagittal; /* voxel location of head origin */
76 float headOrigin_coronal; /* voxel location of head origin */
77 float headOrigin_axial; /* voxel location of head origin */
78
79 /* Euler angles to align MR to head */
80 /* coordinate system (angles in degrees!) */
81 float rotate_coronal; /* 1. rotate in coronal plane by this angle */
82 float rotate_sagittal; /* 2. rotate in sagittal plane by this angle */
83 float rotate_axial; /* 3. rotate in axial plane by this angle */
84
85 short orthogonalFlag; /* 1 if image is orthogonalized to head frame */
86 short interpolatedFlag; /* 1 if slices interpolated during conversion */
87 float originalSliceThickness;
88 float transformMatrix[4][4]; /* 4x4 transformation matrix (head to mri) */
89 unsigned char unused[202]; /* pad header to 1028 bytes */
90 } Version_2_Header;
91
92 /*---------------------------------------------------------------*/
93 /*! Swap the 4 bytes pointed to by ppp: abcd -> dcba. */
94
swap_4(void * ppp)95 static void swap_4(void *ppp)
96 {
97 unsigned char *pntr = (unsigned char *) ppp ;
98 unsigned char b0, b1, b2, b3;
99
100 b0 = *pntr; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
101 *pntr = b3; *(pntr+1) = b2; *(pntr+2) = b1; *(pntr+3) = b0;
102 }
103
104 /*---------------------------------------------------------------*/
105 /*! Swap the 8 bytes pointed to by ppp: abcdefgh -> hgfedcba. */
106
swap_8(void * ppp)107 static void swap_8(void *ppp)
108 {
109 unsigned char *pntr = (unsigned char *) ppp ;
110 unsigned char b0, b1, b2, b3;
111 unsigned char b4, b5, b6, b7;
112
113 b0 = *pntr ; b1 = *(pntr+1); b2 = *(pntr+2); b3 = *(pntr+3);
114 b4 = *(pntr+4); b5 = *(pntr+5); b6 = *(pntr+6); b7 = *(pntr+7);
115
116 *pntr = b7; *(pntr+1) = b6; *(pntr+2) = b5; *(pntr+3) = b4;
117 *(pntr+4) = b3; *(pntr+5) = b2; *(pntr+6) = b1; *(pntr+7) = b0;
118 }
119
120 /*---------------------------------------------------------------*/
121 /*! Swap the 2 bytes pointed to by ppp: ab -> ba. */
122
swap_2(void * ppp)123 static void swap_2(void *ppp)
124 {
125 unsigned char *pntr = (unsigned char *) ppp ;
126 unsigned char b0, b1;
127
128 b0 = *pntr; b1 = *(pntr+1);
129 *pntr = b1; *(pntr+1) = b0;
130 }
131
132 /*-------------------------*/
133 /*! Macro for bad return. */
134
135 #undef BADBAD
136 #define BADBAD(s) \
137 do{ fprintf(stderr,"** THD_open_ctfmri(%s): %s\n",fname,s); \
138 RETURN(NULL); \
139 } while(0)
140
141 /*-----------------------------------------------------------------*/
142 /*! Function to count slices like CTF does. */
143
CTF_count(double start,double end,double delta)144 int CTF_count( double start, double end , double delta )
145 {
146 int nn=0 ; double cc ;
147 for( cc=start ; cc <= (end+1.0e-6) ; cc += delta ) nn++ ;
148 return nn ;
149 }
150
151 /*-----------------------------------------------------------------*/
152 /*! Open a CTF .mri file as an unpopulated AFNI dataset.
153 It will be populated later, in THD_load_ctfmri().
154 -------------------------------------------------------------------*/
155
THD_open_ctfmri(char * fname)156 THD_3dim_dataset * THD_open_ctfmri( char *fname )
157 {
158 FILE *fp ;
159 Version_2_Header hh ;
160 int ii,nn , swap ;
161 THD_3dim_dataset *dset=NULL ;
162 char prefix[THD_MAX_PREFIX] , *ppp , tname[12] , ori[4] ;
163 THD_ivec3 nxyz , orixyz ;
164 THD_fvec3 dxyz , orgxyz ;
165 int iview ;
166 int ngood , length , datum_type=0 , datum_len , oxx,oyy,ozz ;
167 int nx,ny,nz ;
168 float dx,dy,dz=0.0, xorg,yorg,zorg=0.0;
169
170 ENTRY("THD_open_ctfmri") ;
171
172 /* open input file */
173
174 if( fname == NULL || *fname == '\0' ) BADBAD("bad input filename");
175 fp = fopen( fname , "rb" ) ;
176 if( fp == NULL ) BADBAD("can't open input file");
177
178 /* read 1028 byte header */
179
180 nn = fread( &hh , sizeof(hh) , 1 , fp ) ;
181 fclose(fp) ;
182 if( nn <= 0 ) BADBAD("can't read input file");
183
184 /* check if header string matches what we want */
185
186 hh.identifierString[31] = '\0' ; /* make sure is terminated */
187 if( strcmp(hh.identifierString,VERSION_22_STR) != 0 ) BADBAD("bad version string");
188 if( hh.imageSize == 0 ) BADBAD("bad imageSize") ;
189
190 /* determine if must swap header */
191
192 swap = (hh.imageSize != 256) ;
193
194 if( swap ){ /* swap bytes in various header fields */
195 swap_2(&hh.imageSize ) ;
196 swap_2(&hh.dataSize ) ;
197 swap_2(&hh.clippingRange ) ;
198 swap_2(&hh.imageOrientation ) ;
199 swap_4(&hh.mmPerPixel_sagittal ) ;
200 swap_4(&hh.mmPerPixel_coronal ) ;
201 swap_4(&hh.mmPerPixel_axial ) ;
202 swap_4(&hh.headOrigin_sagittal ) ;
203 swap_4(&hh.headOrigin_coronal ) ;
204 swap_4(&hh.headOrigin_axial ) ;
205 swap_4(&hh.rotate_coronal ) ;
206 swap_4(&hh.rotate_sagittal ) ;
207 swap_4(&hh.rotate_axial ) ;
208 swap_2(&hh.orthogonalFlag ) ;
209 swap_2(&hh.interpolatedFlag ) ;
210 swap_4(&hh.originalSliceThickness ) ;
211 swap_2(&hh.headModel.Nasion_Sag ) ;
212 swap_2(&hh.headModel.Nasion_Cor ) ;
213 swap_2(&hh.headModel.Nasion_Axi ) ;
214 swap_2(&hh.headModel.LeftEar_Sag ) ;
215 swap_2(&hh.headModel.LeftEar_Cor ) ;
216 swap_2(&hh.headModel.LeftEar_Axi ) ;
217 swap_2(&hh.headModel.RightEar_Sag ) ;
218 swap_2(&hh.headModel.RightEar_Cor ) ;
219 swap_2(&hh.headModel.RightEar_Axi ) ;
220
221 swap_4(&hh.transformMatrix[0][0] ) ; /* this stuff not used yet */
222 swap_4(&hh.transformMatrix[0][1] ) ;
223 swap_4(&hh.transformMatrix[0][2] ) ;
224 swap_4(&hh.transformMatrix[0][3] ) ;
225 swap_4(&hh.transformMatrix[1][0] ) ;
226 swap_4(&hh.transformMatrix[1][1] ) ;
227 swap_4(&hh.transformMatrix[1][2] ) ;
228 swap_4(&hh.transformMatrix[1][3] ) ;
229 swap_4(&hh.transformMatrix[2][0] ) ;
230 swap_4(&hh.transformMatrix[2][1] ) ;
231 swap_4(&hh.transformMatrix[2][2] ) ;
232 swap_4(&hh.transformMatrix[2][3] ) ;
233 swap_4(&hh.transformMatrix[3][0] ) ;
234 swap_4(&hh.transformMatrix[3][1] ) ;
235 swap_4(&hh.transformMatrix[3][2] ) ;
236 swap_4(&hh.transformMatrix[3][3] ) ;
237 }
238
239 /* simple checks on header stuff */
240
241 if( hh.imageSize != 256 ||
242 hh.dataSize < 1 ||
243 hh.dataSize > 2 ||
244 hh.mmPerPixel_sagittal <= 0.0 ||
245 hh.mmPerPixel_coronal <= 0.0 ||
246 hh.mmPerPixel_axial <= 0.0 ) BADBAD("bad header data") ;
247
248 /*- 16 Mar 2005: instead of complaining about negative Origins,
249 just reset them to something semi-reasonable;
250 I just get by with a little help from my friends
251 - Zuxiang Li in this case. -*/
252
253 if( hh.headOrigin_sagittal <= 0.0 ) hh.headOrigin_sagittal = hh.imageSize/2.;
254 if( hh.headOrigin_coronal <= 0.0 ) hh.headOrigin_coronal = hh.imageSize/2.;
255 if( hh.headOrigin_axial <= 0.0 ) hh.headOrigin_axial = hh.imageSize/2.;
256
257 /* debugging code to print header information */
258 #if 0
259 printf("\n") ;
260 printf("*** CTF MRI filename = %s\n",fname ) ;
261 printf("identifierString = %s\n",hh.identifierString ) ;
262 printf("imageSize = %d\n",hh.imageSize ) ;
263 printf("dataSize = %d\n",hh.dataSize ) ;
264 printf("clippingRange = %d\n",hh.clippingRange ) ;
265 printf("imageOrientation = %d\n",hh.imageOrientation ) ;
266 printf("mmPerPixel_sagittal = %f\n",hh.mmPerPixel_sagittal) ;
267 printf("mmPerPixel_coronal = %f\n",hh.mmPerPixel_coronal ) ;
268 printf("mmPerPixel_axial = %f\n",hh.mmPerPixel_axial ) ;
269 printf("headOrigin_sagittal = %f\n",hh.headOrigin_sagittal) ;
270 printf("headOrigin_coronal = %f\n",hh.headOrigin_coronal ) ;
271 printf("headOrigin_axial = %f\n",hh.headOrigin_axial ) ;
272 printf("rotate_coronal = %f\n",hh.rotate_coronal ) ;
273 printf("rotate_sagittal = %f\n",hh.rotate_sagittal ) ;
274 printf("rotate_axial = %f\n",hh.rotate_axial ) ;
275 printf("orthogonalFlag = %d\n",hh.orthogonalFlag ) ;
276 printf("interpolatedFlag = %d\n",hh.interpolatedFlag ) ;
277 printf("originalSliceThickness = %f\n",hh.originalSliceThickness ) ;
278 printf("\n") ;
279 printf("headModel.Nasion_Sag = %d\n",hh.headModel.Nasion_Sag ) ;
280 printf("headModel.Nasion_Cor = %d\n",hh.headModel.Nasion_Cor ) ;
281 printf("headModel.Nasion_Axi = %d\n",hh.headModel.Nasion_Axi ) ;
282 printf("headModel.LeftEar_Sag = %d\n",hh.headModel.LeftEar_Sag ) ;
283 printf("headModel.LeftEar_Cor = %d\n",hh.headModel.LeftEar_Cor ) ;
284 printf("headModel.LeftEar_Axi = %d\n",hh.headModel.LeftEar_Axi ) ;
285 printf("headModel.RightEar_Sag = %d\n",hh.headModel.RightEar_Sag) ;
286 printf("headModel.RightEar_Cor = %d\n",hh.headModel.RightEar_Cor) ;
287 printf("headModel.RightEar_Axi = %d\n",hh.headModel.RightEar_Axi) ;
288 printf("\n") ;
289 printf("transformMatrix:\n"
290 " [ %9.4f %9.4f %9.4f %9.4f ]\n"
291 " [ %9.4f %9.4f %9.4f %9.4f ]\n"
292 " [ %9.4f %9.4f %9.4f %9.4f ]\n"
293 " [ %9.4f %9.4f %9.4f %9.4f ]\n" ,
294 hh.transformMatrix[0][0] , hh.transformMatrix[0][1] ,
295 hh.transformMatrix[0][2] , hh.transformMatrix[0][3] ,
296 hh.transformMatrix[1][0] , hh.transformMatrix[1][1] ,
297 hh.transformMatrix[1][2] , hh.transformMatrix[1][3] ,
298 hh.transformMatrix[2][0] , hh.transformMatrix[2][1] ,
299 hh.transformMatrix[2][2] , hh.transformMatrix[2][3] ,
300 hh.transformMatrix[3][0] , hh.transformMatrix[3][1] ,
301 hh.transformMatrix[3][2] , hh.transformMatrix[3][3] ) ;
302 #endif
303
304 /* determine if file is big enough to hold all data it claims */
305
306 nn = THD_filesize(fname) ;
307 if( nn < hh.dataSize*hh.imageSize*hh.imageSize*hh.imageSize )
308 BADBAD("input file too small") ;
309
310 /*** from here, a lot of code is adapted from thd_analyzeread.c ***/
311
312 datum_len = hh.dataSize ;
313 switch( datum_len ){ /* the only 2 cases */
314 case 1: datum_type = MRI_byte ; break ;
315 case 2: datum_type = MRI_short; break ;
316 }
317 nx = ny = nz = hh.imageSize ; /* volumes are cubes! */
318
319 /* set orientation:
320 for now, assume (based on 1 sample) that data is stored in ASL or ASR order */
321
322 ori[0] = 'A' ; /* x is A-P */
323 ori[1] = 'S' ; /* y is S-I */
324
325 /* determine if z is L-R or R-L from position of markers */
326
327 ori[2] = (hh.headModel.LeftEar_Sag <= hh.headModel.RightEar_Sag) ? 'L' : 'R' ;
328
329 oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
330 if( !OR3OK(oxx,oyy,ozz) ){
331 oxx = ORI_A2P_TYPE; oyy = ORI_S2I_TYPE; ozz = ORI_L2R_TYPE; /** ASL? **/
332 }
333
334 /* now set grid size, keeping in mind that
335 A-P is positive and P-A is negative,
336 R-L is positive and L-R is negative,
337 I-S is positive and S-I is negative. */
338
339 switch( ori[0] ){
340 case 'A': dx = hh.mmPerPixel_coronal ; xorg = hh.headOrigin_coronal ; break ;
341 case 'P': dx = -hh.mmPerPixel_coronal ; xorg = hh.headOrigin_coronal ; break ;
342 case 'R': dx = hh.mmPerPixel_sagittal; xorg = hh.headOrigin_sagittal; break ;
343 case 'L': dx = -hh.mmPerPixel_sagittal; xorg = hh.headOrigin_sagittal; break ;
344 case 'I': dx = hh.mmPerPixel_axial ; xorg = hh.headOrigin_axial ; break ;
345 case 'S': dx = -hh.mmPerPixel_axial ; xorg = hh.headOrigin_axial ; break ;
346 }
347 switch( ori[1] ){
348 case 'A': dy = hh.mmPerPixel_coronal ; yorg = hh.headOrigin_coronal ; break ;
349 case 'P': dy = -hh.mmPerPixel_coronal ; yorg = hh.headOrigin_coronal ; break ;
350 case 'R': dy = hh.mmPerPixel_sagittal; yorg = hh.headOrigin_sagittal; break ;
351 case 'L': dy = -hh.mmPerPixel_sagittal; yorg = hh.headOrigin_sagittal; break ;
352 case 'I': dy = hh.mmPerPixel_axial ; yorg = hh.headOrigin_axial ; break ;
353 case 'S': dy = -hh.mmPerPixel_axial ; yorg = hh.headOrigin_axial ; break ;
354 }
355 switch( ori[2] ){
356 case 'A': dz = hh.mmPerPixel_coronal ; zorg = hh.headOrigin_coronal ; break ;
357 case 'P': dz = -hh.mmPerPixel_coronal ; zorg = hh.headOrigin_coronal ; break ;
358 case 'R': dz = hh.mmPerPixel_sagittal; zorg = hh.headOrigin_sagittal; break ;
359 case 'L': dz = -hh.mmPerPixel_sagittal; zorg = hh.headOrigin_sagittal; break ;
360 case 'I': dz = hh.mmPerPixel_axial ; zorg = hh.headOrigin_axial ; break ;
361 case 'S': dz = -hh.mmPerPixel_axial ; zorg = hh.headOrigin_axial ; break ;
362 }
363
364 /* At this point, (xorg,yorg,zorg) are voxel indices;
365 now, translate them into shifts such that if voxel
366 index ii is the location of x=0, then xorg+ii*dx=0. */
367
368 xorg = -dx*xorg ; yorg = -dy*yorg ; zorg = -dz*zorg ;
369
370 /*-- make a dataset --*/
371
372 dset = EDIT_empty_copy(NULL) ;
373
374 dset->idcode.str[0] = 'C' ; /* overwrite 1st 3 bytes */
375 dset->idcode.str[1] = 'T' ;
376 dset->idcode.str[2] = 'F' ;
377
378 MCW_hash_idcode( fname , dset ) ; /* 06 May 2005 */
379
380 ppp = THD_trailname(fname,0) ; /* strip directory */
381 MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
382
383 nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ; /* setup axes lengths and voxel sizes */
384 nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
385 nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
386
387 orixyz.ijk[0] = oxx ; orgxyz.xyz[0] = xorg ;
388 orixyz.ijk[1] = oyy ; orgxyz.xyz[1] = yorg ;
389 orixyz.ijk[2] = ozz ; orgxyz.xyz[2] = zorg ;
390
391 iview = VIEW_ORIGINAL_TYPE ;
392
393 /*-- actually send the values above into the dataset header --*/
394
395 EDIT_dset_items( dset ,
396 ADN_prefix , prefix ,
397 ADN_datum_all , datum_type ,
398 ADN_nxyz , nxyz ,
399 ADN_xyzdel , dxyz ,
400 ADN_xyzorg , orgxyz ,
401 ADN_xyzorient , orixyz ,
402 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
403 ADN_nvals , 1 ,
404 ADN_type , HEAD_ANAT_TYPE ,
405 ADN_view_type , iview ,
406 ADN_func_type , ANAT_MRAN_TYPE ,
407 ADN_none ) ;
408
409 /*-- set byte order (for reading from disk) --*/
410
411 ii = mri_short_order() ;
412 if( swap ) ii = REVERSE_ORDER(ii) ;
413 dset->dblk->diskptr->byte_order = ii ;
414
415 /*-- flag to read data from disk using CTF MRI mode --*/
416
417 dset->dblk->diskptr->storage_mode = STORAGE_BY_CTFMRI ;
418 strcpy( dset->dblk->diskptr->brick_name , fname ) ;
419
420 /*-- for fun, add a set of tags for MEG fiducial points, if present --*/
421
422 if( hh.headModel.LeftEar_Sag != hh.headModel.RightEar_Sag ){
423 THD_usertaglist *tagset = myRwcNew(THD_usertaglist) ;
424 int nas_ii,nas_jj,nas_kk=0 , lft_ii,lft_jj,lft_kk=0 ,
425 rgt_ii,rgt_jj,rgt_kk=0 ;
426 THD_fvec3 fv ; THD_ivec3 iv ;
427
428 tagset->num = 3 ;
429 TAGLIST_SETLABEL( tagset , "CTF MEG Fiducials" ) ;
430
431 /* load voxel indexes into dataset of the 3 tag points;
432 note we have to permute these into the dataset axes order */
433
434 switch( ori[0] ){
435 case 'P':
436 case 'A': nas_ii = hh.headModel.Nasion_Cor ;
437 lft_ii = hh.headModel.LeftEar_Cor ;
438 rgt_ii = hh.headModel.RightEar_Cor ; break ;
439 case 'R':
440 case 'L': nas_ii = hh.headModel.Nasion_Sag ;
441 lft_ii = hh.headModel.LeftEar_Sag ;
442 rgt_ii = hh.headModel.RightEar_Sag ; break ;
443 case 'I':
444 case 'S': nas_ii = hh.headModel.Nasion_Axi ;
445 lft_ii = hh.headModel.LeftEar_Axi ;
446 rgt_ii = hh.headModel.RightEar_Axi ; break ;
447 }
448 switch( ori[1] ){
449 case 'P':
450 case 'A': nas_jj = hh.headModel.Nasion_Cor ;
451 lft_jj = hh.headModel.LeftEar_Cor ;
452 rgt_jj = hh.headModel.RightEar_Cor ; break ;
453 case 'R':
454 case 'L': nas_jj = hh.headModel.Nasion_Sag ;
455 lft_jj = hh.headModel.LeftEar_Sag ;
456 rgt_jj = hh.headModel.RightEar_Sag ; break ;
457 case 'I':
458 case 'S': nas_jj = hh.headModel.Nasion_Axi ;
459 lft_jj = hh.headModel.LeftEar_Axi ;
460 rgt_jj = hh.headModel.RightEar_Axi ; break ;
461 }
462 switch( ori[2] ){
463 case 'P':
464 case 'A': nas_kk = hh.headModel.Nasion_Cor ;
465 lft_kk = hh.headModel.LeftEar_Cor ;
466 rgt_kk = hh.headModel.RightEar_Cor ; break ;
467 case 'R':
468 case 'L': nas_kk = hh.headModel.Nasion_Sag ;
469 lft_kk = hh.headModel.LeftEar_Sag ;
470 rgt_kk = hh.headModel.RightEar_Sag ; break ;
471 case 'I':
472 case 'S': nas_kk = hh.headModel.Nasion_Axi ;
473 lft_kk = hh.headModel.LeftEar_Axi ;
474 rgt_kk = hh.headModel.RightEar_Axi ; break ;
475 }
476
477 TAG_SETLABEL( tagset->tag[0] , "Nasion" ) ;
478 LOAD_IVEC3( iv , nas_ii,nas_jj,nas_kk ) ; /* compute DICOM */
479 fv = THD_3dind_to_3dmm( dset , iv ) ; /* coordinates of */
480 fv = THD_3dmm_to_dicomm( dset , fv ) ; /* this point */
481 UNLOAD_FVEC3( fv , tagset->tag[0].x , tagset->tag[0].y , tagset->tag[0].z ) ;
482 tagset->tag[0].val = 0.0 ;
483 tagset->tag[0].ti = 0 ;
484 tagset->tag[0].set = 1 ;
485
486 TAG_SETLABEL( tagset->tag[1] , "Left Ear" ) ;
487 LOAD_IVEC3( iv , lft_ii,lft_jj,lft_kk ) ;
488 fv = THD_3dind_to_3dmm( dset , iv ) ;
489 fv = THD_3dmm_to_dicomm( dset , fv ) ;
490 UNLOAD_FVEC3( fv , tagset->tag[1].x , tagset->tag[1].y , tagset->tag[1].z ) ;
491 tagset->tag[1].val = 0.0 ;
492 tagset->tag[1].ti = 0 ;
493 tagset->tag[1].set = 1 ;
494
495 TAG_SETLABEL( tagset->tag[2] , "Right Ear" ) ;
496 LOAD_IVEC3( iv , rgt_ii,rgt_jj,rgt_kk ) ;
497 fv = THD_3dind_to_3dmm( dset , iv ) ;
498 fv = THD_3dmm_to_dicomm( dset , fv ) ;
499 UNLOAD_FVEC3( fv , tagset->tag[2].x , tagset->tag[2].y , tagset->tag[2].z ) ;
500 tagset->tag[2].val = 0.0 ;
501 tagset->tag[2].ti = 0 ;
502 tagset->tag[2].set = 1 ;
503
504 dset->tagset = tagset ;
505 }
506
507 RETURN(dset) ;
508 }
509
510 /*------------------------------------------------------------------*/
511 /*! Actually load data from a CTF MRI file into a dataset.
512 Adapted from THD_load_analyze().
513 --------------------------------------------------------------------*/
514
THD_load_ctfmri(THD_datablock * dblk)515 void THD_load_ctfmri( THD_datablock *dblk )
516 {
517 THD_diskptr *dkptr ;
518 int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr=0,nbad ;
519 FILE *fp ;
520 void *ptr ;
521
522 ENTRY("THD_load_ctfmri") ;
523
524 /*-- check inputs --*/
525
526 if( !ISVALID_DATABLOCK(dblk) ||
527 dblk->diskptr->storage_mode != STORAGE_BY_CTFMRI ||
528 dblk->brick == NULL ) EXRETURN ;
529
530 dkptr = dblk->diskptr ;
531
532 /* open and position file at start of data (after header) */
533
534 fp = fopen( dkptr->brick_name , "rb" ) ; /* .mri file */
535 if( fp == NULL ) EXRETURN ;
536
537 /*-- allocate space for data --*/
538
539 nx = dkptr->dimsizes[0] ;
540 ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
541 nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
542 nv = dkptr->nvals ; nxyzv = nxyz * nv ;
543
544 /* 26 Feb 2005: seek backwards from end,
545 instead of forwards from start */
546
547 #if 0
548 fseek( fp , sizeof(Version_2_Header) , SEEK_SET ) ; /* old */
549 #else
550 switch( DBLK_BRICK_TYPE(dblk,0) ){
551 default: ERROR_exit("Unrecognized type in CTF file") ; break ;
552 case MRI_float: ibr = sizeof(float) ; break ; /* illegal */
553 case MRI_short: ibr = sizeof(short) ; break ;
554 case MRI_byte: ibr = sizeof(byte) ; break ;
555 }
556 fseek( fp , -ibr*nxyzv , SEEK_END ) ; /* new */
557 #endif
558
559 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
560
561 /*-- malloc space for each brick separately --*/
562
563 for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
564 if( DBLK_ARRAY(dblk,ibr) == NULL ){
565 ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
566 mri_fix_data_pointer( ptr , DBLK_BRICK(dblk,ibr) ) ;
567 if( ptr == NULL ) nbad++ ;
568 }
569 }
570
571 /*-- if couldn't get them all, take our ball and go home in a snit --*/
572
573 if( nbad > 0 ){
574 fprintf(stderr,
575 "\n** failed to malloc %d CTR MRI bricks out of %d\n\a",nbad,nv);
576 for( ibr=0 ; ibr < nv ; ibr++ ){
577 if( DBLK_ARRAY(dblk,ibr) != NULL ){
578 free(DBLK_ARRAY(dblk,ibr)) ;
579 mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
580 }
581 }
582 fclose(fp) ; EXRETURN ;
583 }
584
585
586 /*-- read data from .img file into sub-brick arrays! --*/
587
588 for( ibr=0 ; ibr < nv ; ibr++ )
589 fread( DBLK_ARRAY(dblk,ibr), 1, DBLK_BRICK_BYTES(dblk,ibr), fp ) ;
590
591 fclose(fp) ;
592
593 /*-- swap bytes? --*/
594
595 if( dkptr->byte_order != mri_short_order() ){
596 for( ibr=0 ; ibr < nv ; ibr++ ){
597 switch( DBLK_BRICK_TYPE(dblk,ibr) ){
598 default: break ;
599 case MRI_short:
600 mri_swap2( DBLK_BRICK_NVOX(dblk,ibr) , DBLK_ARRAY(dblk,ibr) ) ;
601 break ;
602 }
603 }
604 }
605
606 EXRETURN ;
607 }
608
609 /****************************************************************************
610 SAM static image files are structured as follows:
611
612 char Identity[8] = "SAMIMAGE"; // uniquely identifies image file
613 SAM_HDR SAMHeader; // SAM header
614 double Voxel[0]; // 1st SAM voxel (units=A-m, (A-m)^2,
615 double Voxel[1]; // 2nd SAM voxel Z, T, F, or P)
616 "
617 "
618 double Voxel[V]; // last SAM voxel
619
620 Coefficients & image voxels are ordered in X,Y,Z sequence, with Z the least
621 significant index (most rapidly changing), Y is next, and then X.
622 Coordinate indices always advance in the positive direction. This implies
623 that Voxel[0] is in the right, posterior, inferior position relative to
624 the region of interest (bounding box of image).
625
626 RWCox: the data storage order seems to be IRP, based on the above
627 comments, and on the CTF head coordinates system being PRI.
628 *****************************************************************************/
629
630 #define COV_VERSION 1 /* this is version 1 -- got it? */
631 #define SAM_VERSION 1 /* this, too! */
632
633 /** SAM file types **/
634
635 #define SAM_TYPE_IMAGE 0 /* flags file as a SAM static image file */
636 #define SAM_TYPE_WT_ARRAY 1 /* flags file as SAM coefficients for regular target array */
637 #define SAM_TYPE_WT_LIST 2 /* flags file as SAM coefficients for target list */
638
639 /** define SAM unit types **/
640
641 #define SAM_UNIT_COEFF 0 /* SAM coefficients A-m/T */
642 #define SAM_UNIT_MOMENT 1 /* SAM source (or noise) strength A-m */
643 #define SAM_UNIT_POWER 2 /* SAM source (or noise) power (A-m)^2 */
644 #define SAM_UNIT_SPMZ 3 /* SAM z-deviate */
645 #define SAM_UNIT_SPMF 4 /* SAM F-statistic */
646 #define SAM_UNIT_SPMT 5 /* SAM T-statistic */
647 #define SAM_UNIT_SPMP 6 /* SAM probability */
648 #define SAM_UNIT_MUSIC 7 /* MUSIC metric */
649
650 /* 'SAM_HDR' is to be used for both SAM coefficients & SAM static images */
651 typedef struct {
652 int Version; /* file version number (should be 1) */
653 char SetName[256]; /* name of parent dataset */
654 int NumChans; /* # of channels used by SAM */
655 int NumWeights; /* # of SAM virtual channels (0=static image) */
656 int pad_bytes1; /* ** align next double on 8 byte boundary */
657 double XStart; /* x-start coordinate (m) */
658 double XEnd; /* x-end coordinate (m) */
659 double YStart; /* y-start coordinate (m) */
660 double YEnd; /* y-end coordinate (m) */
661 double ZStart; /* z-start coordinate (m) */
662 double ZEnd; /* z-end coordinate (m) */
663 double StepSize; /* voxel step size (m) */
664 double HPFreq; /* highpass frequency (Hz) */
665 double LPFreq; /* lowpass frequency (Hz) */
666 double BWFreq; /* bandwidth of filters (Hz) */
667 double MeanNoise; /* mean primary sensor noise (T) */
668 char MriName[256]; /* MRI image file name */
669 int Nasion[3]; /* MRI voxel index for nasion */
670 int RightPA[3]; /* MRI voxel index for right pre-auricular */
671 int LeftPA[3]; /* MRI voxel index for left pre-auricular */
672 int SAMType; /* SAM file type */
673 int SAMUnit; /* SAM units (a bit redundant, but may be useful) */
674 int pad_bytes2; /* ** align end of structure on 8 byte boundary */
675 } SAM_HDR;
676
677 /*** 26 Feb 2005: version 2 of the SAM header ***/
678
679 #if 0
680 typedef struct {
681 int Version; /* file version number (should be 2) */
682 char SetName[256]; /* name of parent dataset */
683 int NumChans; /* # of channels used by SAM */
684 int NumWeights; /* # of SAM virtual channels (0=static image) */
685 int pad_bytes1; /* ** align next double on 8 byte boundary */
686 double XStart; /* x-start coordinate (m) */
687 double XEnd; /* x-end coordinate (m) */
688 double YStart; /* y-start coordinate (m) */
689 double YEnd; /* y-end coordinate (m) */
690 double ZStart; /* z-start coordinate (m) */
691 double ZEnd; /* z-end coordinate (m) */
692 double StepSize; /* voxel step size (m) */
693 double HPFreq; /* highpass frequency (Hz) */
694 double LPFreq; /* lowpass frequency (Hz) */
695 double BWFreq; /* bandwidth of filters (Hz) */
696 double MeanNoise; /* mean primary sensor noise (T) */
697 char MriName[256]; /* MRI image file name */
698 int Nasion[3]; /* MRI voxel index for nasion */
699 int RightPA[3]; /* MRI voxel index for right pre-auricular */
700 int LeftPA[3]; /* MRI voxel index for left pre-auricular */
701 int SAMType; /* SAM file type */
702 int SAMUnit; /* SAM units (a bit redundant, but may be useful) */
703 int pad_bytes2; /* ** align end of structure on 8 byte boundary */
704 double MegNasion[3]; /* MEG dewar coordinates for nasion (m) */
705 double MegRightPA[3]; /* MEG dewar coordinates for R pre-auricular (m) */
706 double MegLeftPA[3]; /* MEG dewar coordinates for L pre-auricular (m) */
707 char SAMUnitName[32]; /* SAM units (redundant, but useful too!) */
708 } SAM_HDR_v2;
709 #endif
710
711
712 /*-------------------------*/
713 /*! Macro for bad return. */
714
715 #undef BADBAD
716 #define BADBAD(s) \
717 do{ fprintf(stderr,"** THD_open_ctfsam(%s): %s\n",fname,s); \
718 RETURN(NULL); \
719 } while(0)
720
721 /*-----------------------------------------------------------------*/
722 /*! Open a CTF .svl (SAM) file as an unpopulated AFNI dataset.
723 It will be populated later, in THD_load_ctfsam().
724 -------------------------------------------------------------------*/
725
THD_open_ctfsam(char * fname)726 THD_3dim_dataset * THD_open_ctfsam( char *fname )
727 {
728 FILE *fp ;
729 SAM_HDR hh ;
730 char Identity[9] ;
731 int ii,nn , swap ;
732 THD_3dim_dataset *dset=NULL ;
733 char prefix[THD_MAX_PREFIX] , *ppp , tname[12] , ori[4] ;
734 THD_ivec3 nxyz , orixyz ;
735 THD_fvec3 dxyz , orgxyz ;
736 int iview ;
737 int ngood , length , datum_type , datum_len , oxx,oyy,ozz ;
738 int nx,ny,nz ;
739 float dx,dy,dz , xorg,yorg,zorg ;
740
741 ENTRY("THD_open_ctfsam") ;
742
743 /* open input file */
744
745 if( fname == NULL || *fname == '\0' ) BADBAD("bad input filename");
746 fp = fopen( fname , "rb" ) ;
747 if( fp == NULL ) BADBAD("can't open input file");
748
749 /* read header [1st 8 bytes are "SAMIMAGE"] */
750
751 fread( Identity , 1,8 , fp ) ; Identity[8] = '\0' ;
752 fread( &hh , sizeof(hh) , 1 , fp ) ;
753 fclose(fp) ;
754
755 if( strcmp(Identity,"SAMIMAGE") != 0 ) BADBAD("Identity != SAMIMAGE") ;
756 if( hh.Version == 0 ) BADBAD("bad header Version") ;
757
758 swap = (hh.Version < 0) || (hh.Version > 3) ; /* byte swap? */
759
760 if( swap ){ /* swap various header fields */
761 swap_4( &hh.Version ) ;
762 swap_4( &hh.NumChans ) ;
763 swap_4( &hh.NumWeights ) ;
764 swap_8( &hh.XStart ) ;
765 swap_8( &hh.XEnd ) ;
766 swap_8( &hh.YStart ) ;
767 swap_8( &hh.YEnd ) ;
768 swap_8( &hh.ZStart ) ;
769 swap_8( &hh.ZEnd ) ;
770 swap_8( &hh.StepSize ) ;
771 swap_8( &hh.HPFreq ) ;
772 swap_8( &hh.LPFreq ) ;
773 swap_8( &hh.BWFreq ) ;
774 swap_8( &hh.MeanNoise ) ;
775 swap_4( &hh.Nasion[0] ) ;
776 swap_4( &hh.RightPA[0] ) ;
777 swap_4( &hh.LeftPA[0] ) ;
778 swap_4( &hh.Nasion[1] ) ;
779 swap_4( &hh.RightPA[1] ) ;
780 swap_4( &hh.LeftPA[1] ) ;
781 swap_4( &hh.Nasion[2] ) ;
782 swap_4( &hh.RightPA[2] ) ;
783 swap_4( &hh.LeftPA[2] ) ;
784 swap_4( &hh.SAMType ) ;
785 swap_4( &hh.SAMUnit ) ;
786 }
787
788 /* simple checks on header values */
789
790 if( hh.Version < 0 ||
791 hh.Version > 3 || /* 26 Feb 2005 */
792 hh.XStart >= hh.XEnd ||
793 hh.YStart >= hh.YEnd ||
794 hh.ZStart >= hh.ZEnd ||
795 hh.StepSize <= 0.0 ) BADBAD("bad header data") ;
796
797 #if 0
798 printf("\n") ;
799 printf("**CTF SAM : %s\n",fname) ;
800 printf("Version = %d\n",hh.Version) ;
801 printf("NumChans = %d\n",hh.NumChans) ;
802 printf("NumWeights= %d\n",hh.NumWeights) ;
803 printf("XStart = %g\n",hh.XStart) ;
804 printf("Xend = %g\n",hh.XEnd) ;
805 printf("YStart = %g\n",hh.YStart) ;
806 printf("YEnd = %g\n",hh.YEnd) ;
807 printf("ZStart = %g\n",hh.ZStart) ;
808 printf("Zend = %g\n",hh.ZEnd) ;
809 printf("StepSize = %g\n",hh.StepSize) ;
810 printf("HPFreq = %g\n",hh.HPFreq) ;
811 printf("LPFreq = %g\n",hh.LPFreq) ;
812 printf("BWFreq = %g\n",hh.BWFreq) ;
813 printf("MeanNoise = %g\n",hh.MeanNoise) ;
814 printf("Nasion[0] = %d\n",hh.Nasion[0]) ;
815 printf("Nasion[1] = %d\n",hh.Nasion[1]) ;
816 printf("Nasion[2] = %d\n",hh.Nasion[2]) ;
817 printf("RightPA[0]= %d\n",hh.RightPA[0]) ;
818 printf("RightPA[1]= %d\n",hh.RightPA[1]) ;
819 printf("RightPA[2]= %d\n",hh.RightPA[2]) ;
820 printf("LeftPA[0] = %d\n",hh.LeftPA[0]) ;
821 printf("LeftPA[1] = %d\n",hh.LeftPA[1]) ;
822 printf("LeftPA[2] = %d\n",hh.LeftPA[2]) ;
823 printf("SAMtype = %d\n",hh.SAMType) ;
824 printf("SAMunit = %d\n",hh.SAMUnit) ;
825 printf("SetName = %s\n",hh.SetName) ;
826 printf("MriName = %s\n",hh.MriName) ;
827 printf("headersize= %d\n",sizeof(hh)+8) ;
828 #endif
829
830 hh.StepSize *= 1000.0 ; /* convert distances from m to mm */
831 hh.XStart *= 1000.0 ; /* (who the hell uses meters for brain imaging?) */
832 hh.YStart *= 1000.0 ; /* (blue whales? elephants?) */
833 hh.ZStart *= 1000.0 ;
834 hh.XEnd *= 1000.0 ;
835 hh.YEnd *= 1000.0 ;
836 hh.ZEnd *= 1000.0 ;
837
838 dx = dy = dz = hh.StepSize ; /* will be altered below */
839
840 #if 0
841 nx = (int)((hh.ZEnd - hh.ZStart)/dz + 0.99999); /* dataset is stored in Z,Y,X order */
842 ny = (int)((hh.YEnd - hh.YStart)/dy + 0.99999); /* but AFNI calls these x,y,z */
843 nz = (int)((hh.XEnd - hh.XStart)/dx + 0.99999);
844 #else
845 nx = CTF_count( hh.ZStart , hh.ZEnd , hh.StepSize ) ;
846 ny = CTF_count( hh.YStart , hh.YEnd , hh.StepSize ) ;
847 nz = CTF_count( hh.XStart , hh.XEnd , hh.StepSize ) ;
848 #endif
849
850 /* determine if file is big enough to hold all data it claims */
851
852 nn = THD_filesize(fname) ;
853 if( nn < sizeof(double)*nx*ny*nz ) BADBAD("input file too small") ;
854
855 datum_type = MRI_float ; /* actually is double, but AFNI doesn't grok that */
856 /* will be converted to floats when reading data */
857
858 /* set orientation = IRP = xyz ordering */
859
860 ori[0] = 'I'; ori[1] = 'R'; ori[2] = 'P';
861
862 oxx = ORCODE(ori[0]); oyy = ORCODE(ori[1]); ozz = ORCODE(ori[2]);
863 if( !OR3OK(oxx,oyy,ozz) ){
864 oxx = ORI_I2S_TYPE; oyy = ORI_R2L_TYPE; ozz = ORI_P2A_TYPE; /** IRP? **/
865 }
866
867 orixyz.ijk[0] = oxx ; orixyz.ijk[1] = oyy ; orixyz.ijk[2] = ozz ;
868
869 /* now set grid size, keeping in mind that
870 A-P is positive and P-A is negative,
871 R-L is positive and L-R is negative,
872 I-S is positive and S-I is negative. */
873
874 switch( ori[0] ){
875 case 'A': dx = hh.StepSize ; xorg = -hh.XStart ; break ;
876 case 'P': dx = -hh.StepSize ; xorg = -hh.XStart ; break ;
877 case 'R': dx = hh.StepSize ; xorg = hh.YStart ; break ;
878 case 'L': dx = -hh.StepSize ; xorg = hh.YStart ; break ;
879 case 'I': dx = hh.StepSize ; xorg = hh.ZStart ; break ;
880 case 'S': dx = -hh.StepSize ; xorg = hh.ZStart ; break ;
881 }
882 switch( ori[1] ){
883 case 'A': dy = hh.StepSize ; yorg = -hh.XStart ; break ;
884 case 'P': dy = -hh.StepSize ; yorg = -hh.XStart ; break ;
885 case 'R': dy = hh.StepSize ; yorg = hh.YStart ; break ;
886 case 'L': dy = -hh.StepSize ; yorg = hh.YStart ; break ;
887 case 'I': dy = hh.StepSize ; yorg = hh.ZStart ; break ;
888 case 'S': dy = -hh.StepSize ; yorg = hh.ZStart ; break ;
889 }
890 switch( ori[2] ){
891 case 'A': dz = hh.StepSize ; zorg = -hh.XStart ; break ;
892 case 'P': dz = -hh.StepSize ; zorg = -hh.XStart ; break ;
893 case 'R': dz = hh.StepSize ; zorg = hh.YStart ; break ;
894 case 'L': dz = -hh.StepSize ; zorg = hh.YStart ; break ;
895 case 'I': dz = hh.StepSize ; zorg = hh.ZStart ; break ;
896 case 'S': dz = -hh.StepSize ; zorg = hh.ZStart ; break ;
897 }
898
899 /*-- make a dataset --*/
900
901 dset = EDIT_empty_copy(NULL) ;
902
903 dset->idcode.str[0] = 'C' ; /* overwrite 1st 3 bytes */
904 dset->idcode.str[1] = 'T' ;
905 dset->idcode.str[2] = 'F' ;
906
907 MCW_hash_idcode( fname , dset ) ; /* 06 May 2005 */
908
909 ppp = THD_trailname(fname,0) ; /* strip directory */
910 MCW_strncpy( prefix , ppp , THD_MAX_PREFIX ) ; /* to make prefix */
911
912 nxyz.ijk[0] = nx ; dxyz.xyz[0] = dx ; /* setup axes lengths and voxel sizes */
913 nxyz.ijk[1] = ny ; dxyz.xyz[1] = dy ;
914 nxyz.ijk[2] = nz ; dxyz.xyz[2] = dz ;
915
916 orixyz.ijk[0] = oxx ; orgxyz.xyz[0] = xorg ;
917 orixyz.ijk[1] = oyy ; orgxyz.xyz[1] = yorg ;
918 orixyz.ijk[2] = ozz ; orgxyz.xyz[2] = zorg ;
919
920 iview = VIEW_ORIGINAL_TYPE ;
921
922 /*-- actually send the values above into the dataset header --*/
923
924 EDIT_dset_items( dset ,
925 ADN_prefix , prefix ,
926 ADN_datum_all , datum_type ,
927 ADN_nxyz , nxyz ,
928 ADN_xyzdel , dxyz ,
929 ADN_xyzorg , orgxyz ,
930 ADN_xyzorient , orixyz ,
931 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
932 ADN_nvals , 1 ,
933 ADN_type , HEAD_FUNC_TYPE ,
934 ADN_view_type , iview ,
935 ADN_func_type , FUNC_FIM_TYPE ,
936 ADN_none ) ;
937
938 /*-- set byte order (for reading from disk) --*/
939
940 ii = mri_short_order() ;
941 if( swap ) ii = REVERSE_ORDER(ii) ;
942 dset->dblk->diskptr->byte_order = ii ;
943
944 /*-- flag to read data from disk using CTF SAM mode --*/
945
946 dset->dblk->diskptr->storage_mode = STORAGE_BY_CTFSAM ;
947 strcpy( dset->dblk->diskptr->brick_name , fname ) ;
948
949 RETURN(dset) ;
950 }
951
952 /*------------------------------------------------------------------*/
953 /*! Actually load data from a CTF SAM file into a dataset.
954 --------------------------------------------------------------------*/
955
THD_load_ctfsam(THD_datablock * dblk)956 void THD_load_ctfsam( THD_datablock *dblk )
957 {
958 THD_diskptr *dkptr ;
959 int nx,ny,nz,nv , nxy,nxyz,nxyzv , ibr,nbad , ii,swap ;
960 FILE *fp ;
961 void *ptr ;
962 double *dbar ;
963 float *ftar ;
964
965 ENTRY("THD_load_ctfsam") ;
966
967 /*-- check inputs --*/
968
969 if( !ISVALID_DATABLOCK(dblk) ||
970 dblk->diskptr->storage_mode != STORAGE_BY_CTFSAM ||
971 dblk->brick == NULL ) EXRETURN ;
972
973 dkptr = dblk->diskptr ;
974
975 /*-- allocate space for data --*/
976
977 nx = dkptr->dimsizes[0] ;
978 ny = dkptr->dimsizes[1] ; nxy = nx * ny ;
979 nz = dkptr->dimsizes[2] ; nxyz = nxy * nz ;
980 nv = dkptr->nvals ; nxyzv = nxyz * nv ;
981
982 /* position file 8*nxyzv bytes before end of file */
983
984 fp = fopen( dkptr->brick_name , "rb" ) ; /* .svl file */
985 if( fp == NULL ) EXRETURN ;
986
987 /* 26 Feb 2005: instead of skipping the header,
988 whose size varies with the SAM version number,
989 just seek backwards from the end to the correct size */
990
991 fseek( fp , -sizeof(double)*nxyzv , SEEK_END ) ;
992
993 dblk->malloc_type = DATABLOCK_MEM_MALLOC ;
994
995 /*-- malloc space for each brick separately --*/
996
997 for( nbad=ibr=0 ; ibr < nv ; ibr++ ){
998 if( DBLK_ARRAY(dblk,ibr) == NULL ){
999 ptr = AFMALL(void, DBLK_BRICK_BYTES(dblk,ibr) ) ;
1000 mri_fix_data_pointer( ptr , DBLK_BRICK(dblk,ibr) ) ;
1001 if( ptr == NULL ) nbad++ ;
1002 }
1003 }
1004
1005 /*-- if couldn't get them all, take our ball and go home in a snit --*/
1006
1007 if( nbad > 0 ){
1008 fprintf(stderr,
1009 "\n** failed to malloc %d CTR MRI bricks out of %d\n\a",nbad,nv);
1010 for( ibr=0 ; ibr < nv ; ibr++ ){
1011 if( DBLK_ARRAY(dblk,ibr) != NULL ){
1012 free(DBLK_ARRAY(dblk,ibr)) ;
1013 mri_fix_data_pointer( NULL , DBLK_BRICK(dblk,ibr) ) ;
1014 }
1015 }
1016 fclose(fp) ; EXRETURN ;
1017 }
1018
1019 /*-- SAM data is stored as doubles,
1020 but we have to store it in AFNI as floats --*/
1021
1022 dbar = (double *) calloc(sizeof(double),nxyz) ; /* workspace */
1023 swap = ( dkptr->byte_order != mri_short_order() ) ;
1024
1025 for( ibr=0 ; ibr < nv ; ibr++ ){ /* loop over sub-bricks */
1026 fread( dbar, 1, sizeof(double)*nxyz, fp ) ; /* read data to workspace */
1027 ftar = DBLK_ARRAY(dblk,ibr) ; /* float array in dataset */
1028 for( ii=0 ; ii < nxyz ; ii++ ){ /* loop over voxels */
1029 if( swap ) swap_8(dbar+ii) ; /* swap it */
1030 ftar[ii] = dbar[ii] ; /* save it as a float */
1031 }
1032 }
1033
1034 fclose(fp) ; free(dbar) ; /* toss out the trash */
1035 EXRETURN ;
1036 }
1037