1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 * filename: m-anlz.c *
3 * *
4 * UTIL C-source: Medical Image Conversion Utility *
5 * *
6 * purpose : read and write ANALYZE files *
7 * *
8 * project : (X)MedCon by Erik Nolf *
9 * *
10 * Functions : MdcCheckANLZ() - Check ANALYZE format *
11 * MdcReadANLZ() - Read ANALYZE file *
12 * MdcWriteANLZ() - Write ANALYZE file *
13 * MdcWriteHeaderKey() - Write Header Key to file *
14 * MdcWriteImageDimension() - Write Image Dimension to file *
15 * MdcWriteDataHistory() - Write Data History to file *
16 * MdcWriteImagesData() - Write the images to file *
17 * MdcGetSpmOpt() - Get specific SPM options *
18 * *
19 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
20 /*
21 */
22
23 /*
24 Copyright (C) 1997-2021 by Erik Nolf
25
26 This program is free software; you can redistribute it and/or modify it
27 under the terms of the GNU General Public License as published by the
28 Free Software Foundation; either version 2, or (at your option) any later
29 version.
30
31 This program is distributed in the hope that it will be useful, but
32 WITHOUT ANY WARRANTY; without even the implied warranty of
33 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
34 Public License for more details.
35
36 You should have received a copy of the GNU General Public License along
37 with this program; if not, write to the Free Software Foundation, Inc.,
38 59 Place - Suite 330, Boston, MA 02111-1307, USA. */
39
40 /****************************************************************************
41 H E A D E R S
42 ****************************************************************************/
43
44 #include "m-depend.h"
45
46 #include <stdio.h>
47 #ifdef HAVE_STDLIB_H
48 #include <stdlib.h>
49 #endif
50 #ifdef HAVE_STRING_H
51 #include <string.h>
52 #endif
53 #ifdef HAVE_STRINGS_H
54 #ifndef _WIN32
55 #include <strings.h>
56 #endif
57 #endif
58 #ifdef HAVE_UNISTD_H
59 #include <unistd.h>
60 #endif
61
62 #include "medcon.h"
63
64 /****************************************************************************
65 D E F I N E S
66 *****************************************************************************/
67
68 static Int8 INIT_SPMOPT = MDC_YES;
69 static MDC_SPMOPT spmopt;
70
71 #define MDC_ALWAYS_SET_4D 1 /* 0/1 disable/enable always set 4 dims */
72
73 /****************************************************************************
74 F U N C T I O N S
75 *****************************************************************************/
76
MdcCheckANLZ(FILEINFO * fi)77 int MdcCheckANLZ(FILEINFO *fi)
78 {
79 MDC_ANLZ_HEADER_KEY hk;
80 int check=2, FORMAT=MDC_FRMT_NONE;
81
82 if (fread((char *)&hk,1,MDC_ANLZ_HK_SIZE,fi->ifp) != MDC_ANLZ_HK_SIZE)
83 return(MDC_BAD_READ);
84
85 MDC_FILE_ENDIAN = MDC_HOST_ENDIAN;
86
87 while (check--) {
88
89 if ( (hk.sizeof_hdr==348 || hk.sizeof_hdr==148 || hk.sizeof_hdr==228 ||
90 hk.sizeof_hdr==384)
91 && (hk.regular == MDC_ANLZ_SIG ) ) {
92 FORMAT = MDC_FRMT_ANLZ;
93 break;
94 }
95 MDC_FILE_ENDIAN = !MDC_HOST_ENDIAN;
96
97 MdcSWAP(hk.sizeof_hdr);
98
99 }
100
101 return(FORMAT);
102
103 }
104
MdcReadANLZ(FILEINFO * fi)105 const char *MdcReadANLZ(FILEINFO *fi)
106 {
107 MDC_SPMOPT *opt = &spmopt;
108 FILE *fp=fi->ifp;
109 MDC_ANLZ_HEADER_KEY hk;
110 MDC_ANLZ_IMAGE_DIMS imd;
111 MDC_ANLZ_DATA_HIST dh;
112 IMG_DATA *id=NULL;
113 DYNAMIC_DATA *dd=NULL;
114 Uint32 bytes, i, plane, f, number;
115 Uint8 *img8=NULL;
116 Int8 WAS_COMPRESSED = MDC_NO;
117 char *origpath=NULL;
118 const char *err=NULL;
119
120 if (MDC_FILE_STDIN == MDC_YES)
121 return("ANLZ File input from stdin unsupported");
122
123 if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Reading Analyze:");
124
125 if (MDC_VERBOSE) MdcPrntMesg("ANLZ Reading <%s> ...",fi->ifname);
126
127 /* get endian of the file in MDC_FILE_ENDIAN */
128 i=MdcCheckANLZ(fi); fseek(fp,0,SEEK_SET);
129 if (i != MDC_FRMT_ANLZ) {
130 if (MDC_FALLBACK_FRMT == MDC_FRMT_ANLZ) {
131 /* set host endian to try analyze */
132 MDC_FILE_ENDIAN = MDC_HOST_ENDIAN;
133 }else{
134 /* bail out */
135 return("ANLZ Endian check failed");
136 }
137 }
138
139 memset(&hk,0,MDC_ANLZ_HK_SIZE);
140 memset(&imd,0,MDC_ANLZ_IMD_SIZE);
141 memset(&dh,0,MDC_ANLZ_DH_SIZE);
142
143 /* put some default we use */
144 fi->reconstructed = MDC_YES;
145 fi->acquisition_type = MDC_ACQUISITION_TOMO;
146
147 dh.orient=(char)0xff;
148
149 if (fread((char *)&hk,1,MDC_ANLZ_HK_SIZE,fp) != MDC_ANLZ_HK_SIZE)
150 return("ANLZ Bad read HeadKey struct");
151
152
153 fi->endian=MDC_FILE_ENDIAN;
154
155 MdcSWAP(hk.sizeof_hdr); MdcSWAP(hk.extents); MdcSWAP(hk.session_error);
156
157 if (MDC_INFO) {
158 MdcPrntScrn("\nMDC_ANLZ_HEADER_KEY (%d bytes)\n",MDC_ANLZ_HK_SIZE);
159 MdcPrintLine('-',MDC_HALF_LENGTH);
160 MdcPrntScrn("sizeof_hdr : %d\n",hk.sizeof_hdr);
161 strncpy(mdcbufr,hk.data_type,10); mdcbufr[10]='\0';
162 MdcPrntScrn("data_type : "); MdcPrintStr(mdcbufr);
163 strncpy(mdcbufr,hk.db_name,18); mdcbufr[18]='\0';
164 MdcPrntScrn("db_name : "); MdcPrintStr(mdcbufr);
165 MdcPrntScrn("extents : %d\n",hk.extents);
166 MdcPrntScrn("session_error : %hd\n",hk.session_error);
167 MdcPrntScrn("regular : "); MdcPrintChar(hk.regular);
168 MdcPrntScrn("\n");
169 MdcPrntScrn("hkey_un0 : "); MdcPrintChar(hk.hkey_un0);
170 MdcPrntScrn("\n");
171 }
172
173 if (MDC_INFO) {
174 MdcPrntScrn("\nIMAGE_DIMENSION (%d bytes)\n",MDC_ANLZ_IMD_SIZE);
175 MdcPrintLine('-',MDC_HALF_LENGTH);
176 }
177 if (fread((char *)&imd,1,MDC_ANLZ_IMD_SIZE,fp)!=MDC_ANLZ_IMD_SIZE)
178 return("ANLZ Bad read ImageDimensions struct");
179
180 for (i=0; i<MDC_ANLZ_MAX_DIMS; i++) {
181 MdcSWAP(imd.dim[i]); MdcSWAP(imd.pixdim[i]);
182 }
183 MdcSWAP(imd.unused1);
184 MdcSWAP(imd.avw_vox_offset);
185 MdcSWAP(imd.spm_pix_rescale);
186 MdcSWAP(imd.funused1); MdcSWAP(imd.funused2);
187 MdcSWAP(imd.avw_cal_max); MdcSWAP(imd.avw_cal_min);
188 MdcSWAP(imd.datatype); MdcSWAP(imd.bitpix); MdcSWAP(imd.dim_un0);
189 MdcSWAP(imd.compressed); MdcSWAP(imd.verified);
190 MdcSWAP(imd.glmax); MdcSWAP(imd.glmin);
191
192 if ((imd.dim[0] < 3) || (imd.dim[0] > (MDC_ANLZ_MAX_DIMS - 1))) {
193 if (MDC_FALLBACK_FRMT == MDC_FRMT_ANLZ) {
194 /* force reading, set to an acceptable value */
195 /* hope unused dims were initialized to zero */
196 for (i = 3; i < MDC_ANLZ_MAX_DIMS; i++) if (imd.dim[i] <= 0) break;
197 imd.dim[0] = i-1;
198 MdcPrntWarn("ANLZ Bad header value in dim[0] dimension");
199 }else{
200 /* bail out safely */
201 return("ANLZ Bad header value in dim[0] dimension");
202 }
203 }
204
205 if (MDC_INFO) {
206 for (i=0; i<MDC_ANLZ_MAX_DIMS; i++)
207 MdcPrntScrn("dim[%d] : %hd\n",i,imd.dim[i]);
208 strncpy(mdcbufr,imd.avw_vox_units,4); mdcbufr[4]='\0';
209 MdcPrntScrn("avw_vox_units : "); MdcPrintStr(mdcbufr);
210 strncpy(mdcbufr,imd.avw_cal_units,8); mdcbufr[8]='\0';
211 MdcPrntScrn("avw_cal_units : "); MdcPrintStr(mdcbufr);
212 MdcPrntScrn("unused1 : %hd\n",imd.unused1);
213 MdcPrntScrn("datatype : ");
214 switch (imd.datatype) {
215 case MDC_ANLZ_DT_UNKNOWN : MdcPrntScrn("Unknown");
216 break;
217 case MDC_ANLZ_DT_BINARY : MdcPrntScrn("Binary");
218 break;
219 case MDC_ANLZ_DT_UNSIGNED_CHAR: MdcPrntScrn("Unsigned character");
220 break;
221 case MDC_ANLZ_DT_SIGNED_SHORT : MdcPrntScrn("Signed short integer");
222 break;
223 case MDC_ANLZ_DT_SIGNED_INT : MdcPrntScrn("Signed integer");
224 break;
225 case MDC_ANLZ_DT_FLOAT : MdcPrntScrn("Floating point");
226 break;
227 case MDC_ANLZ_DT_COMPLEX : MdcPrntScrn("Complex");
228 break;
229 case MDC_ANLZ_DT_DOUBLE : MdcPrntScrn("Double precision");
230 break;
231 case MDC_ANLZ_DT_RGB : MdcPrntScrn("Coloured RGB");
232 break;
233 case MDC_ANLZ_DT_ALL : MdcPrntScrn("All");
234 break;
235 default: MdcPrntScrn("Unknown");
236 }
237 MdcPrntScrn("\n");
238 MdcPrntScrn("bitpix : %hd\n",imd.bitpix);
239 for (i=0; i<MDC_ANLZ_MAX_DIMS; i++) {
240 MdcPrntScrn("pixdim[%d] : %-6e ",i,imd.pixdim[i]);
241 if ( i<4 ) MdcPrntScrn("[mm]\n");
242 else MdcPrntScrn("[ms]\n");
243 }
244 MdcPrntScrn("avw_vox_offset : %-6e\n",imd.avw_vox_offset);
245 MdcPrntScrn("spm_pix_rescale: %-6e\n",imd.spm_pix_rescale);
246 MdcPrntScrn("funused1 : %-6e\n",imd.funused1);
247 MdcPrntScrn("funused2 : %-6e\n",imd.funused2);
248 MdcPrntScrn("avw_cal_max : %-6e\n",imd.avw_cal_max);
249 MdcPrntScrn("avw_cal_min : %-6e\n",imd.avw_cal_min);
250 MdcPrntScrn("compressed : %6f\n",imd.compressed);
251 MdcPrntScrn("verified : %6f\n",imd.verified);
252 MdcPrntScrn("glmax : %d\n",imd.glmax);
253 MdcPrntScrn("glmin : %d\n",imd.glmin);
254
255 }
256
257 if (MDC_INFO) {
258 MdcPrntScrn("\nMDC_ANLZ_DATA_HISTORY (%d bytes)\n", MDC_ANLZ_DH_SIZE);
259 MdcPrintLine('-',MDC_HALF_LENGTH);
260 }
261 if (fread((char *)&dh,1,MDC_ANLZ_DH_SIZE,fp) != MDC_ANLZ_DH_SIZE) {
262 MdcPrntWarn("Reading <%s> - ANLZ Truncated header",fi->ifname);
263 }
264
265 memcpy(&opt->origin_x,&dh.originator[0],2); MdcSWAP(opt->origin_x);
266 memcpy(&opt->origin_y,&dh.originator[2],2); MdcSWAP(opt->origin_y);
267 memcpy(&opt->origin_z,&dh.originator[4],2); MdcSWAP(opt->origin_z);
268
269
270 MdcSWAP(dh.views); MdcSWAP(dh.vols_added);
271 MdcSWAP(dh.start_field); MdcSWAP(dh.field_skip);
272 MdcSWAP(dh.omax); MdcSWAP(dh.omin);
273 MdcSWAP(dh.smax); MdcSWAP(dh.smin);
274
275 if (MDC_INFO) {
276 strncpy(mdcbufr,dh.descrip,80); mdcbufr[80]='\0';
277 MdcPrntScrn("description : "); MdcPrintStr(mdcbufr);
278 strncpy(mdcbufr,dh.aux_file,24); mdcbufr[24]='\0';
279 MdcPrntScrn("aux_file : "); MdcPrintStr(mdcbufr);
280 MdcPrntScrn("orient : ");
281 switch (dh.orient) {
282 case MDC_ANLZ_TRANS_UNFLIPPED: MdcPrntScrn("transverse unflipped");
283 break;
284 case MDC_ANLZ_CORON_UNFLIPPED: MdcPrntScrn("coronal unflipped");
285 break;
286 case MDC_ANLZ_SAGIT_UNFLIPPED: MdcPrntScrn("sagittal unflipped");
287 break;
288 case MDC_ANLZ_TRANS_FLIPPED : MdcPrntScrn("transverse flipped");
289 break;
290 case MDC_ANLZ_CORON_FLIPPED : MdcPrntScrn("coronal flipped");
291 break;
292 case MDC_ANLZ_SAGIT_FLIPPED : MdcPrntScrn("sagittal flipped");
293 break;
294 default: MdcPrntScrn("Unknown");
295 }
296 MdcPrntScrn("\n");
297 strncpy(mdcbufr,dh.originator,10); mdcbufr[10]='\0';
298 MdcPrntScrn("originator : "); MdcPrintStr(mdcbufr);
299 strncpy(mdcbufr,dh.generated,10); mdcbufr[10]='\0';
300 MdcPrntScrn("generated : "); MdcPrintStr(mdcbufr);
301 strncpy(mdcbufr,dh.scannum,10); mdcbufr[10]='\0';
302 MdcPrntScrn("scannum : "); MdcPrintStr(mdcbufr);
303 strncpy(mdcbufr,dh.patient_id,10); mdcbufr[10]='\0';
304 MdcPrntScrn("patient_id : "); MdcPrintStr(mdcbufr);
305 strncpy(mdcbufr,dh.exp_date,10); mdcbufr[10]='\0';
306 MdcPrntScrn("exp_date : "); MdcPrintStr(mdcbufr);
307 strncpy(mdcbufr,dh.exp_time,10); mdcbufr[10]='\0';
308 MdcPrntScrn("exp_time : "); MdcPrintStr(mdcbufr);
309 MdcPrntScrn("views : %d\n",dh.views);
310 MdcPrntScrn("vols_added : %d\n",dh.vols_added);
311 MdcPrntScrn("start_field : %d\n",dh.start_field);
312 MdcPrntScrn("omax : %d\n",dh.omax);
313 MdcPrntScrn("omin : %d\n",dh.omin);
314 MdcPrntScrn("smax : %d\n",dh.smax);
315 MdcPrntScrn("smin : %d\n",dh.smin);
316 }
317
318
319 if (MDC_INFO) {
320 MdcPrntScrn("\n");
321 MdcPrintLine('=',MDC_FULL_LENGTH);
322 MdcPrntScrn("SPM - HEADER INTERPRETATION\n");
323 MdcPrintLine('-',MDC_HALF_LENGTH);
324 MdcPrntScrn("image {x} : %hd\n",imd.dim[1]);
325 MdcPrntScrn("image {y} : %hd\n",imd.dim[2]);
326 MdcPrntScrn("image {z} : %hd\n",imd.dim[3]);
327 MdcPrntScrn("voxel {x} : %+e\n",imd.pixdim[1]);
328 MdcPrntScrn("voxel {y} : %+e\n",imd.pixdim[2]);
329 MdcPrntScrn("voxel {z} : %+e\n",imd.pixdim[3]);
330 MdcPrntScrn("scaling : %+e\n",imd.spm_pix_rescale);
331 MdcPrntScrn("data type : %hd\n",imd.datatype);
332 MdcPrntScrn("offset : %+e\n",imd.avw_vox_offset);
333 MdcPrntScrn("origin : %hd %hd %hd\n",opt->origin_x
334 ,opt->origin_y
335 ,opt->origin_z);
336 MdcPrntScrn("description : "); MdcPrintStr(dh.descrip);
337 MdcPrintLine('=',MDC_FULL_LENGTH);
338 }
339
340 /* save the offset, valid for AVW / SPM / MRIcro Analyze files */
341 opt->offset = imd.avw_vox_offset;
342
343 /* update our FILEINFO structure */
344 MdcStringCopy(fi->study_descr,dh.descrip,80);
345 MdcStringCopy(fi->patient_id,dh.patient_id,10);
346 MdcStringCopy(fi->study_id,dh.scannum,10);
347
348 if (MDC_ECHO_ALIAS == MDC_YES) {
349 MdcEchoAliasName(fi); return(NULL);
350 }
351
352 memcpy(fi->dim,imd.dim,sizeof(imd.dim));
353 memcpy(fi->pixdim,imd.pixdim,sizeof(imd.pixdim));
354
355 fi->mwidth = (Uint32) imd.dim[1]; fi->mheight = (Uint32) imd.dim[2];
356 for ( number=1, i=3; i<=imd.dim[0]; i++)
357 number*=imd.dim[i];
358
359 if (number == 0) return("ANLZ No valid images specified");
360
361 fi->bits = imd.bitpix;
362
363 switch (imd.datatype) {
364 case MDC_ANLZ_DT_BINARY : fi->type=BIT1; fi->bits=8; break;
365 case MDC_ANLZ_DT_UNSIGNED_CHAR: fi->type=BIT8_U; fi->bits=8; break;
366 case MDC_ANLZ_DT_SIGNED_SHORT : fi->type=BIT16_S; fi->bits=16; break;
367 case MDC_ANLZ_DT_SIGNED_INT : fi->type=BIT32_S; fi->bits=32; break;
368 case MDC_ANLZ_DT_FLOAT : fi->type=FLT32; fi->bits=32; break;
369 case MDC_ANLZ_DT_COMPLEX :
370 return("ANLZ Datatype `complex' unsupported"); break;
371 case MDC_ANLZ_DT_DOUBLE : fi->type=FLT64; fi->bits=64; break;
372 case MDC_ANLZ_DT_RGB : fi->type=COLRGB; fi->bits=24;
373 fi->map=MDC_MAP_PRESENT; break;
374 case MDC_ANLZ_DT_ALL :
375 return("ANLZ Datatype `All' unsupported"); break;
376 default :
377 switch (fi->bits) {
378 case 1: fi->type=BIT1; break;
379 case 8: fi->type=BIT8_U; break;
380 case 16: fi->type=BIT16_S; break;
381 case 32: fi->type=BIT32_S; break; /* could be FLT32 as well */
382 default: MdcPrntWarn("ANLZ Unknown datatype");
383 }
384 }
385
386
387 /* preserve original path */
388 MdcMergePath(fi->ipath,fi->idir,fi->ifname);
389 if ((origpath=malloc(strlen(fi->ipath) + 1)) == NULL)
390 return("ANLZ Couldn't allocate original path");
391
392 strcpy(origpath,fi->ipath);
393 MdcSplitPath(fi->ipath,fi->idir,fi->ifname);
394
395
396 /* read the image file */
397 MdcCloseFile(fi->ifp);
398 MdcMergePath(fi->ipath,fi->idir,fi->ifname);
399 MdcSetExt(fi->ipath,"img");
400
401 /* check for compressed image file */
402 if (MdcFileExists(fi->ipath) == MDC_NO) {
403 MdcAddCompressionExt(fi->compression,fi->ipath);
404 if (MdcDecompressFile(fi->ipath) != MDC_OK) {
405 MdcFree(origpath); return("ANLZ Decompression image file failed");
406 }
407 WAS_COMPRESSED = MDC_YES;
408 }
409
410 if ( (fi->ifp=fopen(fi->ipath,"rb")) == NULL ) {
411 MdcFree(origpath); return("ANLZ Couldn't open image file");
412 }
413
414 if (WAS_COMPRESSED == MDC_YES) {
415 unlink(fi->ipath);
416 if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Reading Analyze:");
417 }
418
419 MdcSplitPath(fi->ipath,fi->idir,fi->ifname);
420
421 if (MDC_ANLZ_SPM == MDC_YES) {
422 /* interpret offset value from the header but we keep */
423 /* our precautions ... badly initialized headers */
424 long offsetu = (long)opt->offset;
425 if ((float)offsetu == opt->offset) fseek(fi->ifp,offsetu,SEEK_SET);
426 }
427
428 if (!MdcGetStructID(fi,number)) {
429 MdcFree(origpath); return("ANLZ Bad malloc IMG_DATA structs");
430 }
431
432
433 /* attempt to fill in orientation information */
434 switch (dh.orient) {
435 /* flipped, what's the meaning of flipped here ? */
436 case MDC_ANLZ_TRANS_UNFLIPPED:
437 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_TRANSAXIAL; break;
438 case MDC_ANLZ_CORON_UNFLIPPED:
439 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_CORONAL; break;
440 case MDC_ANLZ_SAGIT_UNFLIPPED:
441 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_SAGITTAL; break;
442 case MDC_ANLZ_TRANS_FLIPPED:
443 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_TRANSAXIAL; break;
444 case MDC_ANLZ_CORON_FLIPPED:
445 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_CORONAL; break;
446 case MDC_ANLZ_SAGIT_FLIPPED:
447 fi->pat_slice_orient=MDC_SUPINE_HEADFIRST_SAGITTAL; break;
448 }
449 strcpy(fi->pat_pos,MdcGetStrPatPos(fi->pat_slice_orient));
450 strcpy(fi->pat_orient,MdcGetStrPatOrient(fi->pat_slice_orient));
451
452 for ( i=0; i < fi->number; i++) {
453
454 if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL);
455
456 plane = i % fi->dim[3];
457
458 id = &fi->image[i];
459
460 id->width = fi->mwidth;
461 id->height = fi->mheight;
462 id->bits = fi->bits;
463 id->type = fi->type;
464 if (MDC_ANLZ_SPM) {
465 /* consider the scaling factor */
466 if (imd.spm_pix_rescale > 0.0 ) id->quant_scale = imd.spm_pix_rescale;
467 }
468 if (fi->pixdim[0] == 3.0 ) {
469 id->pixel_xsize = fi->pixdim[1];
470 id->pixel_ysize = fi->pixdim[2];
471 id->slice_width = fi->pixdim[3];
472 }else if (fi->pixdim[0] == 4.0 ) {
473 id->pixel_xsize = fi->pixdim[1];
474 id->pixel_ysize = fi->pixdim[2];
475 id->slice_width = fi->pixdim[3];
476 }else if ( (fi->pixdim[1] > 0.0) &&
477 (fi->pixdim[2] > 0.0) &&
478 (fi->pixdim[3] > 0.0) ) { /* we will try it anyway */
479 /* some don't fill in pixdim[0] */
480 /* for example PMOD (11-Apr-2000) */
481 id->pixel_xsize = fi->pixdim[1];
482 id->pixel_ysize = fi->pixdim[2];
483 id->slice_width = fi->pixdim[3];
484 fi->pixdim[0] = 3.0;
485 }else {
486 id->pixel_xsize = 1.0;
487 id->pixel_ysize = 1.0;
488 id->slice_width = 1.0;
489 }
490 id->slice_spacing = id->slice_width;
491 MdcFillImgPos(fi,i,plane,0.0);
492 MdcFillImgOrient(fi,i);
493
494 bytes = MdcPixels2Bytes(fi->mwidth*fi->mheight*fi->bits);
495 if ( (id->buf=MdcGetImgBuffer(bytes)) == NULL ) {
496 MdcFree(img8); MdcFree(origpath); return("ANLZ Bad malloc image buffer");
497 }
498
499 if (img8 != NULL) {
500 /* image from buffer */
501 memcpy(id->buf,img8+i*bytes,bytes);
502 }else{
503 /* image from file */
504 if (fread(id->buf,1,bytes,fi->ifp) != bytes ) {
505 err=MdcHandleTruncated(fi, i+1,MDC_YES);
506 if (err != NULL) { MdcFree(origpath); return(err); }
507 }
508
509 if (fi->truncated) break;
510 }
511
512 }
513
514 MdcFree(img8);
515
516 MdcCloseFile(fi->ifp);
517
518 /* check some final FILEINFO entries */
519 if (fi->dim[4] > 1) {
520 fi->acquisition_type = MDC_ACQUISITION_DYNAMIC;
521 /* fill in dynamic data struct */
522 if (!MdcGetStructDD(fi,(unsigned)fi->dim[4]))
523 return("ANLZ Couldn't malloc DYNAMIC_DATA structs");
524 for (f=0; f < fi->dynnr; f++) {
525 dd = &fi->dyndata[f];
526 dd->nr_of_slices = fi->dim[3];
527 dd->time_frame_duration = fi->pixdim[4];
528 dd->time_frame_start = f * dd->time_frame_duration;
529 }
530 }
531
532 /* restore original filename */
533 strcpy(fi->ipath,origpath);
534 MdcSplitPath(fi->ipath,fi->idir,fi->ifname);
535 MdcFree(origpath);
536
537 if (fi->truncated) return("ANLZ Truncated image file");
538
539 return NULL;
540 }
541
MdcWriteHeaderKey(FILEINFO * fi)542 int MdcWriteHeaderKey(FILEINFO *fi)
543 {
544 MDC_ANLZ_HEADER_KEY hk;
545 char *p = NULL;
546
547 memset(&hk,0,MDC_ANLZ_HK_SIZE);
548
549 hk.sizeof_hdr = MDC_ANLZ_HK_SIZE + MDC_ANLZ_IMD_SIZE + MDC_ANLZ_DH_SIZE;
550 sprintf(hk.data_type,"dsr");
551 MdcSplitPath(fi->opath,fi->odir,fi->ofname);
552 p = strrchr(fi->ofname,'.');
553 if (p != NULL) *p = '\0'; /* remove extension */
554 sprintf(hk.db_name,"%.17s",fi->ofname);
555 if (p != NULL) *p = '.'; /* add extension */
556 MdcMergePath(fi->opath,fi->odir,fi->ofname);
557 hk.extents=16384;
558 hk.session_error=0;
559 hk.regular='r';
560
561 MdcSWAP(hk.sizeof_hdr); MdcSWAP(hk.extents); MdcSWAP(hk.session_error);
562
563 fwrite((char *)&hk,1,MDC_ANLZ_HK_SIZE,fi->ofp);
564
565 if (ferror(fi->ofp)) return(MDC_NO);
566
567
568 return(MDC_YES);
569
570 }
571
MdcWriteImageDimension(FILEINFO * fi,MDC_SPMOPT * opt)572 int MdcWriteImageDimension(FILEINFO *fi, MDC_SPMOPT *opt)
573 {
574 MDC_ANLZ_IMAGE_DIMS imd;
575 float glmax=0., glmin=0.;
576 int i;
577
578 memset(&imd,0,MDC_ANLZ_IMD_SIZE);
579
580 strcpy(imd.avw_vox_units,"mm");
581
582 for (i=0; i <= fi->dim[0]; i++) imd.dim[i] = fi->dim[i];
583 for (i=0; i <= fi->pixdim[0]; i++) imd.pixdim[i] = fi->pixdim[i];
584
585 #if MDC_ALWAYS_SET_4D
586 /* set dummy 4th dimension (time) */
587 if (imd.dim[0] == 3) {
588 imd.dim[0] = 4;
589 imd.dim[4] = 1;
590 }
591 if (imd.pixdim[0] == 3.) {
592 imd.pixdim[0] = 4.;
593 imd.pixdim[4] = 0.;
594 }
595 #endif
596
597 #ifdef MDC_USE_SLICE_SPACING
598 if (fi->number > 1) imd.pixdim[3] = fi->image[0].slice_spacing;
599 #endif
600
601 imd.dim[1] = (Int16) fi->mwidth;
602 imd.dim[2] = (Int16) fi->mheight;
603
604 if (fi->map == MDC_MAP_PRESENT) {
605 /* colored */
606 imd.datatype = MDC_ANLZ_DT_RGB;
607 imd.bitpix = 24;
608 }else{
609 /* grayscale */
610 if (MDC_FORCE_INT != MDC_NO) {
611 switch (MDC_FORCE_INT) {
612 case BIT8_U : imd.datatype = MDC_ANLZ_DT_UNSIGNED_CHAR;
613 imd.bitpix = 8;
614 break;
615 case BIT16_S: imd.datatype = MDC_ANLZ_DT_SIGNED_SHORT;
616 imd.bitpix = 16;
617 break;
618 default : imd.datatype = MDC_ANLZ_DT_SIGNED_SHORT;
619 imd.bitpix = 16;
620 }
621 }else if (!(MDC_QUANTIFY || MDC_CALIBRATE)) {
622 if ( fi->diff_type ) {
623 imd.datatype = MDC_ANLZ_DT_SIGNED_SHORT;
624 imd.bitpix = 16;
625 }else{
626 switch ( fi->type ) {
627 case BIT8_U:
628 case BIT8_S: imd.datatype = MDC_ANLZ_DT_UNSIGNED_CHAR;
629 imd.bitpix = 8;
630 break;
631 case BIT16_U:
632 case BIT16_S: imd.datatype = MDC_ANLZ_DT_SIGNED_SHORT;
633 imd.bitpix = 16;
634 break;
635 case BIT64_U:
636 case BIT64_S:
637 case BIT32_U:
638 case BIT32_S: imd.datatype = MDC_ANLZ_DT_SIGNED_INT;
639 imd.bitpix = 32;
640 break;
641 case FLT32: imd.datatype = MDC_ANLZ_DT_FLOAT;
642 imd.bitpix = 32;
643 break;
644 case FLT64: imd.datatype = MDC_ANLZ_DT_DOUBLE;
645 imd.bitpix = 64;
646 break;
647 }
648 }
649 }else{
650 if (MDC_ANLZ_SPM == MDC_YES) { /* BIT16_S with scaling factor */
651 imd.datatype = MDC_ANLZ_DT_SIGNED_SHORT;
652 imd.bitpix = 16;
653 }else{
654 imd.datatype = MDC_ANLZ_DT_FLOAT;
655 imd.bitpix = 32;
656 }
657 }
658 }
659
660 /* find and set max/min values */
661 for (i = 0; i < fi->number; i++) {
662 IMG_DATA *id = &fi->image[i];
663 if (id->rescaled == MDC_YES) {
664 if (i == 0) { /* init values */
665 glmax = id->rescaled_max;
666 glmin = id->rescaled_min;
667 }else{ /* get max/min */
668 glmax = (id->rescaled_max > glmax) ? id->rescaled_max : glmax;
669 glmin = (id->rescaled_min < glmin) ? id->rescaled_min : glmin;
670 }
671 }else{
672 if (i == 0) { /* init values */
673 glmax = id->max;
674 glmin = id->min;
675 }else{ /* get max/min */
676 glmax = (id->max > glmax) ? id->max : glmax;
677 glmin = (id->min < glmin) ? id->min : glmin;
678 }
679 }
680 }
681
682 imd.glmax = (Int32) glmax;
683 imd.glmin = (Int32) glmin;
684
685 imd.avw_cal_max = fi->qglmax;
686 imd.avw_cal_min = fi->qglmin;
687
688 /* thinking about SPM */
689 if (imd.pixdim[0] <= 0.0 || imd.pixdim[0] >= (float)MDC_ANLZ_MAX_DIMS) {
690 imd.pixdim[0]=3.;
691 imd.pixdim[1]=1.;
692 imd.pixdim[2]=1.;
693 imd.pixdim[3]=1.;
694 }
695
696 if (opt != NULL) imd.avw_vox_offset = opt->offset;
697
698 if (MDC_ANLZ_SPM == MDC_YES) { /* the scaling factor */
699 if (fi->image[0].rescaled)
700 imd.spm_pix_rescale=(float)fi->image[0].rescaled_fctr;
701 /* did rescale over all images -> all images same factor */
702 }else{
703 imd.spm_pix_rescale=1.;
704 }
705
706 /* swap the data if necessary */
707 for (i=0; i<MDC_ANLZ_MAX_DIMS; i++) {
708 MdcSWAP(imd.dim[i]); MdcSWAP(imd.pixdim[i]);
709 }
710 MdcSWAP(imd.unused1);
711 MdcSWAP(imd.avw_vox_offset);
712 MdcSWAP(imd.spm_pix_rescale);
713 MdcSWAP(imd.funused1); MdcSWAP(imd.funused2);
714 MdcSWAP(imd.avw_cal_max); MdcSWAP(imd.avw_cal_min);
715 MdcSWAP(imd.datatype); MdcSWAP(imd.bitpix);
716 MdcSWAP(imd.compressed); MdcSWAP(imd.verified);
717 MdcSWAP(imd.glmax); MdcSWAP(imd.glmin);
718
719 fwrite((char *)&imd,1,MDC_ANLZ_IMD_SIZE,fi->ofp);
720
721 if (ferror(fi->ofp)) return(MDC_NO);
722
723 return(MDC_YES);
724
725 }
726
MdcWriteDataHistory(FILEINFO * fi,MDC_SPMOPT * opt)727 int MdcWriteDataHistory(FILEINFO *fi, MDC_SPMOPT *opt)
728 {
729 MDC_ANLZ_DATA_HIST dh;
730
731 memset(&dh,0,MDC_ANLZ_DH_SIZE);
732
733 sprintf(dh.descrip,"%.35s",fi->study_descr);
734 sprintf(dh.scannum,"%.9s",fi->study_id);
735 sprintf(dh.patient_id,"%.9s",fi->patient_id);
736 sprintf(dh.generated,"%.9s",MDC_PRGR);
737
738 switch (fi->pat_slice_orient) {
739 case MDC_SUPINE_HEADFIRST_TRANSAXIAL :
740 case MDC_PRONE_HEADFIRST_TRANSAXIAL :
741 case MDC_DECUBITUS_RIGHT_HEADFIRST_TRANSAXIAL:
742 case MDC_DECUBITUS_LEFT_HEADFIRST_TRANSAXIAL :
743 case MDC_SUPINE_FEETFIRST_TRANSAXIAL :
744 case MDC_PRONE_FEETFIRST_TRANSAXIAL :
745 case MDC_DECUBITUS_RIGHT_FEETFIRST_TRANSAXIAL:
746 case MDC_DECUBITUS_LEFT_FEETFIRST_TRANSAXIAL :
747 dh.orient = MDC_ANLZ_TRANS_UNFLIPPED;
748 break;
749 case MDC_SUPINE_HEADFIRST_CORONAL :
750 case MDC_PRONE_HEADFIRST_CORONAL :
751 case MDC_DECUBITUS_RIGHT_HEADFIRST_CORONAL :
752 case MDC_DECUBITUS_LEFT_HEADFIRST_CORONAL :
753 case MDC_SUPINE_FEETFIRST_CORONAL :
754 case MDC_PRONE_FEETFIRST_CORONAL :
755 case MDC_DECUBITUS_RIGHT_FEETFIRST_CORONAL :
756 case MDC_DECUBITUS_LEFT_FEETFIRST_CORONAL :
757 dh.orient = MDC_ANLZ_CORON_UNFLIPPED;
758 break;
759
760 case MDC_SUPINE_HEADFIRST_SAGITTAL :
761 case MDC_PRONE_HEADFIRST_SAGITTAL :
762 case MDC_DECUBITUS_RIGHT_HEADFIRST_SAGITTAL :
763 case MDC_DECUBITUS_LEFT_HEADFIRST_SAGITTAL :
764 case MDC_SUPINE_FEETFIRST_SAGITTAL :
765 case MDC_PRONE_FEETFIRST_SAGITTAL :
766 case MDC_DECUBITUS_RIGHT_FEETFIRST_SAGITTAL :
767 case MDC_DECUBITUS_LEFT_FEETFIRST_SAGITTAL :
768 dh.orient = MDC_ANLZ_SAGIT_UNFLIPPED;
769 break;
770 }
771
772 if (opt != NULL) {
773 MdcSWAP(opt->origin_x); memcpy(&dh.originator[0],&opt->origin_x,2);
774 MdcSWAP(opt->origin_y); memcpy(&dh.originator[2],&opt->origin_y,2);
775 MdcSWAP(opt->origin_z); memcpy(&dh.originator[4],&opt->origin_z,2);
776 }
777
778 fwrite((char *)&dh,1,MDC_ANLZ_DH_SIZE,fi->ofp);
779
780 if (ferror(fi->ofp)) return(MDC_NO);
781
782 return(MDC_YES);
783
784 }
785
MdcWriteImagesData(FILEINFO * fi)786 char *MdcWriteImagesData(FILEINFO *fi)
787 {
788 double pval;
789 Uint8 grval;
790 Uint32 i, FREE;
791 Uint32 size, n, nr;
792 Uint16 type;
793 Uint8 *buf, *maxbuf;
794 Int8 saved_norm_over_frames=MDC_NORM_OVER_FRAMES;
795 IMG_DATA *id;
796
797 for (i=fi->number; i>0; i-- ) {
798
799 if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL);
800
801 nr = fi->number - i; /* normal planes */
802
803 id = &fi->image[nr];
804
805 buf = id->buf; FREE = MDC_NO;
806 type = id->type;
807
808 if (fi->map != MDC_MAP_PRESENT) {
809 /* grayscale */
810 if (MDC_FORCE_INT != MDC_NO) {
811 if (MDC_ANLZ_SPM) MDC_NORM_OVER_FRAMES = MDC_NO;
812 switch (MDC_FORCE_INT) {
813 case BIT8_U : buf = MdcGetImgBIT8_U(fi,nr);
814 type=BIT8_U; FREE=MDC_YES;
815 break;
816 case BIT16_S: buf = MdcGetImgBIT16_S(fi,nr);
817 type=BIT16_S; FREE=MDC_YES;
818 break;
819 default : buf = MdcGetImgBIT16_S(fi,nr);
820 type=BIT16_S; FREE=MDC_YES;
821 }
822 if (MDC_ANLZ_SPM) MDC_NORM_OVER_FRAMES = saved_norm_over_frames;
823 }else if (!(MDC_QUANTIFY || MDC_CALIBRATE)) {
824 if ( fi->diff_type ) {
825 switch (id->type) {
826 case BIT16_S: buf = id->buf;
827 type = BIT16_S; FREE=MDC_NO;
828 break;
829 default : buf = MdcGetImgBIT16_S(fi,nr);
830 type = BIT16_S; FREE=MDC_YES;
831 break;
832 }
833 }else{
834 switch (id->type) {
835 case BIT8_S: buf = MdcGetImgBIT8_U(fi,nr);
836 type=BIT8_U ; FREE=MDC_YES;
837 break;
838 case BIT16_U: buf = MdcGetImgBIT16_S(fi,nr);
839 type=BIT16_S; FREE=MDC_YES;
840 break;
841 case BIT32_U: buf = MdcGetImgBIT32_S(fi,nr);
842 type=BIT32_S; FREE=MDC_YES;
843 break;
844 case BIT64_S:
845 case BIT64_U: buf = MdcGetImgBIT32_S(fi,nr);
846 type=BIT32_S; FREE=MDC_YES;
847 break;
848 }
849 }
850 }else{
851 if (MDC_ANLZ_SPM == MDC_YES) {
852 /* using the global scale factor <=> normalize over ALL images! */
853 /* so all images have the same scale factor */
854 MDC_NORM_OVER_FRAMES=MDC_NO;
855 buf = MdcGetImgBIT16_S(fi,nr);
856 type = BIT16_S; FREE=MDC_YES;
857 MDC_NORM_OVER_FRAMES=saved_norm_over_frames;
858 }else{
859 buf = MdcGetImgFLT32(fi,nr);
860 type=FLT32; FREE=MDC_YES;
861 }
862 }
863 }
864
865 if (buf == NULL) return("ANLZ Bad malloc image buffer");
866
867
868 if (fi->diff_size) {
869
870 maxbuf = MdcGetResizedImage(fi, buf, type, nr);
871
872 if (maxbuf == NULL) return("ANLZ Bad malloc maxbuf");
873
874 if (FREE) MdcFree(buf);
875
876 FREE = MDC_YES;
877
878 }else{
879
880 maxbuf = buf;
881
882 }
883
884 size = fi->mwidth * fi->mheight * MdcType2Bytes(type);
885
886 if (fi->type == COLRGB) {
887 /* true color */
888 if (fwrite(maxbuf,1,size,fi->ofp) != size)
889 return("ANLZ Bad write RGB buffer");
890 }else{
891 for (n=0; n < size; n += MdcType2Bytes(type)) {
892 /* indexed */
893 pval = MdcGetDoublePixel((Uint8 *)&maxbuf[n],type);
894 if (fi->map == MDC_MAP_PRESENT) {
895 /* colored */
896 grval = (Uint8)pval;
897 fwrite(&fi->palette[grval * 3 + 0], 1, 1, fi->ofp); /* red */
898 fwrite(&fi->palette[grval * 3 + 1], 1, 1, fi->ofp); /* green */
899 fwrite(&fi->palette[grval * 3 + 2], 1, 1, fi->ofp); /* blue */
900 if (ferror(fi->ofp))
901 return("ANLZ Bad write colored pixel");
902 }else{
903 /* grayscale */
904 if (!MdcWriteDoublePixel(pval,type,fi->ofp))
905 return("ANLZ Bad write image pixel");
906 }
907 }
908 }
909
910 if (FREE) MdcFree(maxbuf);
911
912 if (ferror(fi->ofp)) return("ANLZ Bad writing of images");
913
914 }
915
916 return NULL;
917
918 }
919
MdcGetSpmOpt(FILEINFO * fi,MDC_SPMOPT * opt)920 void MdcGetSpmOpt(FILEINFO *fi, MDC_SPMOPT *opt)
921 {
922 if (INIT_SPMOPT == MDC_YES) {
923 opt->origin_x = 0;
924 opt->origin_y = 0;
925 opt->origin_z = 0;
926 opt->offset = 0.;
927 INIT_SPMOPT = MDC_NO;
928 }
929
930 if (MDC_FILE_STDIN == MDC_YES) return; /* stdin already in use */
931
932 MdcPrintLine('-',MDC_FULL_LENGTH);
933 MdcPrntScrn("\tSPM OPTIONS\t\tORIG FILE: %s\n",fi->ifname);
934 MdcPrintLine('-',MDC_FULL_LENGTH);
935 MdcPrntScrn("\n\tThe origin values must be an Int16 value");
936 MdcPrntScrn("\n\tThere is NO check performed on the input!\n");
937 MdcPrntScrn("\n\tOrigin X [%d]? ",opt->origin_x);
938 if (!MdcPutDefault(mdcbufr)) opt->origin_x = (Int16)atoi(mdcbufr);
939 MdcPrntScrn("\n\tOrigin Y [%d]? ",opt->origin_y);
940 if (!MdcPutDefault(mdcbufr)) opt->origin_y = (Int16)atoi(mdcbufr);
941 MdcPrntScrn("\n\tOrigin Z [%d]? ",opt->origin_z);
942 if (!MdcPutDefault(mdcbufr)) opt->origin_z = (Int16)atoi(mdcbufr);
943 /* MARK: skip asking about offset */
944 /* MdcPrntScrn("\n\tOffset [%+e]? ",opt->offset); */
945 /* if (!MdcPutDefault(mdcbufr)) opt->offset = (float)atof(mdcbufr); */
946 MdcPrntScrn("\n");
947 MdcPrintLine('-',MDC_FULL_LENGTH);
948 }
949
MdcWriteANLZ(FILEINFO * fi)950 const char *MdcWriteANLZ(FILEINFO *fi)
951 {
952 MDC_SPMOPT *opt = &spmopt;
953 char tmpfname[MDC_MAX_PATH + 1];
954 const char *msg;
955
956 MDC_FILE_ENDIAN = MDC_WRITE_ENDIAN;
957
958 /* user wanted to supply some parameters */
959 if ((MDC_ANLZ_OPTIONS == MDC_YES) && (XMDC_GUI == MDC_NO)) {
960 MdcGetSpmOpt(fi,opt);
961 }else {
962 /* set default origin to image centre of middle slice */
963 opt->origin_x = (Int16)((fi->dim[1] + 1)/2);
964 opt->origin_y = (Int16)((fi->dim[2] + 1)/2);
965 opt->origin_z = (Int16)((fi->dim[3] + 1)/2);
966 opt->offset = 0.;
967 }
968
969 /* header and image separate, rescaled stuff very important */
970 /* so we will write the images first ! */
971
972 /* get filename: no longer with truncation */
973 /* SPM, PMOD etc don't rely on db_name[18]) */
974 if (XMDC_GUI == MDC_YES) {
975 strcpy(tmpfname,fi->opath);
976 }else{
977 if (MDC_ALIAS_NAME == MDC_YES) {
978 MdcAliasName(fi,tmpfname);
979 }else{
980 strcpy(tmpfname,fi->ifname);
981 }
982 MdcDefaultName(fi,MDC_FRMT_ANLZ,fi->ofname,tmpfname);
983 }
984
985 if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Writing Analyze:");
986
987 if (MDC_VERBOSE) MdcPrntMesg("ANLZ Writing <%s> & <.img> ...",fi->ofname);
988
989 /* writing images */
990 if (XMDC_GUI == MDC_YES) {
991 fi->ofname[0]='\0'; MdcNewExt(fi->ofname,tmpfname,"img");
992 }else{
993 MdcNewName(fi->ofname,tmpfname,"img");
994 }
995
996 if (MDC_FILE_STDOUT == MDC_YES) {
997 /* send image data to stdout (1>stdout) */
998 fi->ofp = stdout;
999 }else{
1000 if (MdcKeepFile(fi->ofname))
1001 return("ANLZ Image file exists!!");
1002 if ( (fi->ofp=fopen(fi->ofname,"wb")) == NULL )
1003 return ("ANLZ Couldn't open image file");
1004 }
1005
1006 msg = MdcWriteImagesData(fi);
1007 if (msg != NULL) return(msg);
1008
1009 MdcCloseFile(fi->ofp);
1010
1011 /* writing header with rescaled stuff */
1012 if (XMDC_GUI == MDC_YES) {
1013 strcpy(fi->ofname,tmpfname);
1014 }else{
1015 MdcDefaultName(fi,MDC_FRMT_ANLZ,fi->ofname,tmpfname);
1016 }
1017
1018 if (MDC_FILE_STDOUT == MDC_YES) {
1019 /* send header to stderr (2>stderr) */
1020 fi->ofp = stderr;
1021 }else{
1022 if (MdcKeepFile(fi->ofname))
1023 return("ANLZ Header file exists!!");
1024 if ( (fi->ofp=fopen(fi->ofname,"wb")) == NULL )
1025 return("ANLZ Couldn't open header file");
1026 }
1027
1028 if ( !MdcWriteHeaderKey(fi) )
1029 return("ANLZ Bad write HeaderKey struct");
1030
1031 if ( !MdcWriteImageDimension(fi, opt) )
1032 return("ANLZ Bad write ImageDimension struct");
1033
1034 if ( !MdcWriteDataHistory(fi, opt) )
1035 return("ANLZ Bad write DataHistory struct");
1036
1037 MdcCheckQuantitation(fi);
1038
1039 MdcCloseFile(fi->ofp);
1040
1041 return(NULL);
1042
1043 }
1044