1 #define _NIFTI1_IO_C_
2
3 #include "nifti1_io.h" /* typedefs, prototypes, macros, etc. */
4
5 /*****===================================================================*****/
6 /***** Sample functions to deal with NIFTI-1 and ANALYZE files *****/
7 /*****...................................................................*****/
8 /***** This code is released to the public domain. *****/
9 /*****...................................................................*****/
10 /***** Author: Robert W Cox, SSCC/DIRP/NIMH/NIH/DHHS/USA/EARTH *****/
11 /***** Date: August 2003 *****/
12 /*****...................................................................*****/
13 /***** Neither the National Institutes of Health (NIH), nor any of its *****/
14 /***** employees imply any warranty of usefulness of this software for *****/
15 /***** any purpose, and do not assume any liability for damages, *****/
16 /***** incidental or otherwise, caused by any use of this document. *****/
17 /*****===================================================================*****/
18
19 /** \file nifti1_io.c
20 \brief main collection of nifti1 i/o routines
21 - written by Bob Cox, SSCC NIMH
22 - revised by Mark Jenkinson, FMRIB
23 - revised by Rick Reynolds, SSCC, NIMH
24 - revised by Kate Fissell, University of Pittsburgh
25
26 The library history can be viewed via "nifti_tool -nifti_hist".
27 <br>The library version can be viewed via "nifti_tool -nifti_ver".
28 */
29
30 /*! global history and version strings, for printing */
31 static char * gni_history[] =
32 {
33 "----------------------------------------------------------------------\n"
34 "history (of nifti library changes):\n"
35 "\n",
36 "0.0 August, 2003 [rwcox]\n"
37 " (Robert W Cox of the National Institutes of Health, SSCC/DIRP/NIMH)\n"
38 " - initial version\n"
39 "\n",
40 "0.1 July/August, 2004 [Mark Jenkinson]\n"
41 " (FMRIB Centre, University of Oxford, UK)\n"
42 " - Mainly adding low-level IO and changing things to allow gzipped\n"
43 " files to be read and written\n"
44 " - Full backwards compatability should have been maintained\n"
45 "\n",
46 "0.2 16 Nov 2004 [rickr]\n"
47 " (Rick Reynolds of the National Institutes of Health, SSCC/DIRP/NIMH)\n"
48 " - included Mark's changes in the AFNI distribution (including znzlib/)\n"
49 " (HAVE_ZLIB is commented out for the standard distribution)\n"
50 " - modified nifti_validfilename() and nifti_makebasename()\n"
51 " - added nifti_find_file_extension()\n"
52 "\n",
53 "0.3 3 Dec 2004 [rickr]\n"
54 " - note: header extensions are not yet checked for\n"
55 " - added formatted history as global string, for printing\n"
56 " - added nifti_disp_lib_hist(), to display the nifti library history\n"
57 " - added nifti_disp_lib_version(), to display the nifti library history\n",
58 " - re-wrote nifti_findhdrname()\n"
59 " o used nifti_find_file_extension()\n"
60 " o changed order of file tests (default is .nii, depends on input)\n"
61 " o free hdrname on failure\n"
62 " - made similar changes to nifti_findimgname()\n"
63 " - check for NULL return from nifti_findhdrname() calls\n",
64 " - removed most of ERREX() macros\n"
65 " - modified nifti_image_read()\n"
66 " o added debug info and error checking (on gni_debug > 0, only)\n"
67 " o fail if workingname is NULL\n"
68 " o check for failure to open header file\n"
69 " o free workingname on failure\n"
70 " o check for failure of nifti_image_load()\n"
71 " o check for failure of nifti_convert_nhdr2nim()\n",
72 " - changed nifti_image_load() to int, and check nifti_read_buffer return\n"
73 " - changed nifti_read_buffer() to fail on short read, and to count float\n"
74 " fixes (to print on debug)\n"
75 " - changed nifti_image_infodump to print to stderr\n"
76 " - updated function header comments, or moved comments above header\n"
77 " - removed const keyword\n"
78 " - added LNI_FERR() macro for error reporting on input files\n"
79 "\n",
80 "0.4 10 Dec 2004 [rickr] - added header extensions\n"
81 " - in nifti1_io.h:\n"
82 " o added num_ext and ext_list to the definition of nifti_image\n"
83 " o made many functions static (more to follow)\n"
84 " o added LNI_MAX_NIA_EXT_LEN, for max nifti_type 3 extension length\n",
85 " - added __DATE__ to version output in nifti_disp_lib_version()\n"
86 " - added nifti_disp_matrix_orient() to print orientation information\n"
87 " - added '.nia' as a valid file extension in nifti_find_file_extension()\n"
88 " - added much more debug output\n"
89 " - in nifti_image_read(), in the case of an ASCII header, check for\n"
90 " extensions after the end of the header\n",
91 " - added nifti_read_extensions() function\n"
92 " - added nifti_read_next_extension() function\n"
93 " - added nifti_add_exten_to_list() function\n"
94 " - added nifti_check_extension() function\n"
95 " - added nifti_write_extensions() function\n"
96 " - added nifti_extension_size() function\n"
97 " - in nifti_set_iname_offest():\n"
98 " o adjust offset by the extension size and the extender size\n",
99 " o fixed the 'ceiling modulo 16' computation\n"
100 " - in nifti_image_write_hdr_img2(): \n"
101 " o added extension writing\n"
102 " o check for NULL return from nifti_findimgname()\n"
103 " - include number of extensions in nifti_image_to_ascii() output\n"
104 " - in nifti_image_from_ascii():\n"
105 " o return bytes_read as a parameter, computed from the final spos\n"
106 " o extract num_ext from ASCII header\n"
107 "\n",
108 "0.5 14 Dec 2004 [rickr] - added sub-brick reading functions\n"
109 " - added nifti_brick_list type to nifti1_io.h, along with new prototypes\n"
110 " - added main nifti_image_read_bricks() function, with description\n"
111 " - added nifti_image_load_bricks() - library function (requires nim)\n"
112 " - added valid_nifti_brick_list() - library function\n"
113 " - added free_NBL() - library function\n",
114 " - added update_nifti_image_for_brick_list() for dimension update\n"
115 " - added nifti_load_NBL_bricks(), nifti_alloc_NBL_mem(),\n"
116 " nifti_copynsort() and force_positive() (static functions)\n"
117 " - in nifti_image_read(), check for failed load only if read_data is set\n"
118 " - broke most of nifti_image_load() into nifti_image_load_prep()\n"
119 "\n",
120 "0.6 15 Dec 2004 [rickr] - added sub-brick writing functionality\n"
121 " - in nifti1_io.h, removed znzlib directory from include - all nifti\n"
122 " library files are now under the nifti directory\n"
123 " - nifti_read_extensions(): print no offset warning for nifti_type 3\n"
124 " - nifti_write_all_data():\n"
125 " o pass nifti_brick_list * NBL, for optional writing\n"
126 " o if NBL, write each sub-brick, sequentially\n",
127 " - nifti_set_iname_offset(): case 1 must have sizeof() cast to int\n"
128 " - pass NBL to nifti_image_write_hdr_img2(), and allow NBL or data\n"
129 " - added nifti_image_write_bricks() wrapper for ...write_hdr_img2()\n"
130 " - included compression abilities\n"
131 "\n",
132 "0.7 16 Dec 2004 [rickr] - minor changes to extension reading\n"
133 "\n",
134 "0.8 21 Dec 2004 [rickr] - restrict extension reading, and minor changes\n"
135 " - in nifti_image_read(), compute bytes for extensions (see remaining)\n"
136 " - in nifti_read_extensions(), pass 'remain' as space for extensions,\n"
137 " pass it to nifti_read_next_ext(), and update for each one read \n"
138 " - in nifti_check_extension(), require (size <= remain)\n",
139 " - in update_nifti_image_brick_list(), update nvox\n"
140 " - in nifti_image_load_bricks(), make explicit check for nbricks <= 0\n"
141 " - in int_force_positive(), check for (!list)\n"
142 " - in swap_nifti_header(), swap sizeof_hdr, and reorder to struct order\n"
143 " - change get_filesize functions to signed ( < 0 is no file or error )\n",
144 " - in nifti_validfilename(), lose redundant (len < 0) check\n"
145 " - make print_hex_vals() static\n"
146 " - in disp_nifti_1_header, restrict string field widths\n"
147 "\n",
148 "0.9 23 Dec 2004 [rickr] - minor changes\n"
149 " - broke ASCII header reading out of nifti_image_read(), into new\n"
150 " functions has_ascii_header() and read_ascii_image()\n",
151 " - check image_read failure and znzseek failure\n"
152 " - altered some debug output\n"
153 " - nifti_write_all_data() now returns an int\n"
154 "\n",
155 "0.10 29 Dec 2004 [rickr]\n"
156 " - renamed nifti_valid_extension() to nifti_check_extension()\n"
157 " - added functions nifti_makehdrname() and nifti_makeimgname()\n"
158 " - added function valid_nifti_extensions()\n"
159 " - in nifti_write_extensions(), check for validity before writing\n",
160 " - rewrote nifti_image_write_hdr_img2():\n"
161 " o set write_data and leave_open flags from write_opts\n"
162 " o add debug print statements\n"
163 " o use nifti_write_ascii_image() for the ascii case\n"
164 " o rewrote the logic of all cases to be easier to follow\n",
165 " - broke out code as nifti_write_ascii_image() function\n"
166 " - added debug to top-level write functions, and free the znzFile\n"
167 " - removed unused internal function nifti_image_open()\n"
168 "\n",
169 "0.11 30 Dec 2004 [rickr] - small mods\n"
170 " - moved static function prototypes from header to C file\n"
171 " - free extensions in nifti_image_free()\n"
172 "\n",
173 "1.0 07 Jan 2005 [rickr] - INITIAL RELEASE VERSION\n"
174 " - added function nifti_set_filenames()\n"
175 " - added function nifti_read_header()\n"
176 " - added static function nhdr_looks_good()\n"
177 " - added static function need_nhdr_swap()\n"
178 " - exported nifti_add_exten_to_list symbol\n",
179 " - fixed #bytes written in nifti_write_extensions()\n"
180 " - only modify offset if it is too small (nifti_set_iname_offset)\n"
181 " - added nifti_type 3 to nifti_makehdrname and nifti_makeimgname\n"
182 " - added function nifti_set_filenames()\n"
183 "\n",
184 "1.1 07 Jan 2005 [rickr]\n"
185 " - in nifti_read_header(), swap if needed\n"
186 "\n",
187 "1.2 07 Feb 2005 [kate fissell c/o rickr] \n"
188 " - nifti1.h: added doxygen comments for main struct and #define groups\n"
189 " - nifti1_io.h: added doxygen comments for file and nifti_image struct\n"
190 " - nifti1_io.h: added doxygen comments for file and some functions\n"
191 " - nifti1_io.c: changed nifti_copy_nim_info to use memcpy\n"
192 "\n",
193 "1.3 09 Feb 2005 [rickr]\n"
194 " - nifti1.h: added doxygen comments for extension structs\n"
195 " - nifti1_io.h: put most #defines in #ifdef _NIFTI1_IO_C_ block\n"
196 " - added a doxygen-style description to every exported function\n"
197 " - added doxygen-style comments within some functions\n"
198 " - re-exported many znzFile functions that I had made static\n"
199 " - re-added nifti_image_open (sorry, Mark)\n"
200 " - every exported function now has 'nifti' in the name (19 functions)\n",
201 " - made sure every alloc() has a failure test\n"
202 " - added nifti_copy_extensions function, for use in nifti_copy_nim_info\n"
203 " - nifti_is_gzfile: added initial strlen test\n"
204 " - nifti_set_filenames: added set_byte_order parameter option\n"
205 " (it seems appropriate to set the BO when new files are associated)\n"
206 " - disp_nifti_1_header: prints to stdout (a.o.t. stderr), with fflush\n"
207 "\n",
208 "1.4 23 Feb 2005 [rickr] - sourceforge merge\n"
209 " - merged into the nifti_io CVS directory structure at sourceforge.net\n"
210 " - merged in 4 changes by Mark, and re-added his const keywords\n"
211 " - cast some pointers to (void *) for -pedantic compile option\n"
212 " - added nifti_free_extensions()\n"
213 "\n",
214 "1.5 02 Mar 2005 [rickr] - started nifti global options\n"
215 " - gni_debug is now g_opts.debug\n"
216 " - added validity check parameter to nifti_read_header\n"
217 " - need_nhdr_swap no longer does test swaps on the stack\n"
218 "\n",
219 "1.6 05 April 2005 [rickr] - validation and collapsed_image_read\n"
220 " - added nifti_read_collapsed_image(), an interface for reading partial\n"
221 " datasets, specifying a subset of array indices\n"
222 " - for read_collapsed_image, added static functions: rci_read_data(),\n"
223 " rci_alloc_mem(), and make_pivot_list()\n",
224 " - added nifti_nim_is_valid() to check for consistency (more to do)\n"
225 " - added nifti_nim_has_valid_dims() to do many dimensions tests\n"
226 "\n",
227 "1.7 08 April 2005 [rickr]\n"
228 " - added nifti_update_dims_from_array() - to update dimensions\n"
229 " - modified nifti_makehdrname() and nifti_makeimgname():\n"
230 " if prefix has a valid extension, use it (else make one up)\n"
231 " - added nifti_get_intlist - for making an array of ints\n"
232 " - fixed init of NBL->bsize in nifti_alloc_NBL_mem() {thanks, Bob}\n"
233 "\n",
234 "1.8 14 April 2005 [rickr]\n"
235 " - added nifti_set_type_from_names(), for nifti_set_filenames()\n"
236 " (only updates type if number of files does not match it)\n"
237 " - added is_valid_nifti_type(), just to be sure\n"
238 " - updated description of nifti_read_collapsed_image() for *data change\n"
239 " (if *data is already set, assume memory exists for results)\n"
240 " - modified rci_alloc_mem() to allocate only if *data is NULL\n"
241 "\n",
242 "1.9 19 April 2005 [rickr]\n"
243 " - added extension codes NIFTI_ECODE_COMMENT and NIFTI_ECODE_XCEDE\n"
244 " - added nifti_type codes NIFTI_MAX_ECODE and NIFTI_MAX_FTYPE\n"
245 " - added nifti_add_extension() {exported}\n"
246 " - added nifti_fill_extension() as a static function\n"
247 " - added nifti_is_valid_ecode() {exported}\n",
248 " - nifti_type values are now NIFTI_FTYPE_* file codes\n"
249 " - in nifti_read_extensions(), decrement 'remain' by extender size, 4\n"
250 " - in nifti_set_iname_offset(), case 1, update if offset differs\n"
251 " - only output '-d writing nifti file' if debug > 1\n"
252 "\n",
253 "1.10 10 May 2005 [rickr]\n"
254 " - files are read using ZLIB only if they end in '.gz'\n"
255 "\n",
256 "1.11 12 August 2005 [kate fissell]\n"
257 " - Kate's 0.2 release packaging, for sourceforge\n"
258 "\n",
259 "1.12 17 August 2005 [rickr] - comment (doxygen) updates\n"
260 " - updated comments for most functions (2 updates from Cinly Ooi)\n"
261 " - added nifti_type_and_names_match()\n"
262 "\n",
263 "1.12a 24 August 2005 [rickr] - remove all tabs from Clibs/*/*.[ch]\n",
264 "1.12b 25 August 2005 [rickr] - changes by Hans Johnson\n",
265 "1.13 25 August 2005 [rickr]\n",
266 " - finished changes by Hans for Insight\n"
267 " - added const in all appropraite parameter locations (30-40)\n"
268 " (any pointer referencing data that will not change)\n"
269 " - shortened all string constants below 509 character limit\n"
270 "1.14 28 October 2005 [HJohnson]\n",
271 " - use nifti_set_filenames() in nifti_convert_nhdr2nim()\n"
272 "1.15 02 November 2005 [rickr]\n",
273 " - added skip_blank_ext to nifti_global_options\n"
274 " - added nifti_set_skip_blank_ext(), to set option\n"
275 " - if skip_blank_ext and no extensions, do not read/write extender\n"
276 "1.16 18 November 2005 [rickr]\n",
277 " - removed any test or access of dim[i], i>dim[0]\n"
278 " - do not set pixdim for collapsed dims to 1.0, leave them as they are\n"
279 " - added magic and dim[i] tests in nifti_hdr_looks_good()\n"
280 " - added 2 size_t casts\n"
281 "1.17 22 November 2005 [rickr]\n",
282 " - in hdr->nim, for i > dim[0], pass 0 or 1, else set to 1\n"
283 "1.18 02 March 2006 [rickr]\n",
284 " - in nifti_alloc_NBL_mem(), fixed nt=0 case from 1.17 change\n"
285 "1.19 23 May 2006 [HJohnson,rickr]\n",
286 " - nifti_write_ascii_image(): free(hstr)\n"
287 " - nifti_copy_extensions(): clear num_ext and ext_list\n"
288 "1.20 27 Jun 2006 [rickr]\n",
289 " - nifti_findhdrname(): fixed assign of efirst to match stated logic\n"
290 " (problem found by Atle Bjørnerud)\n"
291 "1.21 5 Sep 2006 [rickr] update for nifticlib-0.4 release\n",
292 " - was reminded to actually add nifti_set_skip_blank_ext()\n"
293 " - init g_opts.skip_blank_ext to 0\n"
294 "----------------------------------------------------------------------\n"
295 };
296 static char gni_version[] = "nifti library version 1.21 5 Sep, 2006)";
297
298 /*! global nifti options structure */
299 static nifti_global_options g_opts = { 1, 0 };
300
301 /*---------------------------------------------------------------------------*/
302 /* prototypes for internal functions - not part of exported library */
303
304 /* extension routines */
305 static int nifti_read_extensions( nifti_image *nim, znzFile fp, int remain );
306 static int nifti_read_next_extension( nifti1_extension * nex, nifti_image *nim, int remain, znzFile fp );
307 static int nifti_check_extension(nifti_image *nim, int size,int code, int rem);
308 static void update_nifti_image_for_brick_list(nifti_image * nim , int nbricks);
309 static int nifti_add_exten_to_list(nifti1_extension * new_ext,
310 nifti1_extension ** list, int new_length);
311 static int nifti_fill_extension(nifti1_extension * ext, const char * data,
312 int len, int ecode);
313
314 /* NBL routines */
315 static int nifti_load_NBL_bricks(nifti_image * nim , int * slist, int * sindex, nifti_brick_list * NBL, znzFile fp );
316 static int nifti_alloc_NBL_mem( nifti_image * nim, int nbricks,
317 nifti_brick_list * nbl);
318 static int nifti_copynsort(int nbricks, const int *blist, int **slist,
319 int **sindex);
320
321 /* for nifti_read_collapsed_image: */
322 static int rci_read_data(nifti_image *nim, int *pivots, int *prods, int nprods,
323 const int dims[], char *data, znzFile fp, int base_offset);
324 static int rci_alloc_mem(void ** data, int prods[8], int nprods, int nbyper );
325 static int make_pivot_list(nifti_image * nim, const int dims[], int pivots[],
326 int prods[], int * nprods );
327
328 /* misc */
329 static int need_nhdr_swap (short dim0, int hdrsize);
330 static int print_hex_vals (const char * data, int nbytes, FILE * fp);
331 static int unescape_string (char *str); /* string utility functions */
332 static char *escapize_string (const char *str);
333
334 /* internal I/O routines */
335 static znzFile nifti_image_load_prep( nifti_image *nim );
336 static int has_ascii_header(znzFile fp);
337 /*---------------------------------------------------------------------------*/
338
339
340 /* for calling from some main program */
341
342 /*----------------------------------------------------------------------*/
343 /*! display the nifti library module history (via stdout)
344 *//*--------------------------------------------------------------------*/
nifti_disp_lib_hist(void)345 void nifti_disp_lib_hist( void )
346 {
347 int c, len = sizeof(gni_history)/sizeof(char *);
348 for( c = 0; c < len; c++ )
349 fputs(gni_history[c], stdout);
350 }
351
352 /*----------------------------------------------------------------------*/
353 /*! display the nifti library version (via stdout)
354 *//*--------------------------------------------------------------------*/
nifti_disp_lib_version(void)355 void nifti_disp_lib_version( void )
356 {
357 printf("%s, compiled %s\n", gni_version, __DATE__);
358 }
359
360
361 /*----------------------------------------------------------------------*/
362 /*! nifti_image_read_bricks - read nifti data as array of bricks
363 *
364 * 13 Dec 2004 [rickr]
365 *
366 * \param hname - filename of dataset to read (must be valid)
367 * \param nbricks - number of sub-bricks to read
368 * (if blist is valid, nbricks must be > 0)
369 * \param blist - list of sub-bricks to read
370 * (can be NULL; if NULL, read complete dataset)
371 * \param NBL - pointer to empty nifti_brick_list struct
372 * (must be a valid pointer)
373 *
374 * \return
375 * <br> nim - same as nifti_image_read, but nim->data will be NULL
376 * <br> NBL - filled with data
377 *
378 * By default, this function will read the nifti dataset and break the data
379 * into a list of nt*nu*nv*nw sub-bricks, each having size nx*ny*nz elements.
380 * That is to say, instead of reading the entire dataset as a single array,
381 * break it up into sub-bricks, each of size nx*ny*nz elements.
382 *
383 * If 'blist' is valid, it is taken to be a list of sub-bricks, of length
384 * 'nbricks'. The data will still be separated into sub-bricks of size
385 * nx*ny*nz elements, but now 'nbricks' sub-bricks will be returned, of the
386 * caller's choosing via 'blist'.
387 *
388 * E.g. consider a dataset with 12 sub-bricks (numbered 0..11), and the
389 * following code:
390 *
391 * <pre>
392 * { nifti_brick_list NB_orig, NB_select;
393 * nifti_image * nim_orig, * nim_select;
394 * int blist[5] = { 7, 0, 5, 5, 9 };
395 *
396 * nim_orig = nifti_image_read_bricks("myfile.nii", 0, NULL, &NB_orig);
397 * nim_select = nifti_image_read_bricks("myfile.nii", 5, blist, &NB_select);
398 * }
399 * </pre>
400 *
401 * Here, nim_orig gets the entire dataset, where NB_orig.nbricks = 11. But
402 * nim_select has NB_select.nbricks = 5.
403 *
404 * Note that the first case is not quite the same as just calling the
405 * nifti_image_read function, as here the data is separated into sub-bricks.
406 *
407 * Note that valid blist elements are in [0..nt*nu*nv*nw-1],
408 * or written [ 0 .. (dim[4]*dim[5]*dim[6]*dim[7] - 1) ].
409 *
410 * Note that, as is the case with all of the reading functions, the
411 * data will be allocated, read in, and properly byte-swapped, if
412 * necessary.
413 *
414 * \sa nifti_image_load_bricks, nifti_free_NBL, valid_nifti_brick_list,
415 nifti_image_read
416 *//*----------------------------------------------------------------------*/
nifti_image_read_bricks(const char * hname,int nbricks,const int * blist,nifti_brick_list * NBL)417 nifti_image *nifti_image_read_bricks(const char * hname, int nbricks,
418 const int * blist, nifti_brick_list * NBL)
419 {
420 nifti_image * nim;
421
422 if( !hname || !NBL ){
423 fprintf(stderr,"** nifti_image_read_bricks: bad params (%p,%p)\n",
424 hname, (void *)NBL);
425 return NULL;
426 }
427
428 if( blist && nbricks <= 0 ){
429 fprintf(stderr,"** nifti_image_read_bricks: bad nbricks, %d\n", nbricks);
430 return NULL;
431 }
432
433 nim = nifti_image_read(hname, 0); /* read header, but not data */
434
435 if( !nim ) return NULL; /* errors were already printed */
436
437 /* if we fail, free image and return */
438 if( nifti_image_load_bricks(nim, nbricks, blist, NBL) <= 0 ){
439 nifti_image_free(nim);
440 return NULL;
441 }
442
443 if( blist ) update_nifti_image_for_brick_list(nim, nbricks);
444
445 return nim;
446 }
447
448
449 /*----------------------------------------------------------------------
450 * update_nifti_image_for_brick_list - update nifti_image
451 *
452 * When loading a specific brick list, the distinction between
453 * nt, nu, nv and nw is lost. So put everything in t, and set
454 * dim[0] = 4.
455 *----------------------------------------------------------------------*/
update_nifti_image_for_brick_list(nifti_image * nim,int nbricks)456 static void update_nifti_image_for_brick_list( nifti_image * nim , int nbricks )
457 {
458 int ndim;
459
460 if( g_opts.debug > 2 ){
461 fprintf(stderr,"+d updating image dimensions for %d bricks in list\n",
462 nbricks);
463 fprintf(stderr," ndim = %d\n",nim->ndim);
464 fprintf(stderr," nx,ny,nz,nt,nu,nv,nw: (%d,%d,%d,%d,%d,%d,%d)\n",
465 nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
466 }
467
468 nim->nt = nbricks;
469 nim->nu = nim->nv = nim->nw = 1;
470 nim->dim[4] = nbricks;
471 nim->dim[5] = nim->dim[6] = nim->dim[7] = 1;
472
473 /* compute nvox */
474 /* do not rely on dimensions above dim[0] 16 Nov 2005 [rickr] */
475 for( nim->nvox = 1, ndim = 1; ndim <= nim->dim[0]; ndim++ )
476 nim->nvox *= nim->dim[ndim];
477
478 /* update the dimensions to 4 or lower */
479 for( ndim = 4; (ndim > 1) && (nim->dim[ndim] <= 1); ndim-- )
480 ;
481
482 if( g_opts.debug > 2 ){
483 fprintf(stderr,"+d ndim = %d -> %d\n",nim->ndim, ndim);
484 fprintf(stderr," --> (%d,%d,%d,%d,%d,%d,%d)\n",
485 nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
486 }
487
488 nim->dim[0] = nim->ndim = ndim;
489 }
490
491
492 /*----------------------------------------------------------------------*/
493 /*! nifti_update_dims_from_array - update nx, ny, ... from nim->dim[]
494
495 Fix all the dimension information, based on a new nim->dim[].
496
497 Note: we assume that dim[0] will not increase.
498
499 Check for updates to pixdim[], dx,..., nx,..., nvox, ndim, dim[0].
500 *//*--------------------------------------------------------------------*/
nifti_update_dims_from_array(nifti_image * nim)501 int nifti_update_dims_from_array( nifti_image * nim )
502 {
503 int c, ndim;
504
505 if( !nim ){
506 fprintf(stderr,"** update_dims: missing nim\n");
507 return 1;
508 }
509
510 if( g_opts.debug > 2 ){
511 fprintf(stderr,"+d updating image dimensions given nim->dim:");
512 for( c = 0; c < 8; c++ ) fprintf(stderr," %d", nim->dim[c]);
513 fputc('\n',stderr);
514 }
515
516 /* verify dim[0] first */
517 if(nim->dim[0] < 1 || nim->dim[0] > 7){
518 fprintf(stderr,"** invalid dim[0], dim[] = ");
519 for( c = 0; c < 8; c++ ) fprintf(stderr," %d", nim->dim[c]);
520 fputc('\n',stderr);
521 return 1;
522 }
523
524 /* set nx, ny ..., dx, dy, ..., one by one */
525
526 /* less than 1, set to 1, else copy */
527 if(nim->dim[1] < 1) nim->nx = nim->dim[1] = 1;
528 else nim->nx = nim->dim[1];
529 nim->dx = nim->pixdim[1];
530
531 /* if undefined, or less than 1, set to 1 */
532 if(nim->dim[0] < 2 || (nim->dim[0] >= 2 && nim->dim[2] < 1))
533 nim->ny = nim->dim[2] = 1;
534 else
535 nim->ny = nim->dim[2];
536 /* copy delta values, in any case */
537 nim->dy = nim->pixdim[2];
538
539 if(nim->dim[0] < 3 || (nim->dim[0] >= 3 && nim->dim[3] < 1))
540 nim->nz = nim->dim[3] = 1;
541 else /* just copy vals from arrays */
542 nim->nz = nim->dim[3];
543 nim->dz = nim->pixdim[3];
544
545 if(nim->dim[0] < 4 || (nim->dim[0] >= 4 && nim->dim[4] < 1))
546 nim->nt = nim->dim[4] = 1;
547 else /* just copy vals from arrays */
548 nim->nt = nim->dim[4];
549 nim->dt = nim->pixdim[4];
550
551 if(nim->dim[0] < 5 || (nim->dim[0] >= 5 && nim->dim[5] < 1))
552 nim->nu = nim->dim[5] = 1;
553 else /* just copy vals from arrays */
554 nim->nu = nim->dim[5];
555 nim->du = nim->pixdim[5];
556
557 if(nim->dim[0] < 6 || (nim->dim[0] >= 6 && nim->dim[6] < 1))
558 nim->nv = nim->dim[6] = 1;
559 else /* just copy vals from arrays */
560 nim->nv = nim->dim[6];
561 nim->dv = nim->pixdim[6];
562
563 if(nim->dim[0] < 7 || (nim->dim[0] >= 7 && nim->dim[7] < 1))
564 nim->nw = nim->dim[7] = 1;
565 else /* just copy vals from arrays */
566 nim->nw = nim->dim[7];
567 nim->dw = nim->pixdim[7];
568
569 for( c = 1, nim->nvox = 1; c <= nim->dim[0]; c++ )
570 nim->nvox *= nim->dim[c];
571
572 /* compute ndim, assuming it can be no larger than the old one */
573 for( ndim = nim->dim[0]; (ndim > 1) && (nim->dim[ndim] <= 1); ndim-- )
574 ;
575
576 if( g_opts.debug > 2 ){
577 fprintf(stderr,"+d ndim = %d -> %d\n",nim->ndim, ndim);
578 fprintf(stderr," --> (%d,%d,%d,%d,%d,%d,%d)\n",
579 nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
580 }
581
582 nim->dim[0] = nim->ndim = ndim;
583
584 return 0;
585 }
586
587
588 /*----------------------------------------------------------------------*/
589 /*! Load the image data from disk into an already-prepared image struct.
590 *
591 * \param nim - initialized nifti_image, without data
592 * \param nbricks - the length of blist (must be 0 if blist is NULL)
593 * \param blist - an array of xyz volume indices to read (can be NULL)
594 * \param NBL - pointer to struct where resulting data will be stored
595 *
596 * If blist is NULL, read all sub-bricks.
597 *
598 * \return the number of loaded bricks (NBL->nbricks),
599 * 0 on failure, < 0 on error
600 *
601 * NOTE: it is likely that another function will copy the data pointers
602 * out of NBL, in which case the only pointer the calling function
603 * will want to free is NBL->bricks (not each NBL->bricks[i]).
604 *//*--------------------------------------------------------------------*/
nifti_image_load_bricks(nifti_image * nim,int nbricks,const int * blist,nifti_brick_list * NBL)605 int nifti_image_load_bricks( nifti_image * nim , int nbricks,
606 const int * blist, nifti_brick_list * NBL )
607 {
608 int * slist = NULL, * sindex = NULL, rv;
609 znzFile fp;
610
611 /* we can have blist == NULL */
612 if( !nim || !NBL ){
613 fprintf(stderr,"** nifti_image_load_bricks, bad params (%p,%p)\n",
614 (void *)nim, (void *)NBL);
615 return -1;
616 }
617
618 if( blist && nbricks <= 0 ){
619 if( g_opts.debug > 1 )
620 fprintf(stderr,"-d load_bricks: received blist with nbricks = %d,"
621 "ignoring blist\n", nbricks);
622 blist = NULL; /* pretend nothing was passed */
623 }
624
625 if( blist && ! valid_nifti_brick_list( nim, nbricks, blist, g_opts.debug>0 ) )
626 return -1;
627
628 /* for efficiency, let's read the file in order */
629 if( blist && nifti_copynsort( nbricks, blist, &slist, &sindex ) != 0 )
630 return -1;
631
632 /* open the file and position the FILE pointer */
633 fp = nifti_image_load_prep( nim );
634 if( !fp ){
635 if( g_opts.debug > 0 )
636 fprintf(stderr,"** nifti_image_load_bricks, failed load_prep\n");
637 if( blist ){ free(slist); free(sindex); }
638 return -1;
639 }
640
641 /* this will flag to allocate defaults */
642 if( !blist ) nbricks = 0;
643 if( nifti_alloc_NBL_mem( nim, nbricks, NBL ) != 0 ){
644 if( blist ){ free(slist); free(sindex); }
645 znzclose(fp);
646 return -1;
647 }
648
649 rv = nifti_load_NBL_bricks(nim, slist, sindex, NBL, fp);
650
651 if( rv != 0 ){
652 nifti_free_NBL( NBL ); /* failure! */
653 NBL->nbricks = 0; /* repetative, but clear */
654 }
655
656 if( slist ){ free(slist); free(sindex); }
657
658 znzclose(fp);
659
660 return NBL->nbricks;
661 }
662
663
664 /*----------------------------------------------------------------------*/
665 /*! nifti_free_NBL - free all pointers and clear structure
666 *
667 * note: this does not presume to free the structure pointer
668 *//*--------------------------------------------------------------------*/
nifti_free_NBL(nifti_brick_list * NBL)669 void nifti_free_NBL( nifti_brick_list * NBL )
670 {
671 int c;
672
673 if( NBL->bricks ){
674 for( c = 0; c < NBL->nbricks; c++ )
675 if( NBL->bricks[c] ) free(NBL->bricks[c]);
676 free(NBL->bricks);
677 NBL->bricks = NULL;
678 }
679
680 NBL->nbricks = NBL->bsize = 0;
681 }
682
683
684 /*----------------------------------------------------------------------
685 * nifti_load_NBL_bricks - read the file data into the NBL struct
686 *
687 * return 0 on success, -1 on failure
688 *----------------------------------------------------------------------*/
nifti_load_NBL_bricks(nifti_image * nim,int * slist,int * sindex,nifti_brick_list * NBL,znzFile fp)689 static int nifti_load_NBL_bricks( nifti_image * nim , int * slist, int * sindex,
690 nifti_brick_list * NBL, znzFile fp )
691 {
692 int c, rv;
693 int oposn, fposn; /* orig and current file positions */
694 int prev, isrc, idest; /* previous and current sub-brick, and new index */
695
696 oposn = znztell(fp); /* store current file position */
697 fposn = oposn;
698 if( fposn < 0 ){
699 fprintf(stderr,"** load bricks: ztell failed??\n");
700 return -1;
701 }
702
703 /* first, handle the default case, no passed blist */
704 if( !slist ){
705 for( c = 0; c < NBL->nbricks; c++ ) {
706 rv = nifti_read_buffer(fp, NBL->bricks[c], NBL->bsize, nim);
707 if( rv != NBL->bsize ){
708 fprintf(stderr,"** load bricks: cannot read brick %d from '%s'\n",
709 c, nim->iname ? nim->iname : nim->fname);
710 return -1;
711 }
712 }
713 if( g_opts.debug > 1 )
714 fprintf(stderr,"+d read %d default bricks from file %s\n",
715 NBL->nbricks, nim->iname ? nim->iname : nim->fname );
716 return 0;
717 }
718
719 if( !sindex ){
720 fprintf(stderr,"** load_NBL_bricks: missing index list\n");
721 return -1;
722 }
723
724 prev = -1; /* use prev for previous sub-brick */
725 for( c = 0; c < NBL->nbricks; c++ ){
726 isrc = slist[c]; /* this is original brick index (c is new one) */
727 idest = sindex[c]; /* this is the destination index for this data */
728
729 /* if this sub-brick is not the previous, we must read from disk */
730 if( isrc != prev ){
731
732 /* if we are not looking at the correct sub-brick, scan forward */
733 if( fposn != (oposn + isrc*NBL->bsize) ){
734 fposn = oposn + isrc*NBL->bsize;
735 if( znzseek(fp, fposn, SEEK_SET) < 0 ){
736 fprintf(stderr,"** failed to locate brick %d in file '%s'\n",
737 isrc, nim->iname ? nim->iname : nim->fname);
738 return -1;
739 }
740 }
741
742 /* only 10,000 lines later and we're actually reading something! */
743 rv = nifti_read_buffer(fp, NBL->bricks[idest], NBL->bsize, nim);
744 if( rv != NBL->bsize ){
745 fprintf(stderr,"** failed to read brick %d from file '%s'\n",
746 isrc, nim->iname ? nim->iname : nim->fname);
747 return -1;
748 }
749 fposn += NBL->bsize;
750 } else {
751 /* we have already read this sub-brick, just copy the previous one */
752 /* note that this works because they are sorted */
753 memcpy(NBL->bricks[idest], NBL->bricks[sindex[c-1]], NBL->bsize);
754 }
755
756 prev = isrc; /* in any case, note the now previous sub-brick */
757 }
758
759 return 0;
760 }
761
762
763 /*----------------------------------------------------------------------
764 * nifti_alloc_NBL_mem - allocate memory for bricks
765 *
766 * return 0 on success, -1 on failure
767 *----------------------------------------------------------------------*/
nifti_alloc_NBL_mem(nifti_image * nim,int nbricks,nifti_brick_list * nbl)768 static int nifti_alloc_NBL_mem(nifti_image * nim, int nbricks,
769 nifti_brick_list * nbl)
770 {
771 int c;
772
773 /* if nbricks is not specified, use the default */
774 if( nbricks > 0 ) nbl->nbricks = nbricks;
775 else { /* I missed this one with the 1.17 change 02 Mar 2006 [rickr] */
776 nbl->nbricks = 1;
777 for( c = 4; c <= nim->ndim; c++ )
778 nbl->nbricks *= nim->dim[c];
779 }
780
781 nbl->bsize = nim->nx * nim->ny * nim->nz * nim->nbyper; /* bytes */
782 nbl->bricks = (void **)malloc(nbl->nbricks * sizeof(void *));
783
784 if( ! nbl->bricks ){
785 fprintf(stderr,"** NANM: failed to alloc %d void ptrs\n",nbricks);
786 return -1;
787 }
788
789 for( c = 0; c < nbl->nbricks; c++ ){
790 nbl->bricks[c] = (void *)malloc(nbl->bsize);
791 if( ! nbl->bricks[c] ){
792 fprintf(stderr,"** NANM: failed to alloc %d bytes for brick %d\n",
793 nbl->bsize, c);
794 /* so free and clear everything before returning */
795 while( c > 0 ){
796 c--;
797 free(nbl->bricks[c]);
798 }
799 free(nbl->bricks);
800 nbl->bricks = NULL;
801 nbl->nbricks = nbl->bsize = 0;
802 return -1;
803 }
804 }
805
806 if( g_opts.debug > 2 )
807 fprintf(stderr,"+d NANM: alloc'd %d bricks of %d bytes for NBL\n",
808 nbl->nbricks, nbl->bsize);
809
810 return 0;
811 }
812
813
814 /*----------------------------------------------------------------------
815 * nifti_copynsort - copy int list, and sort with indices
816 *
817 * 1. duplicate the incoming list
818 * 2. create an sindex list, and init with 0..nbricks-1
819 * 3. do a slow insertion sort on the small slist, along with sindex list
820 * 4. check results, just to be positive
821 *
822 * So slist is sorted, and sindex hold original positions.
823 *
824 * return 0 on success, -1 on failure
825 *----------------------------------------------------------------------*/
nifti_copynsort(int nbricks,const int * blist,int ** slist,int ** sindex)826 static int nifti_copynsort(int nbricks, const int * blist, int ** slist,
827 int ** sindex)
828 {
829 int * stmp, * itmp; /* for ease of typing/reading */
830 int c1, c2, spos, tmp;
831
832 *slist = (int *)malloc(nbricks * sizeof(int));
833 *sindex = (int *)malloc(nbricks * sizeof(int));
834
835 if( !*slist || !*sindex ){
836 fprintf(stderr,"** NCS: failed to alloc %d ints for sorting\n",nbricks);
837 if(*slist) free(*slist); /* maybe one succeeded */
838 if(*sindex) free(*sindex);
839 return -1;
840 }
841
842 /* init the lists */
843 memcpy(*slist, blist, nbricks*sizeof(int));
844 for( c1 = 0; c1 < nbricks; c1++ ) (*sindex)[c1] = c1;
845
846 /* now actually sort slist */
847 stmp = *slist;
848 itmp = *sindex;
849 for( c1 = 0; c1 < nbricks-1; c1++ ) {
850 /* find smallest value, init to current */
851 spos = c1;
852 for( c2 = c1+1; c2 < nbricks; c2++ )
853 if( stmp[c2] < stmp[spos] ) spos = c2;
854 if( spos != c1 ) /* swap: fine, don't maintain sub-order, see if I care */
855 {
856 tmp = stmp[c1]; /* first swap the sorting values */
857 stmp[c1] = stmp[spos];
858 stmp[spos] = tmp;
859
860 tmp = itmp[c1]; /* then swap the index values */
861 itmp[c1] = itmp[spos];
862 itmp[spos] = tmp;
863 }
864 }
865
866 if( g_opts.debug > 2 ){
867 fprintf(stderr, "+d sorted indexing list:\n");
868 fprintf(stderr, " orig : ");
869 for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",blist[c1]);
870 fprintf(stderr,"\n new : ");
871 for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",stmp[c1]);
872 fprintf(stderr,"\n indices: ");
873 for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr," %d",itmp[c1]);
874 fputc('\n', stderr);
875 }
876
877 /* check the sort (why not? I've got time...) */
878 for( c1 = 0; c1 < nbricks-1; c1++ ){
879 if( (stmp[c1] > stmp[c1+1]) || (blist[itmp[c1]] != stmp[c1]) ){
880 fprintf(stderr,"** sorting screw-up, way to go, rick!\n");
881 free(stmp); free(itmp); *slist = NULL; *sindex = NULL;
882 return -1;
883 }
884 }
885
886 if( g_opts.debug > 2 ) fprintf(stderr,"-d sorting is okay\n");
887
888 return 0;
889 }
890
891
892 /*----------------------------------------------------------------------*/
893 /*! valid_nifti_brick_list - check sub-brick list for image
894 *
895 * This function verifies that nbricks and blist are appropriate
896 * for use with this nim, based on the dimensions.
897 *
898 * \param nim nifti_image to check against
899 * \param nbricks number of brick indices in blist
900 * \param blist list of brick indices to check in nim
901 * \param disp_error if this flag is set, report errors to user
902 *
903 * \return 1 if valid, 0 if not
904 *//*--------------------------------------------------------------------*/
valid_nifti_brick_list(nifti_image * nim,int nbricks,const int * blist,int disp_error)905 int valid_nifti_brick_list(nifti_image * nim , int nbricks,
906 const int * blist, int disp_error)
907 {
908 int c, nsubs;
909
910 if( !nim ){
911 if( disp_error || g_opts.debug > 0 )
912 fprintf(stderr,"** valid_nifti_brick_list: missing nifti image\n");
913 return 0;
914 }
915
916 if( nbricks <= 0 || !blist ){
917 if( disp_error || g_opts.debug > 1 )
918 fprintf(stderr,"** valid_nifti_brick_list: no brick list to check\n");
919 return 0;
920 }
921
922 /* nsubs sub-brick is nt*nu*nv*nw */
923 for( c = 4, nsubs = 1; c <= nim->dim[0]; c++ )
924 nsubs *= nim->dim[c];
925
926 if( nsubs <= 0 ){
927 fprintf(stderr,"** VNBL warning: bad dim list (%d,%d,%d,%d)\n",
928 nim->dim[4], nim->dim[5], nim->dim[6], nim->dim[7]);
929 return 0;
930 }
931
932 for( c = 0; c < nbricks; c++ )
933 if( (blist[c] < 0) || (blist[c] >= nsubs) ){
934 if( disp_error || g_opts.debug > 1 )
935 fprintf(stderr,
936 "-d ** bad sub-brick chooser %d (#%d), valid range is [0,%d]\n",
937 blist[c], c, nsubs-1);
938 return 0;
939 }
940
941 return 1; /* all is well */
942 }
943
944 #if 0
945 /* set any non-positive values to 1 */
946 static int int_force_positive( int * list, int nel )
947 {
948 int c;
949 if( !list || nel < 0 ){
950 if( g_opts.debug > 0 )
951 fprintf(stderr,"** int_force_positive: bad params (%p,%d)\n",
952 (void *)list,nel);
953 return -1;
954 }
955 for( c = 0; c < nel; c++ )
956 if( list[c] <= 0 ) list[c] = 1;
957 return 0;
958 }
959 #endif
960 /* end of new nifti_image_read_bricks() functionality */
961
962 /*----------------------------------------------------------------------*/
963 /*! display the orientation from the quaternian fields
964 *
965 * \param mesg if non-NULL, display this message first
966 * \param mat the matrix to convert to "nearest" orientation
967 *
968 * \return -1 if results cannot be determined, 0 if okay
969 *//*--------------------------------------------------------------------*/
nifti_disp_matrix_orient(const char * mesg,mat44 mat)970 int nifti_disp_matrix_orient( const char * mesg, mat44 mat )
971 {
972 int i, j, k;
973
974 if ( mesg ) fputs( mesg, stderr ); /* use stdout? */
975
976 nifti_mat44_to_orientation( mat, &i,&j,&k );
977 if ( i <= 0 || j <= 0 || k <= 0 ) return -1;
978
979 /* so we have good codes */
980 fprintf(stderr, " i orientation = '%s'\n"
981 " j orientation = '%s'\n"
982 " k orientation = '%s'\n",
983 nifti_orientation_string(i),
984 nifti_orientation_string(j),
985 nifti_orientation_string(k) );
986 return 0;
987 }
988
989
990 /*----------------------------------------------------------------------*/
991 /*! duplicate the given string (alloc length+1)
992 *
993 * \return allocated pointer (or NULL on failure)
994 *//*--------------------------------------------------------------------*/
nifti_strdup(const char * str)995 char *nifti_strdup(const char *str)
996 {
997 char *dup= (char *)malloc( strlen(str)+1 );
998 if (dup) strcpy(dup,str);
999 return dup;
1000 }
1001
1002
1003 /*---------------------------------------------------------------------------*/
1004 /*! Return a pointer to a string holding the name of a NIFTI datatype.
1005
1006 \param dt NIfTI-1 datatype
1007
1008 \return pointer to static string holding the datatype name
1009
1010 \warning Do not free() or modify this string!
1011 It points to static storage.
1012
1013 \sa NIFTI1_DATATYPES group in nifti1.h
1014 *//*-------------------------------------------------------------------------*/
nifti_datatype_string(int dt)1015 char *nifti_datatype_string( int dt )
1016 {
1017 switch( dt ){
1018 case DT_UNKNOWN: return "UNKNOWN" ;
1019 case DT_BINARY: return "BINARY" ;
1020 case DT_INT8: return "INT8" ;
1021 case DT_UINT8: return "UINT8" ;
1022 case DT_INT16: return "INT16" ;
1023 case DT_UINT16: return "UINT16" ;
1024 case DT_INT32: return "INT32" ;
1025 case DT_UINT32: return "UINT32" ;
1026 case DT_INT64: return "INT64" ;
1027 case DT_UINT64: return "UINT64" ;
1028 case DT_FLOAT32: return "FLOAT32" ;
1029 case DT_FLOAT64: return "FLOAT64" ;
1030 case DT_FLOAT128: return "FLOAT128" ;
1031 case DT_COMPLEX64: return "COMPLEX64" ;
1032 case DT_COMPLEX128: return "COMPLEX128" ;
1033 case DT_COMPLEX256: return "COMPLEX256" ;
1034 case DT_RGB24: return "RGB24" ;
1035 }
1036 return "**ILLEGAL**" ;
1037 }
1038
1039 /*----------------------------------------------------------------------*/
1040 /*! Determine if the datatype code dt is an integer type (1=YES, 0=NO).
1041
1042 \return whether the given NIfTI-1 datatype code is valid
1043
1044 \sa NIFTI1_DATATYPES group in nifti1.h
1045 *//*--------------------------------------------------------------------*/
nifti_is_inttype(int dt)1046 int nifti_is_inttype( int dt )
1047 {
1048 switch( dt ){
1049 case DT_UNKNOWN: return 0 ;
1050 case DT_BINARY: return 0 ;
1051 case DT_INT8: return 1 ;
1052 case DT_UINT8: return 1 ;
1053 case DT_INT16: return 1 ;
1054 case DT_UINT16: return 1 ;
1055 case DT_INT32: return 1 ;
1056 case DT_UINT32: return 1 ;
1057 case DT_INT64: return 1 ;
1058 case DT_UINT64: return 1 ;
1059 case DT_FLOAT32: return 0 ;
1060 case DT_FLOAT64: return 0 ;
1061 case DT_FLOAT128: return 0 ;
1062 case DT_COMPLEX64: return 0 ;
1063 case DT_COMPLEX128: return 0 ;
1064 case DT_COMPLEX256: return 0 ;
1065 case DT_RGB24: return 1 ;
1066 }
1067 return 0 ;
1068 }
1069
1070 /*---------------------------------------------------------------------------*/
1071 /*! Return a pointer to a string holding the name of a NIFTI units type.
1072
1073 \param uu NIfTI-1 unit code
1074
1075 \return pointer to static string for the given unit type
1076
1077 \warning Do not free() or modify this string!
1078 It points to static storage.
1079
1080 \sa NIFTI1_UNITS group in nifti1.h
1081 *//*-------------------------------------------------------------------------*/
nifti_units_string(int uu)1082 char *nifti_units_string( int uu )
1083 {
1084 switch( uu ){
1085 case NIFTI_UNITS_METER: return "m" ;
1086 case NIFTI_UNITS_MM: return "mm" ;
1087 case NIFTI_UNITS_MICRON: return "um" ;
1088 case NIFTI_UNITS_SEC: return "s" ;
1089 case NIFTI_UNITS_MSEC: return "ms" ;
1090 case NIFTI_UNITS_USEC: return "us" ;
1091 case NIFTI_UNITS_HZ: return "Hz" ;
1092 case NIFTI_UNITS_PPM: return "ppm" ;
1093 case NIFTI_UNITS_RADS: return "rad/s" ;
1094 }
1095 return "Unknown" ;
1096 }
1097
1098 /*---------------------------------------------------------------------------*/
1099 /*! Return a pointer to a string holding the name of a NIFTI transform type.
1100
1101 \param xx NIfTI-1 xform code
1102
1103 \return pointer to static string describing xform code
1104
1105 \warning Do not free() or modify this string!
1106 It points to static storage.
1107
1108 \sa NIFTI1_XFORM_CODES group in nifti1.h
1109 *//*-------------------------------------------------------------------------*/
nifti_xform_string(int xx)1110 char *nifti_xform_string( int xx )
1111 {
1112 switch( xx ){
1113 case NIFTI_XFORM_SCANNER_ANAT: return "Scanner Anat" ;
1114 case NIFTI_XFORM_ALIGNED_ANAT: return "Aligned Anat" ;
1115 case NIFTI_XFORM_TALAIRACH: return "Talairach" ;
1116 case NIFTI_XFORM_MNI_152: return "MNI_152" ;
1117 }
1118 return "Unknown" ;
1119 }
1120
1121 /*---------------------------------------------------------------------------*/
1122 /*! Return a pointer to a string holding the name of a NIFTI intent type.
1123
1124 \param ii NIfTI-1 intent code
1125
1126 \return pointer to static string describing code
1127
1128 \warning Do not free() or modify this string!
1129 It points to static storage.
1130
1131 \sa NIFTI1_INTENT_CODES group in nifti1.h
1132 *//*-------------------------------------------------------------------------*/
nifti_intent_string(int ii)1133 char *nifti_intent_string( int ii )
1134 {
1135 switch( ii ){
1136 case NIFTI_INTENT_CORREL: return "Correlation statistic" ;
1137 case NIFTI_INTENT_TTEST: return "T-statistic" ;
1138 case NIFTI_INTENT_FTEST: return "F-statistic" ;
1139 case NIFTI_INTENT_ZSCORE: return "Z-score" ;
1140 case NIFTI_INTENT_CHISQ: return "Chi-squared distribution" ;
1141 case NIFTI_INTENT_BETA: return "Beta distribution" ;
1142 case NIFTI_INTENT_BINOM: return "Binomial distribution" ;
1143 case NIFTI_INTENT_GAMMA: return "Gamma distribution" ;
1144 case NIFTI_INTENT_POISSON: return "Poisson distribution" ;
1145 case NIFTI_INTENT_NORMAL: return "Normal distribution" ;
1146 case NIFTI_INTENT_FTEST_NONC: return "F-statistic noncentral" ;
1147 case NIFTI_INTENT_CHISQ_NONC: return "Chi-squared noncentral" ;
1148 case NIFTI_INTENT_LOGISTIC: return "Logistic distribution" ;
1149 case NIFTI_INTENT_LAPLACE: return "Laplace distribution" ;
1150 case NIFTI_INTENT_UNIFORM: return "Uniform distribition" ;
1151 case NIFTI_INTENT_TTEST_NONC: return "T-statistic noncentral" ;
1152 case NIFTI_INTENT_WEIBULL: return "Weibull distribution" ;
1153 case NIFTI_INTENT_CHI: return "Chi distribution" ;
1154 case NIFTI_INTENT_INVGAUSS: return "Inverse Gaussian distribution" ;
1155 case NIFTI_INTENT_EXTVAL: return "Extreme Value distribution" ;
1156 case NIFTI_INTENT_PVAL: return "P-value" ;
1157
1158 case NIFTI_INTENT_LOGPVAL: return "Log P-value" ;
1159 case NIFTI_INTENT_LOG10PVAL: return "Log10 P-value" ;
1160
1161 case NIFTI_INTENT_ESTIMATE: return "Estimate" ;
1162 case NIFTI_INTENT_LABEL: return "Label index" ;
1163 case NIFTI_INTENT_NEURONAME: return "NeuroNames index" ;
1164 case NIFTI_INTENT_GENMATRIX: return "General matrix" ;
1165 case NIFTI_INTENT_SYMMATRIX: return "Symmetric matrix" ;
1166 case NIFTI_INTENT_DISPVECT: return "Displacement vector" ;
1167 case NIFTI_INTENT_VECTOR: return "Vector" ;
1168 case NIFTI_INTENT_POINTSET: return "Pointset" ;
1169 case NIFTI_INTENT_TRIANGLE: return "Triangle" ;
1170 case NIFTI_INTENT_QUATERNION: return "Quaternion" ;
1171
1172 case NIFTI_INTENT_DIMLESS: return "Dimensionless number" ;
1173 }
1174 return "Unknown" ;
1175 }
1176
1177 /*---------------------------------------------------------------------------*/
1178 /*! Return a pointer to a string holding the name of a NIFTI slice_code.
1179
1180 \param ss NIfTI-1 slice order code
1181
1182 \return pointer to static string describing code
1183
1184 \warning Do not free() or modify this string!
1185 It points to static storage.
1186
1187 \sa NIFTI1_SLICE_ORDER group in nifti1.h
1188 *//*-------------------------------------------------------------------------*/
nifti_slice_string(int ss)1189 char *nifti_slice_string( int ss )
1190 {
1191 switch( ss ){
1192 case NIFTI_SLICE_SEQ_INC: return "sequential_increasing" ;
1193 case NIFTI_SLICE_SEQ_DEC: return "sequential_decreasing" ;
1194 case NIFTI_SLICE_ALT_INC: return "alternating_increasing" ;
1195 case NIFTI_SLICE_ALT_DEC: return "alternating_decreasing" ;
1196 case NIFTI_SLICE_ALT_INC2: return "alternating_increasing_2" ;
1197 case NIFTI_SLICE_ALT_DEC2: return "alternating_decreasing_2" ;
1198 }
1199 return "Unknown" ;
1200 }
1201
1202 /*---------------------------------------------------------------------------*/
1203 /*! Return a pointer to a string holding the name of a NIFTI orientation.
1204
1205 \param ii orientation code
1206
1207 \return pointer to static string holding the orientation information
1208
1209 \warning Do not free() or modify the return string!
1210 It points to static storage.
1211
1212 \sa NIFTI_L2R in nifti1_io.h
1213 *//*-------------------------------------------------------------------------*/
nifti_orientation_string(int ii)1214 char *nifti_orientation_string( int ii )
1215 {
1216 switch( ii ){
1217 case NIFTI_L2R: return "Left-to-Right" ;
1218 case NIFTI_R2L: return "Right-to-Left" ;
1219 case NIFTI_P2A: return "Posterior-to-Anterior" ;
1220 case NIFTI_A2P: return "Anterior-to-Posterior" ;
1221 case NIFTI_I2S: return "Inferior-to-Superior" ;
1222 case NIFTI_S2I: return "Superior-to-Inferior" ;
1223 }
1224 return "Unknown" ;
1225 }
1226
1227 /*--------------------------------------------------------------------------*/
1228 /*! Given a datatype code, set number of bytes per voxel and the swapsize.
1229
1230 \param datatype nifti1 datatype code
1231 \param nbyper pointer to return value: number of bytes per voxel
1232 \param swapsize pointer to return value: size of swap blocks
1233
1234 \return appropriate values at nbyper and swapsize
1235
1236 The swapsize is set to 0 if this datatype doesn't ever need swapping.
1237
1238 \sa NIFTI1_DATATYPES in nifti1.h
1239 *//*------------------------------------------------------------------------*/
nifti_datatype_sizes(int datatype,int * nbyper,int * swapsize)1240 void nifti_datatype_sizes( int datatype , int *nbyper, int *swapsize )
1241 {
1242 int nb=0, ss=0 ;
1243 switch( datatype ){
1244 case DT_INT8:
1245 case DT_UINT8: nb = 1 ; ss = 0 ; break ;
1246
1247 case DT_INT16:
1248 case DT_UINT16: nb = 2 ; ss = 2 ; break ;
1249
1250 case DT_RGB24: nb = 3 ; ss = 0 ; break ;
1251
1252 case DT_INT32:
1253 case DT_UINT32:
1254 case DT_FLOAT32: nb = 4 ; ss = 4 ; break ;
1255
1256 case DT_COMPLEX64: nb = 8 ; ss = 4 ; break ;
1257
1258 case DT_FLOAT64:
1259 case DT_INT64:
1260 case DT_UINT64: nb = 8 ; ss = 8 ; break ;
1261
1262 case DT_FLOAT128: nb = 16 ; ss = 16 ; break ;
1263
1264 case DT_COMPLEX128: nb = 16 ; ss = 8 ; break ;
1265
1266 case DT_COMPLEX256: nb = 32 ; ss = 16 ; break ;
1267 }
1268
1269 ASSIF(nbyper,nb) ; ASSIF(swapsize,ss) ; return ;
1270 }
1271
1272 /*---------------------------------------------------------------------------*/
1273 /*! Given the quaternion parameters (etc.), compute a transformation matrix.
1274
1275 See comments in nifti1.h for details.
1276 - qb,qc,qd = quaternion parameters
1277 - qx,qy,qz = offset parameters
1278 - dx,dy,dz = grid stepsizes (non-negative inputs are set to 1.0)
1279 - qfac = sign of dz step (< 0 is negative; >= 0 is positive)
1280
1281 <pre>
1282 If qx=qy=qz=0, dx=dy=dz=1, then the output is a rotation matrix.
1283 For qfac >= 0, the rotation is proper.
1284 For qfac < 0, the rotation is improper.
1285 </pre>
1286
1287 \see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h
1288 \see nifti_mat44_to_quatern, nifti_make_orthog_mat44,
1289 nifti_mat44_to_orientation
1290
1291 *//*-------------------------------------------------------------------------*/
nifti_quatern_to_mat44(float qb,float qc,float qd,float qx,float qy,float qz,float dx,float dy,float dz,float qfac)1292 mat44 nifti_quatern_to_mat44( float qb, float qc, float qd,
1293 float qx, float qy, float qz,
1294 float dx, float dy, float dz, float qfac )
1295 {
1296 mat44 R ;
1297 double a,b=qb,c=qc,d=qd , xd,yd,zd ;
1298
1299 /* last row is always [ 0 0 0 1 ] */
1300
1301 R.m[3][0]=R.m[3][1]=R.m[3][2] = 0.0 ; R.m[3][3]= 1.0 ;
1302
1303 /* compute a parameter from b,c,d */
1304
1305 a = 1.0l - (b*b + c*c + d*d) ;
1306 if( a < 1.e-7l ){ /* special case */
1307 a = 1.0l / sqrt(b*b+c*c+d*d) ;
1308 b *= a ; c *= a ; d *= a ; /* normalize (b,c,d) vector */
1309 a = 0.0l ; /* a = 0 ==> 180 degree rotation */
1310 } else{
1311 a = sqrt(a) ; /* angle = 2*arccos(a) */
1312 }
1313
1314 /* load rotation matrix, including scaling factors for voxel sizes */
1315
1316 xd = (dx > 0.0) ? dx : 1.0l ; /* make sure are positive */
1317 yd = (dy > 0.0) ? dy : 1.0l ;
1318 zd = (dz > 0.0) ? dz : 1.0l ;
1319
1320 if( qfac < 0.0 ) zd = -zd ; /* left handedness? */
1321
1322 R.m[0][0] = (a*a+b*b-c*c-d*d) * xd ;
1323 R.m[0][1] = 2.0l * (b*c-a*d ) * yd ;
1324 R.m[0][2] = 2.0l * (b*d+a*c ) * zd ;
1325 R.m[1][0] = 2.0l * (b*c+a*d ) * xd ;
1326 R.m[1][1] = (a*a+c*c-b*b-d*d) * yd ;
1327 R.m[1][2] = 2.0l * (c*d-a*b ) * zd ;
1328 R.m[2][0] = 2.0l * (b*d-a*c ) * xd ;
1329 R.m[2][1] = 2.0l * (c*d+a*b ) * yd ;
1330 R.m[2][2] = (a*a+d*d-c*c-b*b) * zd ;
1331
1332 /* load offsets */
1333
1334 R.m[0][3] = qx ; R.m[1][3] = qy ; R.m[2][3] = qz ;
1335
1336 return R ;
1337 }
1338
1339 /*---------------------------------------------------------------------------*/
1340 /*! Given the 3x4 upper corner of the matrix R, compute the quaternion
1341 parameters that fit it.
1342
1343 - Any NULL pointer on input won't get assigned (e.g., if you don't want
1344 dx,dy,dz, just pass NULL in for those pointers).
1345 - If the 3 input matrix columns are NOT orthogonal, they will be
1346 orthogonalized prior to calculating the parameters, using
1347 the polar decomposition to find the orthogonal matrix closest
1348 to the column-normalized input matrix.
1349 - However, if the 3 input matrix columns are NOT orthogonal, then
1350 the matrix produced by nifti_quatern_to_mat44 WILL have orthogonal
1351 columns, so it won't be the same as the matrix input here.
1352 This "feature" is because the NIFTI 'qform' transform is
1353 deliberately not fully general -- it is intended to model a volume
1354 with perpendicular axes.
1355 - If the 3 input matrix columns are not even linearly independent,
1356 you'll just have to take your luck, won't you?
1357
1358 \see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h
1359
1360 \see nifti_quatern_to_mat44, nifti_make_orthog_mat44,
1361 nifti_mat44_to_orientation
1362 *//*-------------------------------------------------------------------------*/
nifti_mat44_to_quatern(mat44 R,float * qb,float * qc,float * qd,float * qx,float * qy,float * qz,float * dx,float * dy,float * dz,float * qfac)1363 void nifti_mat44_to_quatern( mat44 R ,
1364 float *qb, float *qc, float *qd,
1365 float *qx, float *qy, float *qz,
1366 float *dx, float *dy, float *dz, float *qfac )
1367 {
1368 double r11,r12,r13 , r21,r22,r23 , r31,r32,r33 ;
1369 double xd,yd,zd , a,b,c,d ;
1370 mat33 P,Q ;
1371
1372 /* offset outputs are read write out of input matrix */
1373
1374 ASSIF(qx,R.m[0][3]) ; ASSIF(qy,R.m[1][3]) ; ASSIF(qz,R.m[2][3]) ;
1375
1376 /* load 3x3 matrix into local variables */
1377
1378 r11 = R.m[0][0] ; r12 = R.m[0][1] ; r13 = R.m[0][2] ;
1379 r21 = R.m[1][0] ; r22 = R.m[1][1] ; r23 = R.m[1][2] ;
1380 r31 = R.m[2][0] ; r32 = R.m[2][1] ; r33 = R.m[2][2] ;
1381
1382 /* compute lengths of each column; these determine grid spacings */
1383
1384 xd = sqrt( r11*r11 + r21*r21 + r31*r31 ) ;
1385 yd = sqrt( r12*r12 + r22*r22 + r32*r32 ) ;
1386 zd = sqrt( r13*r13 + r23*r23 + r33*r33 ) ;
1387
1388 /* if a column length is zero, patch the trouble */
1389
1390 if( xd == 0.0l ){ r11 = 1.0l ; r21 = r31 = 0.0l ; xd = 1.0l ; }
1391 if( yd == 0.0l ){ r22 = 1.0l ; r12 = r32 = 0.0l ; yd = 1.0l ; }
1392 if( zd == 0.0l ){ r33 = 1.0l ; r13 = r23 = 0.0l ; zd = 1.0l ; }
1393
1394 /* assign the output lengths */
1395
1396 ASSIF(dx,xd) ; ASSIF(dy,yd) ; ASSIF(dz,zd) ;
1397
1398 /* normalize the columns */
1399
1400 r11 /= xd ; r21 /= xd ; r31 /= xd ;
1401 r12 /= yd ; r22 /= yd ; r32 /= yd ;
1402 r13 /= zd ; r23 /= zd ; r33 /= zd ;
1403
1404 /* At this point, the matrix has normal columns, but we have to allow
1405 for the fact that the hideous user may not have given us a matrix
1406 with orthogonal columns.
1407
1408 So, now find the orthogonal matrix closest to the current matrix.
1409
1410 One reason for using the polar decomposition to get this
1411 orthogonal matrix, rather than just directly orthogonalizing
1412 the columns, is so that inputting the inverse matrix to R
1413 will result in the inverse orthogonal matrix at this point.
1414 If we just orthogonalized the columns, this wouldn't necessarily hold. */
1415
1416 Q.m[0][0] = r11 ; Q.m[0][1] = r12 ; Q.m[0][2] = r13 ; /* load Q */
1417 Q.m[1][0] = r21 ; Q.m[1][1] = r22 ; Q.m[1][2] = r23 ;
1418 Q.m[2][0] = r31 ; Q.m[2][1] = r32 ; Q.m[2][2] = r33 ;
1419
1420 P = nifti_mat33_polar(Q) ; /* P is orthog matrix closest to Q */
1421
1422 r11 = P.m[0][0] ; r12 = P.m[0][1] ; r13 = P.m[0][2] ; /* unload */
1423 r21 = P.m[1][0] ; r22 = P.m[1][1] ; r23 = P.m[1][2] ;
1424 r31 = P.m[2][0] ; r32 = P.m[2][1] ; r33 = P.m[2][2] ;
1425
1426 /* [ r11 r12 r13 ] */
1427 /* at this point, the matrix [ r21 r22 r23 ] is orthogonal */
1428 /* [ r31 r32 r33 ] */
1429
1430 /* compute the determinant to determine if it is proper */
1431
1432 zd = r11*r22*r33-r11*r32*r23-r21*r12*r33
1433 +r21*r32*r13+r31*r12*r23-r31*r22*r13 ; /* should be -1 or 1 */
1434
1435 if( zd > 0 ){ /* proper */
1436 ASSIF(qfac,1.0) ;
1437 } else { /* improper ==> flip 3rd column */
1438 ASSIF(qfac,-1.0) ;
1439 r13 = -r13 ; r23 = -r23 ; r33 = -r33 ;
1440 }
1441
1442 /* now, compute quaternion parameters */
1443
1444 a = r11 + r22 + r33 + 1.0l ;
1445
1446 if( a > 0.5l ){ /* simplest case */
1447 a = 0.5l * sqrt(a) ;
1448 b = 0.25l * (r32-r23) / a ;
1449 c = 0.25l * (r13-r31) / a ;
1450 d = 0.25l * (r21-r12) / a ;
1451 } else { /* trickier case */
1452 xd = 1.0 + r11 - (r22+r33) ; /* 4*b*b */
1453 yd = 1.0 + r22 - (r11+r33) ; /* 4*c*c */
1454 zd = 1.0 + r33 - (r11+r22) ; /* 4*d*d */
1455 if( xd > 1.0 ){
1456 b = 0.5l * sqrt(xd) ;
1457 c = 0.25l* (r12+r21) / b ;
1458 d = 0.25l* (r13+r31) / b ;
1459 a = 0.25l* (r32-r23) / b ;
1460 } else if( yd > 1.0 ){
1461 c = 0.5l * sqrt(yd) ;
1462 b = 0.25l* (r12+r21) / c ;
1463 d = 0.25l* (r23+r32) / c ;
1464 a = 0.25l* (r13-r31) / c ;
1465 } else {
1466 d = 0.5l * sqrt(zd) ;
1467 b = 0.25l* (r13+r31) / d ;
1468 c = 0.25l* (r23+r32) / d ;
1469 a = 0.25l* (r21-r12) / d ;
1470 }
1471 if( a < 0.0l ){ b=-b ; c=-c ; d=-d; a=-a; }
1472 }
1473
1474 ASSIF(qb,b) ; ASSIF(qc,c) ; ASSIF(qd,d) ;
1475 return ;
1476 }
1477
1478 /*---------------------------------------------------------------------------*/
1479 /*! Compute the inverse of a bordered 4x4 matrix.
1480
1481 <pre>
1482 - Some numerical code fragments were generated by Maple 8.
1483 - If a singular matrix is input, the output matrix will be all zero.
1484 - You can check for this by examining the [3][3] element, which will
1485 be 1.0 for the normal case and 0.0 for the bad case.
1486
1487 The input matrix should have the form:
1488 [ r11 r12 r13 v1 ]
1489 [ r21 r22 r23 v2 ]
1490 [ r31 r32 r33 v3 ]
1491 [ 0 0 0 1 ]
1492 </pre>
1493 *//*-------------------------------------------------------------------------*/
nifti_mat44_inverse(mat44 R)1494 mat44 nifti_mat44_inverse( mat44 R )
1495 {
1496 double r11,r12,r13,r21,r22,r23,r31,r32,r33,v1,v2,v3 , deti ;
1497 mat44 Q ;
1498 /* INPUT MATRIX IS: */
1499 r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 v1 ] */
1500 r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 v2 ] */
1501 r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 v3 ] */
1502 v1 = R.m[0][3]; v2 = R.m[1][3]; v3 = R.m[2][3]; /* [ 0 0 0 1 ] */
1503
1504 deti = r11*r22*r33-r11*r32*r23-r21*r12*r33
1505 +r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
1506
1507 if( deti != 0.0l ) deti = 1.0l / deti ;
1508
1509 Q.m[0][0] = deti*( r22*r33-r32*r23) ;
1510 Q.m[0][1] = deti*(-r12*r33+r32*r13) ;
1511 Q.m[0][2] = deti*( r12*r23-r22*r13) ;
1512 Q.m[0][3] = deti*(-r12*r23*v3+r12*v2*r33+r22*r13*v3
1513 -r22*v1*r33-r32*r13*v2+r32*v1*r23) ;
1514
1515 Q.m[1][0] = deti*(-r21*r33+r31*r23) ;
1516 Q.m[1][1] = deti*( r11*r33-r31*r13) ;
1517 Q.m[1][2] = deti*(-r11*r23+r21*r13) ;
1518 Q.m[1][3] = deti*( r11*r23*v3-r11*v2*r33-r21*r13*v3
1519 +r21*v1*r33+r31*r13*v2-r31*v1*r23) ;
1520
1521 Q.m[2][0] = deti*( r21*r32-r31*r22) ;
1522 Q.m[2][1] = deti*(-r11*r32+r31*r12) ;
1523 Q.m[2][2] = deti*( r11*r22-r21*r12) ;
1524 Q.m[2][3] = deti*(-r11*r22*v3+r11*r32*v2+r21*r12*v3
1525 -r21*r32*v1-r31*r12*v2+r31*r22*v1) ;
1526
1527 Q.m[3][0] = Q.m[3][1] = Q.m[3][2] = 0.0l ;
1528 Q.m[3][3] = (deti == 0.0l) ? 0.0l : 1.0l ; /* failure flag if deti == 0 */
1529
1530 return Q ;
1531 }
1532
1533 /*---------------------------------------------------------------------------*/
1534 /*! Input 9 floats and make an orthgonal mat44 out of them.
1535
1536 Each row is normalized, then nifti_mat33_polar() is used to orthogonalize
1537 them. If row #3 (r31,r32,r33) is input as zero, then it will be taken to
1538 be the cross product of rows #1 and #2.
1539
1540 This function can be used to create a rotation matrix for transforming
1541 an oblique volume to anatomical coordinates. For this application:
1542 - row #1 (r11,r12,r13) is the direction vector along the image i-axis
1543 - row #2 (r21,r22,r23) is the direction vector along the image j-axis
1544 - row #3 (r31,r32,r33) is the direction vector along the slice direction
1545 (if available; otherwise enter it as 0's)
1546
1547 The first 2 rows can be taken from the DICOM attribute (0020,0037)
1548 "Image Orientation (Patient)".
1549
1550 After forming the rotation matrix, the complete affine transformation from
1551 (i,j,k) grid indexes to (x,y,z) spatial coordinates can be computed by
1552 multiplying each column by the appropriate grid spacing:
1553 - column #1 (R.m[0][0],R.m[1][0],R.m[2][0]) by delta-x
1554 - column #2 (R.m[0][1],R.m[1][1],R.m[2][1]) by delta-y
1555 - column #3 (R.m[0][2],R.m[1][2],R.m[2][2]) by delta-z
1556
1557 and by then placing the center (x,y,z) coordinates of voxel (0,0,0) into
1558 the column #4 (R.m[0][3],R.m[1][3],R.m[2][3]).
1559
1560 \sa nifti_quatern_to_mat44, nifti_mat44_to_quatern,
1561 nifti_mat44_to_orientation
1562 *//*-------------------------------------------------------------------------*/
nifti_make_orthog_mat44(float r11,float r12,float r13,float r21,float r22,float r23,float r31,float r32,float r33)1563 mat44 nifti_make_orthog_mat44( float r11, float r12, float r13 ,
1564 float r21, float r22, float r23 ,
1565 float r31, float r32, float r33 )
1566 {
1567 mat44 R ;
1568 mat33 Q , P ;
1569 double val ;
1570
1571 R.m[3][0] = R.m[3][1] = R.m[3][2] = 0.0l ; R.m[3][3] = 1.0l ;
1572
1573 Q.m[0][0] = r11 ; Q.m[0][1] = r12 ; Q.m[0][2] = r13 ; /* load Q */
1574 Q.m[1][0] = r21 ; Q.m[1][1] = r22 ; Q.m[1][2] = r23 ;
1575 Q.m[2][0] = r31 ; Q.m[2][1] = r32 ; Q.m[2][2] = r33 ;
1576
1577 /* normalize row 1 */
1578
1579 val = Q.m[0][0]*Q.m[0][0] + Q.m[0][1]*Q.m[0][1] + Q.m[0][2]*Q.m[0][2] ;
1580 if( val > 0.0l ){
1581 val = 1.0l / sqrt(val) ;
1582 Q.m[0][0] *= val ; Q.m[0][1] *= val ; Q.m[0][2] *= val ;
1583 } else {
1584 Q.m[0][0] = 1.0l ; Q.m[0][1] = 0.0l ; Q.m[0][2] = 0.0l ;
1585 }
1586
1587 /* normalize row 2 */
1588
1589 val = Q.m[1][0]*Q.m[1][0] + Q.m[1][1]*Q.m[1][1] + Q.m[1][2]*Q.m[1][2] ;
1590 if( val > 0.0l ){
1591 val = 1.0l / sqrt(val) ;
1592 Q.m[1][0] *= val ; Q.m[1][1] *= val ; Q.m[1][2] *= val ;
1593 } else {
1594 Q.m[1][0] = 0.0l ; Q.m[1][1] = 1.0l ; Q.m[1][2] = 0.0l ;
1595 }
1596
1597 /* normalize row 3 */
1598
1599 val = Q.m[2][0]*Q.m[2][0] + Q.m[2][1]*Q.m[2][1] + Q.m[2][2]*Q.m[2][2] ;
1600 if( val > 0.0l ){
1601 val = 1.0l / sqrt(val) ;
1602 Q.m[2][0] *= val ; Q.m[2][1] *= val ; Q.m[2][2] *= val ;
1603 } else {
1604 Q.m[2][0] = Q.m[0][1]*Q.m[1][2] - Q.m[0][2]*Q.m[1][1] ; /* cross */
1605 Q.m[2][1] = Q.m[0][2]*Q.m[1][0] - Q.m[0][0]*Q.m[1][2] ; /* product */
1606 Q.m[2][2] = Q.m[0][0]*Q.m[1][1] - Q.m[0][1]*Q.m[1][0] ;
1607 }
1608
1609 P = nifti_mat33_polar(Q) ; /* P is orthog matrix closest to Q */
1610
1611 R.m[0][0] = P.m[0][0] ; R.m[0][1] = P.m[0][1] ; R.m[0][2] = P.m[0][2] ;
1612 R.m[1][0] = P.m[1][0] ; R.m[1][1] = P.m[1][1] ; R.m[1][2] = P.m[1][2] ;
1613 R.m[2][0] = P.m[2][0] ; R.m[2][1] = P.m[2][1] ; R.m[2][2] = P.m[2][2] ;
1614
1615 R.m[0][3] = R.m[1][3] = R.m[2][3] = 0.0 ; return R ;
1616 }
1617
1618 /*----------------------------------------------------------------------*/
1619 /*! compute the inverse of a 3x3 matrix
1620 *//*--------------------------------------------------------------------*/
nifti_mat33_inverse(mat33 R)1621 mat33 nifti_mat33_inverse( mat33 R ) /* inverse of 3x3 matrix */
1622 {
1623 double r11,r12,r13,r21,r22,r23,r31,r32,r33 , deti ;
1624 mat33 Q ;
1625 /* INPUT MATRIX: */
1626 r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 ] */
1627 r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 ] */
1628 r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 ] */
1629
1630 deti = r11*r22*r33-r11*r32*r23-r21*r12*r33
1631 +r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
1632
1633 if( deti != 0.0l ) deti = 1.0l / deti ;
1634
1635 Q.m[0][0] = deti*( r22*r33-r32*r23) ;
1636 Q.m[0][1] = deti*(-r12*r33+r32*r13) ;
1637 Q.m[0][2] = deti*( r12*r23-r22*r13) ;
1638
1639 Q.m[1][0] = deti*(-r21*r33+r31*r23) ;
1640 Q.m[1][1] = deti*( r11*r33-r31*r13) ;
1641 Q.m[1][2] = deti*(-r11*r23+r21*r13) ;
1642
1643 Q.m[2][0] = deti*( r21*r32-r31*r22) ;
1644 Q.m[2][1] = deti*(-r11*r32+r31*r12) ;
1645 Q.m[2][2] = deti*( r11*r22-r21*r12) ;
1646
1647 return Q ;
1648 }
1649
1650 /*----------------------------------------------------------------------*/
1651 /*! compute the determinant of a 3x3 matrix
1652 *//*--------------------------------------------------------------------*/
nifti_mat33_determ(mat33 R)1653 float nifti_mat33_determ( mat33 R ) /* determinant of 3x3 matrix */
1654 {
1655 double r11,r12,r13,r21,r22,r23,r31,r32,r33 ;
1656 /* INPUT MATRIX: */
1657 r11 = R.m[0][0]; r12 = R.m[0][1]; r13 = R.m[0][2]; /* [ r11 r12 r13 ] */
1658 r21 = R.m[1][0]; r22 = R.m[1][1]; r23 = R.m[1][2]; /* [ r21 r22 r23 ] */
1659 r31 = R.m[2][0]; r32 = R.m[2][1]; r33 = R.m[2][2]; /* [ r31 r32 r33 ] */
1660
1661 return r11*r22*r33-r11*r32*r23-r21*r12*r33
1662 +r21*r32*r13+r31*r12*r23-r31*r22*r13 ;
1663 }
1664
1665 /*----------------------------------------------------------------------*/
1666 /*! compute the max row norm of a 3x3 matrix
1667 *//*--------------------------------------------------------------------*/
nifti_mat33_rownorm(mat33 A)1668 float nifti_mat33_rownorm( mat33 A ) /* max row norm of 3x3 matrix */
1669 {
1670 float r1,r2,r3 ;
1671
1672 r1 = fabs(A.m[0][0])+fabs(A.m[0][1])+fabs(A.m[0][2]) ;
1673 r2 = fabs(A.m[1][0])+fabs(A.m[1][1])+fabs(A.m[1][2]) ;
1674 r3 = fabs(A.m[2][0])+fabs(A.m[2][1])+fabs(A.m[2][2]) ;
1675 if( r1 < r2 ) r1 = r2 ;
1676 if( r1 < r3 ) r1 = r3 ;
1677 return r1 ;
1678 }
1679
1680 /*----------------------------------------------------------------------*/
1681 /*! compute the max column norm of a 3x3 matrix
1682 *//*--------------------------------------------------------------------*/
nifti_mat33_colnorm(mat33 A)1683 float nifti_mat33_colnorm( mat33 A ) /* max column norm of 3x3 matrix */
1684 {
1685 float r1,r2,r3 ;
1686
1687 r1 = fabs(A.m[0][0])+fabs(A.m[1][0])+fabs(A.m[2][0]) ;
1688 r2 = fabs(A.m[0][1])+fabs(A.m[1][1])+fabs(A.m[2][1]) ;
1689 r3 = fabs(A.m[0][2])+fabs(A.m[1][2])+fabs(A.m[2][2]) ;
1690 if( r1 < r2 ) r1 = r2 ;
1691 if( r1 < r3 ) r1 = r3 ;
1692 return r1 ;
1693 }
1694
1695 /*----------------------------------------------------------------------*/
1696 /*! multiply 2 3x3 matrices
1697 *//*--------------------------------------------------------------------*/
nifti_mat33_mul(mat33 A,mat33 B)1698 mat33 nifti_mat33_mul( mat33 A , mat33 B ) /* multiply 2 3x3 matrices */
1699 {
1700 mat33 C ; int i,j ;
1701 for( i=0 ; i < 3 ; i++ )
1702 for( j=0 ; j < 3 ; j++ )
1703 C.m[i][j] = A.m[i][0] * B.m[0][j]
1704 + A.m[i][1] * B.m[1][j]
1705 + A.m[i][2] * B.m[2][j] ;
1706 return C ;
1707 }
1708
1709 /*---------------------------------------------------------------------------*/
1710 /*! polar decomposition of a 3x3 matrix
1711
1712 This finds the closest orthogonal matrix to input A
1713 (in both the Frobenius and L2 norms).
1714
1715 Algorithm is that from NJ Higham, SIAM J Sci Stat Comput, 7:1160-1174.
1716 *//*-------------------------------------------------------------------------*/
nifti_mat33_polar(mat33 A)1717 mat33 nifti_mat33_polar( mat33 A )
1718 {
1719 mat33 X , Y , Z ;
1720 float alp,bet,gam,gmi , dif=1.0 ;
1721 int k=0 ;
1722
1723 X = A ;
1724
1725 /* force matrix to be nonsingular */
1726
1727 gam = nifti_mat33_determ(X) ;
1728 while( gam == 0.0 ){ /* perturb matrix */
1729 gam = 0.00001 * ( 0.001 + nifti_mat33_rownorm(X) ) ;
1730 X.m[0][0] += gam ; X.m[1][1] += gam ; X.m[2][2] += gam ;
1731 gam = nifti_mat33_determ(X) ;
1732 }
1733
1734 while(1){
1735 Y = nifti_mat33_inverse(X) ;
1736 if( dif > 0.3 ){ /* far from convergence */
1737 alp = sqrt( nifti_mat33_rownorm(X) * nifti_mat33_colnorm(X) ) ;
1738 bet = sqrt( nifti_mat33_rownorm(Y) * nifti_mat33_colnorm(Y) ) ;
1739 gam = sqrt( bet / alp ) ;
1740 gmi = 1.0 / gam ;
1741 } else {
1742 gam = gmi = 1.0 ; /* close to convergence */
1743 }
1744 Z.m[0][0] = 0.5 * ( gam*X.m[0][0] + gmi*Y.m[0][0] ) ;
1745 Z.m[0][1] = 0.5 * ( gam*X.m[0][1] + gmi*Y.m[1][0] ) ;
1746 Z.m[0][2] = 0.5 * ( gam*X.m[0][2] + gmi*Y.m[2][0] ) ;
1747 Z.m[1][0] = 0.5 * ( gam*X.m[1][0] + gmi*Y.m[0][1] ) ;
1748 Z.m[1][1] = 0.5 * ( gam*X.m[1][1] + gmi*Y.m[1][1] ) ;
1749 Z.m[1][2] = 0.5 * ( gam*X.m[1][2] + gmi*Y.m[2][1] ) ;
1750 Z.m[2][0] = 0.5 * ( gam*X.m[2][0] + gmi*Y.m[0][2] ) ;
1751 Z.m[2][1] = 0.5 * ( gam*X.m[2][1] + gmi*Y.m[1][2] ) ;
1752 Z.m[2][2] = 0.5 * ( gam*X.m[2][2] + gmi*Y.m[2][2] ) ;
1753
1754 dif = fabs(Z.m[0][0]-X.m[0][0])+fabs(Z.m[0][1]-X.m[0][1])
1755 +fabs(Z.m[0][2]-X.m[0][2])+fabs(Z.m[1][0]-X.m[1][0])
1756 +fabs(Z.m[1][1]-X.m[1][1])+fabs(Z.m[1][2]-X.m[1][2])
1757 +fabs(Z.m[2][0]-X.m[2][0])+fabs(Z.m[2][1]-X.m[2][1])
1758 +fabs(Z.m[2][2]-X.m[2][2]) ;
1759
1760 k = k+1 ;
1761 if( k > 100 || dif < 3.e-6 ) break ; /* convergence or exhaustion */
1762 X = Z ;
1763 }
1764
1765 return Z ;
1766 }
1767
1768 /*---------------------------------------------------------------------------*/
1769 /*! compute the (closest) orientation from a 4x4 ijk->xyz tranformation matrix
1770
1771 <pre>
1772 Input: 4x4 matrix that transforms (i,j,k) indexes to (x,y,z) coordinates,
1773 where +x=Right, +y=Anterior, +z=Superior.
1774 (Only the upper-left 3x3 corner of R is used herein.)
1775 Output: 3 orientation codes that correspond to the closest "standard"
1776 anatomical orientation of the (i,j,k) axes.
1777 Method: Find which permutation of (x,y,z) has the smallest angle to the
1778 (i,j,k) axes directions, which are the columns of the R matrix.
1779 Errors: The codes returned will be zero.
1780
1781 For example, an axial volume might get return values of
1782 *icod = NIFTI_R2L (i axis is mostly Right to Left)
1783 *jcod = NIFTI_P2A (j axis is mostly Posterior to Anterior)
1784 *kcod = NIFTI_I2S (k axis is mostly Inferior to Superior)
1785 </pre>
1786
1787 \see "QUATERNION REPRESENTATION OF ROTATION MATRIX" in nifti1.h
1788
1789 \see nifti_quatern_to_mat44, nifti_mat44_to_quatern,
1790 nifti_make_orthog_mat44
1791 *//*-------------------------------------------------------------------------*/
nifti_mat44_to_orientation(mat44 R,int * icod,int * jcod,int * kcod)1792 void nifti_mat44_to_orientation( mat44 R , int *icod, int *jcod, int *kcod )
1793 {
1794 float xi,xj,xk , yi,yj,yk , zi,zj,zk , val,detQ,detP ;
1795 mat33 P , Q , M ;
1796 int i,j,k=0,p,q,r , ibest,jbest,kbest,pbest,qbest,rbest ;
1797 float vbest ;
1798
1799 if( icod == NULL || jcod == NULL || kcod == NULL ) return ; /* bad */
1800
1801 *icod = *jcod = *kcod = 0 ; /* error returns, if sh*t happens */
1802
1803 /* load column vectors for each (i,j,k) direction from matrix */
1804
1805 /*-- i axis --*/ /*-- j axis --*/ /*-- k axis --*/
1806
1807 xi = R.m[0][0] ; xj = R.m[0][1] ; xk = R.m[0][2] ;
1808 yi = R.m[1][0] ; yj = R.m[1][1] ; yk = R.m[1][2] ;
1809 zi = R.m[2][0] ; zj = R.m[2][1] ; zk = R.m[2][2] ;
1810
1811 /* normalize column vectors to get unit vectors along each ijk-axis */
1812
1813 /* normalize i axis */
1814
1815 val = sqrt( xi*xi + yi*yi + zi*zi ) ;
1816 if( val == 0.0 ) return ; /* stupid input */
1817 xi /= val ; yi /= val ; zi /= val ;
1818
1819 /* normalize j axis */
1820
1821 val = sqrt( xj*xj + yj*yj + zj*zj ) ;
1822 if( val == 0.0 ) return ; /* stupid input */
1823 xj /= val ; yj /= val ; zj /= val ;
1824
1825 /* orthogonalize j axis to i axis, if needed */
1826
1827 val = xi*xj + yi*yj + zi*zj ; /* dot product between i and j */
1828 if( fabs(val) > 1.e-4 ){
1829 xj -= val*xi ; yj -= val*yi ; zj -= val*zi ;
1830 val = sqrt( xj*xj + yj*yj + zj*zj ) ; /* must renormalize */
1831 if( val == 0.0 ) return ; /* j was parallel to i? */
1832 xj /= val ; yj /= val ; zj /= val ;
1833 }
1834
1835 /* normalize k axis; if it is zero, make it the cross product i x j */
1836
1837 val = sqrt( xk*xk + yk*yk + zk*zk ) ;
1838 if( val == 0.0 ){ xk = yi*zj-zi*yj; yk = zi*xj-zj*xi ; zk=xi*yj-yi*xj ; }
1839 else { xk /= val ; yk /= val ; zk /= val ; }
1840
1841 /* orthogonalize k to i */
1842
1843 val = xi*xk + yi*yk + zi*zk ; /* dot product between i and k */
1844 if( fabs(val) > 1.e-4 ){
1845 xk -= val*xi ; yk -= val*yi ; zk -= val*zi ;
1846 val = sqrt( xk*xk + yk*yk + zk*zk ) ;
1847 if( val == 0.0 ) return ; /* bad */
1848 xk /= val ; yk /= val ; zk /= val ;
1849 }
1850
1851 /* orthogonalize k to j */
1852
1853 val = xj*xk + yj*yk + zj*zk ; /* dot product between j and k */
1854 if( fabs(val) > 1.e-4 ){
1855 xk -= val*xj ; yk -= val*yj ; zk -= val*zj ;
1856 val = sqrt( xk*xk + yk*yk + zk*zk ) ;
1857 if( val == 0.0 ) return ; /* bad */
1858 xk /= val ; yk /= val ; zk /= val ;
1859 }
1860
1861 Q.m[0][0] = xi ; Q.m[0][1] = xj ; Q.m[0][2] = xk ;
1862 Q.m[1][0] = yi ; Q.m[1][1] = yj ; Q.m[1][2] = yk ;
1863 Q.m[2][0] = zi ; Q.m[2][1] = zj ; Q.m[2][2] = zk ;
1864
1865 /* at this point, Q is the rotation matrix from the (i,j,k) to (x,y,z) axes */
1866
1867 detQ = nifti_mat33_determ( Q ) ;
1868 if( detQ == 0.0 ) return ; /* shouldn't happen unless user is a DUFIS */
1869
1870 /* Build and test all possible +1/-1 coordinate permutation matrices P;
1871 then find the P such that the rotation matrix M=PQ is closest to the
1872 identity, in the sense of M having the smallest total rotation angle. */
1873
1874 /* Despite the formidable looking 6 nested loops, there are
1875 only 3*3*3*2*2*2 = 216 passes, which will run very quickly. */
1876
1877 vbest = -666.0 ; ibest=pbest=qbest=rbest=1 ; jbest=2 ; kbest=3 ;
1878 for( i=1 ; i <= 3 ; i++ ){ /* i = column number to use for row #1 */
1879 for( j=1 ; j <= 3 ; j++ ){ /* j = column number to use for row #2 */
1880 if( i == j ) continue ;
1881 for( k=1 ; k <= 3 ; k++ ){ /* k = column number to use for row #3 */
1882 if( i == k || j == k ) continue ;
1883 P.m[0][0] = P.m[0][1] = P.m[0][2] =
1884 P.m[1][0] = P.m[1][1] = P.m[1][2] =
1885 P.m[2][0] = P.m[2][1] = P.m[2][2] = 0.0 ;
1886 for( p=-1 ; p <= 1 ; p+=2 ){ /* p,q,r are -1 or +1 */
1887 for( q=-1 ; q <= 1 ; q+=2 ){ /* and go into rows #1,2,3 */
1888 for( r=-1 ; r <= 1 ; r+=2 ){
1889 P.m[0][i-1] = p ; P.m[1][j-1] = q ; P.m[2][k-1] = r ;
1890 detP = nifti_mat33_determ(P) ; /* sign of permutation */
1891 if( detP * detQ <= 0.0 ) continue ; /* doesn't match sign of Q */
1892 M = nifti_mat33_mul(P,Q) ;
1893
1894 /* angle of M rotation = 2.0*acos(0.5*sqrt(1.0+trace(M))) */
1895 /* we want largest trace(M) == smallest angle == M nearest to I */
1896
1897 val = M.m[0][0] + M.m[1][1] + M.m[2][2] ; /* trace */
1898 if( val > vbest ){
1899 vbest = val ;
1900 ibest = i ; jbest = j ; kbest = k ;
1901 pbest = p ; qbest = q ; rbest = r ;
1902 }
1903 }}}}}}
1904
1905 /* At this point ibest is 1 or 2 or 3; pbest is -1 or +1; etc.
1906
1907 The matrix P that corresponds is the best permutation approximation
1908 to Q-inverse; that is, P (approximately) takes (x,y,z) coordinates
1909 to the (i,j,k) axes.
1910
1911 For example, the first row of P (which contains pbest in column ibest)
1912 determines the way the i axis points relative to the anatomical
1913 (x,y,z) axes. If ibest is 2, then the i axis is along the y axis,
1914 which is direction P2A (if pbest > 0) or A2P (if pbest < 0).
1915
1916 So, using ibest and pbest, we can assign the output code for
1917 the i axis. Mutatis mutandis for the j and k axes, of course. */
1918
1919 switch( ibest*pbest ){
1920 case 1: i = NIFTI_L2R ; break ;
1921 case -1: i = NIFTI_R2L ; break ;
1922 case 2: i = NIFTI_P2A ; break ;
1923 case -2: i = NIFTI_A2P ; break ;
1924 case 3: i = NIFTI_I2S ; break ;
1925 case -3: i = NIFTI_S2I ; break ;
1926 }
1927
1928 switch( jbest*qbest ){
1929 case 1: j = NIFTI_L2R ; break ;
1930 case -1: j = NIFTI_R2L ; break ;
1931 case 2: j = NIFTI_P2A ; break ;
1932 case -2: j = NIFTI_A2P ; break ;
1933 case 3: j = NIFTI_I2S ; break ;
1934 case -3: j = NIFTI_S2I ; break ;
1935 }
1936
1937 switch( kbest*rbest ){
1938 case 1: k = NIFTI_L2R ; break ;
1939 case -1: k = NIFTI_R2L ; break ;
1940 case 2: k = NIFTI_P2A ; break ;
1941 case -2: k = NIFTI_A2P ; break ;
1942 case 3: k = NIFTI_I2S ; break ;
1943 case -3: k = NIFTI_S2I ; break ;
1944 }
1945
1946 *icod = i ; *jcod = j ; *kcod = k ; return ;
1947 }
1948
1949 /*---------------------------------------------------------------------------*/
1950 /* Routines to swap byte arrays in various ways:
1951 - 2 at a time: ab -> ba [short]
1952 - 4 at a time: abcd -> dcba [int, float]
1953 - 8 at a time: abcdDCBA -> ABCDdcba [long long, double]
1954 - 16 at a time: abcdefghHGFEDCBA -> ABCDEFGHhgfedcba [long double]
1955 -----------------------------------------------------------------------------*/
1956
1957 typedef struct { unsigned char a,b ; } twobytes ;
1958
1959 /*----------------------------------------------------------------------*/
1960 /*! swap each byte pair from the given list of n pairs
1961 *//*--------------------------------------------------------------------*/
nifti_swap_2bytes(int n,void * ar)1962 void nifti_swap_2bytes( int n , void *ar ) /* 2 bytes at a time */
1963 {
1964 register int ii ;
1965 register twobytes *tb = (twobytes *)ar ;
1966 register unsigned char tt ;
1967
1968 for( ii=0 ; ii < n ; ii++ ){
1969 tt = tb[ii].a ; tb[ii].a = tb[ii].b ; tb[ii].b = tt ;
1970 }
1971 return ;
1972 }
1973
1974 /*---------------------------------------------------------------------------*/
1975
1976 typedef struct { unsigned char a,b,c,d ; } fourbytes ;
1977
1978 /*----------------------------------------------------------------------*/
1979 /*! swap 4 bytes at a time from the given list of n sets of 4 bytes
1980 *//*--------------------------------------------------------------------*/
nifti_swap_4bytes(int n,void * ar)1981 void nifti_swap_4bytes( int n , void *ar ) /* 4 bytes at a time */
1982 {
1983 register int ii ;
1984 register fourbytes *tb = (fourbytes *)ar ;
1985 register unsigned char tt ;
1986
1987 for( ii=0 ; ii < n ; ii++ ){
1988 tt = tb[ii].a ; tb[ii].a = tb[ii].d ; tb[ii].d = tt ;
1989 tt = tb[ii].b ; tb[ii].b = tb[ii].c ; tb[ii].c = tt ;
1990 }
1991 return ;
1992 }
1993
1994 /*---------------------------------------------------------------------------*/
1995
1996 typedef struct { unsigned char a,b,c,d , D,C,B,A ; } eightbytes ;
1997
1998 /*----------------------------------------------------------------------*/
1999 /*! swap 8 bytes at a time from the given list of n sets of 8 bytes
2000 *//*--------------------------------------------------------------------*/
nifti_swap_8bytes(int n,void * ar)2001 void nifti_swap_8bytes( int n , void *ar ) /* 8 bytes at a time */
2002 {
2003 register int ii ;
2004 register eightbytes *tb = (eightbytes *)ar ;
2005 register unsigned char tt ;
2006
2007 for( ii=0 ; ii < n ; ii++ ){
2008 tt = tb[ii].a ; tb[ii].a = tb[ii].A ; tb[ii].A = tt ;
2009 tt = tb[ii].b ; tb[ii].b = tb[ii].B ; tb[ii].B = tt ;
2010 tt = tb[ii].c ; tb[ii].c = tb[ii].C ; tb[ii].C = tt ;
2011 tt = tb[ii].d ; tb[ii].d = tb[ii].D ; tb[ii].D = tt ;
2012 }
2013 return ;
2014 }
2015
2016 /*---------------------------------------------------------------------------*/
2017
2018 typedef struct { unsigned char a,b,c,d,e,f,g,h ,
2019 H,G,F,E,D,C,B,A ; } sixteenbytes ;
2020
2021 /*----------------------------------------------------------------------*/
2022 /*! swap 16 bytes at a time from the given list of n sets of 16 bytes
2023 *//*--------------------------------------------------------------------*/
nifti_swap_16bytes(int n,void * ar)2024 void nifti_swap_16bytes( int n , void *ar ) /* 16 bytes at a time */
2025 {
2026 register int ii ;
2027 register sixteenbytes *tb = (sixteenbytes *)ar ;
2028 register unsigned char tt ;
2029
2030 for( ii=0 ; ii < n ; ii++ ){
2031 tt = tb[ii].a ; tb[ii].a = tb[ii].A ; tb[ii].A = tt ;
2032 tt = tb[ii].b ; tb[ii].b = tb[ii].B ; tb[ii].B = tt ;
2033 tt = tb[ii].c ; tb[ii].c = tb[ii].C ; tb[ii].C = tt ;
2034 tt = tb[ii].d ; tb[ii].d = tb[ii].D ; tb[ii].D = tt ;
2035
2036 tt = tb[ii].e ; tb[ii].e = tb[ii].E ; tb[ii].E = tt ;
2037 tt = tb[ii].f ; tb[ii].f = tb[ii].F ; tb[ii].F = tt ;
2038 tt = tb[ii].g ; tb[ii].g = tb[ii].G ; tb[ii].G = tt ;
2039 tt = tb[ii].h ; tb[ii].h = tb[ii].H ; tb[ii].H = tt ;
2040 }
2041 return ;
2042 }
2043
2044 /*---------------------------------------------------------------------------*/
2045
2046 /*----------------------------------------------------------------------*/
2047 /*! based on siz, call the appropriate nifti_swap_Nbytes() function
2048 *//*--------------------------------------------------------------------*/
nifti_swap_Nbytes(int n,int siz,void * ar)2049 void nifti_swap_Nbytes( int n , int siz , void *ar ) /* subsuming case */
2050 {
2051 switch( siz ){
2052 case 2: nifti_swap_2bytes ( n , ar ) ; break ;
2053 case 4: nifti_swap_4bytes ( n , ar ) ; break ;
2054 case 8: nifti_swap_8bytes ( n , ar ) ; break ;
2055 case 16: nifti_swap_16bytes( n , ar ) ; break ;
2056 }
2057 return ;
2058 }
2059
2060
2061 /*-------------------------------------------------------------------------*/
2062 /*! Byte swap NIFTI-1 file header in various places and ways.
2063
2064 If is_nifti is nonzero, will also swap the NIFTI-specific
2065 components of the header; otherwise, only the components
2066 common to NIFTI and ANALYZE will be swapped.
2067 *//*---------------------------------------------------------------------- */
swap_nifti_header(struct nifti_1_header * h,int is_nifti)2068 void swap_nifti_header( struct nifti_1_header *h , int is_nifti )
2069 {
2070
2071 #if 0 /* ANALYZE fields not used by this software */
2072 swap_4(h->sizeof_hdr) ;
2073 swap_4(h->extents) ;
2074 swap_2(h->session_error) ;
2075 swap_4(h->compressed) ;
2076 swap_4(h->glmax) ; swap_4(h->glmin) ;
2077 #endif
2078
2079 /* this stuff is always present, for ANALYZE and NIFTI */
2080
2081 swap_4(h->sizeof_hdr) ;
2082 nifti_swap_2bytes( 8 , h->dim ) ;
2083 nifti_swap_4bytes( 8 , h->pixdim ) ;
2084
2085 swap_2(h->datatype) ;
2086 swap_2(h->bitpix) ;
2087
2088 swap_4(h->vox_offset); swap_4(h->cal_max); swap_4(h->cal_min);
2089
2090 /* this stuff is NIFTI specific */
2091
2092 if( is_nifti ){
2093 swap_4(h->intent_p1); swap_4(h->intent_p2); swap_4(h->intent_p3);
2094 swap_2(h->intent_code);
2095
2096 swap_2(h->slice_start); swap_2(h->slice_end);
2097 swap_4(h->scl_slope); swap_4(h->scl_inter);
2098 swap_4(h->slice_duration); swap_4(h->toffset);
2099
2100 swap_2(h->qform_code); swap_2(h->sform_code);
2101 swap_4(h->quatern_b); swap_4(h->quatern_c); swap_4(h->quatern_d);
2102 swap_4(h->qoffset_x); swap_4(h->qoffset_y); swap_4(h->qoffset_z);
2103 nifti_swap_4bytes(4,h->srow_x);
2104 nifti_swap_4bytes(4,h->srow_y);
2105 nifti_swap_4bytes(4,h->srow_z);
2106 }
2107 return ;
2108 }
2109
2110
2111 #define USE_STAT
2112 #ifdef USE_STAT
2113 /*---------------------------------------------------------------------------*/
2114 /* Return the file length (0 if file not found or has no contents).
2115 This is a Unix-specific function, since it uses stat().
2116 -----------------------------------------------------------------------------*/
2117 #include <sys/types.h>
2118 #include <sys/stat.h>
2119
2120 /*---------------------------------------------------------------------------*/
2121 /*! return the size of a file, in bytes
2122
2123 \return size of file on success, -1 on error or no file
2124
2125 changed to return int, -1 means no file or error 20 Dec 2004 [rickr]
2126 *//*-------------------------------------------------------------------------*/
nifti_get_filesize(const char * pathname)2127 int nifti_get_filesize( const char *pathname )
2128 {
2129 struct stat buf ; int ii ;
2130
2131 if( pathname == NULL || *pathname == '\0' ) return -1 ;
2132 ii = stat( pathname , &buf ); if( ii != 0 ) return -1 ;
2133 return (unsigned int)buf.st_size ;
2134 }
2135
2136 #else /*---------- non-Unix version of the above, less efficient -----------*/
2137
nifti_get_filesize(const char * pathname)2138 int nifti_get_filesize( const char *pathname )
2139 {
2140 znzFile fp ; int len ;
2141
2142 if( pathname == NULL || *pathname == '\0' ) return -1 ;
2143 fp = znzopen(pathname,"rb",0); if( znz_isnull(fp) ) return -1 ;
2144 znzseek(fp,0L,SEEK_END) ; len = znztell(fp) ;
2145 znzclose(fp) ; return len ;
2146 }
2147
2148 #endif /* USE_STAT */
2149
2150
2151 /*----------------------------------------------------------------------*/
2152 /*! return the total volume size, in bytes
2153
2154 This is computed as nvox * nbyper.
2155 *//*--------------------------------------------------------------------*/
nifti_get_volsize(const nifti_image * nim)2156 size_t nifti_get_volsize(const nifti_image *nim)
2157 {
2158 return (size_t)(nim->nbyper) * (size_t)(nim->nvox) ; /* total bytes */
2159 }
2160
2161
2162 /*--------------------------------------------------------------------------*/
2163 /* Support functions for filenames in read and write
2164 - allows for gzipped files
2165 */
2166
2167
2168 /*----------------------------------------------------------------------*/
2169 /*! simple check for file existence
2170
2171 \return 1 on existence, 0 otherwise
2172 *//*--------------------------------------------------------------------*/
nifti_fileexists(const char * fname)2173 int nifti_fileexists(const char* fname)
2174 {
2175 znzFile fp;
2176 fp = znzopen( fname , "rb" , 1 ) ;
2177 if( !znz_isnull(fp) ) { znzclose(fp); return 1; }
2178 return 0; /* fp is NULL */
2179 }
2180
2181 /*----------------------------------------------------------------------*/
2182 /*! return whether the filename is valid
2183
2184 The name is considered valid if the file basename has length greater than
2185 zero, AND one of the valid nifti extensions is provided.
2186 fname input | return |
2187 ===============================
2188 "myimage" | 0 |
2189 "myimage.tif" | 0 |
2190 "myimage.tif.gz" | 0 |
2191 "myimage.nii" | 1 |
2192 ".nii" | 0 |
2193 ".myhiddenimage" | 0 |
2194 ".myhiddenimage.nii" | 1 |
2195 *//*--------------------------------------------------------------------*/
nifti_is_complete_filename(const char * fname)2196 int nifti_is_complete_filename(const char* fname)
2197 {
2198 char * ext;
2199
2200 /* check input file(s) for sanity */
2201 if( fname == NULL || *fname == '\0' ){
2202 if ( g_opts.debug > 1 )
2203 fprintf(stderr,"-- empty filename in nifti_validfilename()\n");
2204 return 0;
2205 }
2206
2207 ext = nifti_find_file_extension(fname);
2208 if ( ext == NULL ) { /*Invalid extension given */
2209 if ( g_opts.debug > 0 )
2210 fprintf(stderr,"-- no nifti valid extension for filename '%s'\n", fname);
2211 return 0;
2212 }
2213
2214 if ( ext && ext == fname ) { /* then no filename prefix */
2215 if ( g_opts.debug > 0 )
2216 fprintf(stderr,"-- no prefix for filename '%s'\n", fname);
2217 return 0;
2218 }
2219 return 1;
2220 }
2221
2222 /*----------------------------------------------------------------------*/
2223 /*! return whether the filename is valid
2224
2225 The name is considered valid if its length is positive, excluding
2226 any nifti filename extension.
2227 fname input | return | result of nifti_makebasename
2228 ====================================================================
2229 "myimage" | 1 | "myimage"
2230 "myimage.tif" | 1 | "myimage.tif"
2231 "myimage.tif.gz" | 1 | "myimage.tif"
2232 "myimage.nii" | 1 | "myimage"
2233 ".nii" | 0 | <ERROR - basename has zero length>
2234 ".myhiddenimage" | 1 | ".myhiddenimage"
2235 ".myhiddenimage.nii | 1 | ".myhiddenimage"
2236 *//*--------------------------------------------------------------------*/
nifti_validfilename(const char * fname)2237 int nifti_validfilename(const char* fname)
2238 {
2239 char * ext;
2240
2241 /* check input file(s) for sanity */
2242 if( fname == NULL || *fname == '\0' ){
2243 if ( g_opts.debug > 1 )
2244 fprintf(stderr,"-- empty filename in nifti_validfilename()\n");
2245 return 0;
2246 }
2247
2248 ext = nifti_find_file_extension(fname);
2249
2250 if ( ext && ext == fname ) { /* then no filename prefix */
2251 if ( g_opts.debug > 0 )
2252 fprintf(stderr,"-- no prefix for filename '%s'\n", fname);
2253 return 0;
2254 }
2255
2256 return 1;
2257 }
2258
2259 /*----------------------------------------------------------------------*/
2260 /*! check the end of the filename for a valid nifti extension
2261
2262 Valid extensions are currently .nii, .hdr, .img, .nia,
2263 or any of them followed by .gz. Note that '.' is part of
2264 the extension.
2265
2266 \return a pointer to the extension (within the filename), or NULL
2267 *//*--------------------------------------------------------------------*/
nifti_find_file_extension(const char * name)2268 char * nifti_find_file_extension( const char * name )
2269 {
2270 char * ext;
2271 int len;
2272
2273 if ( ! name ) return NULL;
2274
2275 len = strlen(name);
2276 if ( len < 4 ) return NULL;
2277
2278 ext = (char *)name + len - 4;
2279
2280 if ( (strcmp(ext, ".hdr") == 0) || (strcmp(ext, ".img") == 0) ||
2281 (strcmp(ext, ".nia") == 0) || (strcmp(ext, ".nii") == 0) )
2282 return ext;
2283
2284 #ifdef HAVE_ZLIB
2285 if ( len < 7 ) return NULL;
2286
2287 ext = (char *)name + len - 7;
2288
2289 if ( (strcmp(ext, ".hdr.gz") == 0) || (strcmp(ext, ".img.gz") == 0) ||
2290 (strcmp(ext, ".nia.gz") == 0) || (strcmp(ext, ".nii.gz") == 0) )
2291 return ext;
2292 #endif
2293
2294 if( g_opts.debug > 1 )
2295 fprintf(stderr,"** find_file_ext: failed for name '%s'\n", name);
2296
2297 return NULL;
2298 }
2299
2300 /*----------------------------------------------------------------------*/
2301 /*! return whether the filename ends in ".gz"
2302 *//*--------------------------------------------------------------------*/
nifti_is_gzfile(const char * fname)2303 int nifti_is_gzfile(const char* fname)
2304 {
2305 /* return true if the filename ends with .gz */
2306 if (fname == NULL) { return 0; }
2307 #ifdef HAVE_ZLIB
2308 { /* just so len doesn't generate compile warning */
2309 int len;
2310 len = strlen(fname);
2311 if (len < 3) return 0; /* so we don't search before the name */
2312 if (strcmp(fname + strlen(fname) - 3,".gz")==0) { return 1; }
2313 }
2314 #endif
2315 return 0;
2316 }
2317
2318
2319 /*----------------------------------------------------------------------*/
2320 /*! duplicate the filename, while clearing any extension
2321
2322 This allocates memory for basename which should eventually be freed.
2323 *//*--------------------------------------------------------------------*/
nifti_makebasename(const char * fname)2324 char * nifti_makebasename(const char* fname)
2325 {
2326 char *basename, *ext;
2327
2328 basename=nifti_strdup(fname);
2329
2330 ext = nifti_find_file_extension(basename);
2331 if ( ext ) *ext = '\0'; /* clear out extension */
2332
2333 return basename; /* in either case */
2334 }
2335
2336 /*----------------------------------------------------------------------*/
2337 /*! set nifti's global debug level, for status reporting
2338
2339 - 0 : quiet, nothing is printed to the terminal, but errors
2340 - 1 : normal execution (the default)
2341 - 2, 3 : more details
2342 *//*--------------------------------------------------------------------*/
nifti_set_debug_level(int level)2343 void nifti_set_debug_level( int level )
2344 {
2345 g_opts.debug = level;
2346 }
2347
2348 /*----------------------------------------------------------------------*/
2349 /*! set nifti's global skip_blank_ext flag 5 Sep 2006 [rickr]
2350
2351 explicitly set to 0 or 1
2352 *//*--------------------------------------------------------------------*/
nifti_set_skip_blank_ext(int skip)2353 void nifti_set_skip_blank_ext( int skip )
2354 {
2355 g_opts.skip_blank_ext = skip ? 1 : 0;
2356 }
2357
2358 /*----------------------------------------------------------------------*/
2359 /*! check current directory for existing header file
2360
2361 \return filename of header on success and NULL if no appropriate file
2362 could be found
2363
2364 NB: it allocates memory for hdrname which should be freed
2365 when no longer required
2366 *//*-------------------------------------------------------------------*/
nifti_findhdrname(const char * fname)2367 char * nifti_findhdrname(const char* fname)
2368 {
2369 char *basename, *hdrname, *ext;
2370 char elist[2][5] = { ".hdr", ".nii" };
2371 int efirst;
2372
2373 /**- check input file(s) for sanity */
2374 if( !nifti_validfilename(fname) ) return NULL;
2375
2376 basename = nifti_makebasename(fname);
2377 if( !basename ) return NULL; /* only on string alloc failure */
2378
2379 /**- return filename if it has a valid extension and exists
2380 (except if it is an .img file (and maybe .gz)) */
2381 ext = nifti_find_file_extension(fname);
2382 if ( ext && nifti_fileexists(fname) ) {
2383 if ( strncmp(ext,".img",4) != 0 ){
2384 hdrname = nifti_strdup(fname);
2385 free(basename);
2386 return hdrname;
2387 }
2388 }
2389
2390 /* So the requested name is a basename, contains .img, or does not exist. */
2391 /* In any case, use basename. */
2392
2393 /**- if .img, look for .hdr, .hdr.gz, .nii, .nii.gz, in that order */
2394 /**- else, look for .nii, .nii.gz, .hdr, .hdr.gz, in that order */
2395
2396 /* if we get more extension choices, this could be a loop */
2397
2398 if ( ext && strncmp(ext,".img",4) == 0 ) efirst = 0;
2399 else efirst = 1;
2400
2401 hdrname = (char *)calloc(sizeof(char),strlen(basename)+8);
2402 if( !hdrname ){
2403 fprintf(stderr,"** nifti_findhdrname: failed to alloc hdrname\n");
2404 free(basename);
2405 return NULL;
2406 }
2407
2408 strcpy(hdrname,basename);
2409 strcat(hdrname,elist[efirst]);
2410 if (nifti_fileexists(hdrname)) { free(basename); return hdrname; }
2411 #ifdef HAVE_ZLIB
2412 strcat(hdrname,".gz");
2413 if (nifti_fileexists(hdrname)) { free(basename); return hdrname; }
2414 #endif
2415
2416 /* okay, try the other possibility */
2417
2418 efirst = 1 - efirst;
2419
2420 strcpy(hdrname,basename);
2421 strcat(hdrname,elist[efirst]);
2422 if (nifti_fileexists(hdrname)) { free(basename); return hdrname; }
2423 #ifdef HAVE_ZLIB
2424 strcat(hdrname,".gz");
2425 if (nifti_fileexists(hdrname)) { free(basename); return hdrname; }
2426 #endif
2427
2428 /**- if nothing has been found, return NULL */
2429 free(basename);
2430 free(hdrname);
2431 return NULL;
2432 }
2433
2434
2435 /*------------------------------------------------------------------------*/
2436 /*! check current directory for existing image file
2437
2438 \param fname filename to check for
2439 \nifti_type nifti_type for dataset - this determines whether to
2440 first check for ".nii" or ".img" (since both may exist)
2441
2442 \return filename of data/img file on success and NULL if no appropriate
2443 file could be found
2444
2445 NB: it allocates memory for the image filename, which should be freed
2446 when no longer required
2447 *//*---------------------------------------------------------------------*/
nifti_findimgname(const char * fname,int nifti_type)2448 char * nifti_findimgname(const char* fname , int nifti_type)
2449 {
2450 char *basename, *imgname, ext[2][5] = { ".nii", ".img" };
2451 int first; /* first extension to use */
2452
2453 /* check input file(s) for sanity */
2454
2455 if( !nifti_validfilename(fname) ) return NULL;
2456
2457 basename = nifti_makebasename(fname);
2458 imgname = (char *)calloc(sizeof(char),strlen(basename)+8);
2459 if( !imgname ){
2460 fprintf(stderr,"** nifti_findimgname: failed to alloc imgname\n");
2461 free(basename);
2462 return NULL;
2463 }
2464
2465 /**- test for .nii and .img (don't assume input type from image type) */
2466 /**- if nifti_type = 1, check for .nii first, else .img first */
2467
2468 /* if we get 3 or more extensions, can make a loop here... */
2469
2470 if (nifti_type == NIFTI_FTYPE_NIFTI1_1) first = 0; /* should match .nii */
2471 else first = 1; /* should match .img */
2472
2473 strcpy(imgname,basename);
2474 strcat(imgname,ext[first]);
2475 if (nifti_fileexists(imgname)) { free(basename); return imgname; }
2476 #ifdef HAVE_ZLIB /* then also check for .gz */
2477 strcat(imgname,".gz");
2478 if (nifti_fileexists(imgname)) { free(basename); return imgname; }
2479 #endif
2480
2481 /* failed to find image file with expected extension, try the other */
2482
2483 strcpy(imgname,basename);
2484 strcat(imgname,ext[1-first]); /* can do this with only 2 choices */
2485 if (nifti_fileexists(imgname)) { free(basename); return imgname; }
2486 #ifdef HAVE_ZLIB /* then also check for .gz */
2487 strcat(imgname,".gz");
2488 if (nifti_fileexists(imgname)) { free(basename); return imgname; }
2489 #endif
2490
2491 /**- if nothing has been found, return NULL */
2492 free(basename);
2493 free(imgname);
2494 return NULL;
2495 }
2496
2497
2498 /*----------------------------------------------------------------------*/
2499 /*! creates a filename for storing the header, based on nifti_type
2500
2501 \param prefix - this will be copied before the suffix is added
2502 \param nifti_type - determines the extension, unless one is in prefix
2503 \param check - check for existence (fail condition)
2504 \param comp - add .gz for compressed name
2505
2506 Note that if prefix provides a file suffix, nifti_type is not used.
2507
2508 NB: this allocates memory which should be freed
2509
2510 \sa nifti_set_filenames
2511 *//*-------------------------------------------------------------------*/
nifti_makehdrname(const char * prefix,int nifti_type,int check,int comp)2512 char * nifti_makehdrname(const char * prefix, int nifti_type, int check,
2513 int comp)
2514 {
2515 char * iname, * ext;
2516
2517 if( !nifti_validfilename(prefix) ) return NULL;
2518
2519 /* add space for extension, optional ".gz", and null char */
2520 iname = (char *)calloc(sizeof(char),strlen(prefix)+8);
2521 if( !iname ){ fprintf(stderr,"** small malloc failure!\n"); return NULL; }
2522 strcpy(iname, prefix);
2523
2524 /* use any valid extension */
2525 if( (ext = nifti_find_file_extension(iname)) != NULL ){
2526 if( strncmp(ext,".img",4) == 0 )
2527 memcpy(ext,".hdr",4); /* then convert img name to hdr */
2528 }
2529 /* otherwise, make one up an */
2530 else if( nifti_type == NIFTI_FTYPE_NIFTI1_1 ) strcat(iname, ".nii");
2531 else if( nifti_type == NIFTI_FTYPE_ASCII ) strcat(iname, ".nia");
2532 else strcat(iname, ".hdr");
2533
2534 #ifdef HAVE_ZLIB /* then also check for .gz */
2535 if( comp && (!ext || !strstr(iname,".gz")) ) strcat(iname,".gz");
2536 #endif
2537
2538 /* check for existence failure */
2539 if( check && nifti_fileexists(iname) ){
2540 fprintf(stderr,"** failure: header file '%s' already exists\n",iname);
2541 free(iname);
2542 return NULL;
2543 }
2544
2545 if(g_opts.debug > 2) fprintf(stderr,"+d made header filename '%s'\n", iname);
2546
2547 return iname;
2548 }
2549
2550
2551 /*----------------------------------------------------------------------*/
2552 /*! creates a filename for storing the image, based on nifti_type
2553
2554 \param prefix - this will be copied before the suffix is added
2555 \param nifti_type - determines the extension, unless provided by prefix
2556 \param check - check for existence (fail condition)
2557 \param comp - add .gz for compressed name
2558
2559 Note that if prefix provides a file suffix, nifti_type is not used.
2560
2561 NB: it allocates memory which should be freed
2562
2563 \sa nifti_set_filenames
2564 *//*-------------------------------------------------------------------*/
nifti_makeimgname(const char * prefix,int nifti_type,int check,int comp)2565 char * nifti_makeimgname(const char * prefix, int nifti_type, int check,
2566 int comp)
2567 {
2568 char * iname, * ext;
2569
2570 if( !nifti_validfilename(prefix) ) return NULL;
2571
2572 /* add space for extension, optional ".gz", and null char */
2573 iname = (char *)calloc(sizeof(char),strlen(prefix)+8);
2574 if( !iname ){ fprintf(stderr,"** small malloc failure!\n"); return NULL; }
2575 strcpy(iname, prefix);
2576
2577 /* use any valid extension */
2578 if( (ext = nifti_find_file_extension(iname)) != NULL ){
2579 if( strncmp(ext,".hdr",4) == 0 )
2580 memcpy(ext,".img",4); /* then convert hdr name to img */
2581 }
2582 /* otherwise, make one up */
2583 else if( nifti_type == NIFTI_FTYPE_NIFTI1_1 ) strcat(iname, ".nii");
2584 else if( nifti_type == NIFTI_FTYPE_ASCII ) strcat(iname, ".nia");
2585 else strcat(iname, ".img");
2586
2587 #ifdef HAVE_ZLIB /* then also check for .gz */
2588 if( comp && (!ext || !strstr(iname,".gz")) ) strcat(iname,".gz");
2589 #endif
2590
2591 /* check for existence failure */
2592 if( check && nifti_fileexists(iname) ){
2593 fprintf(stderr,"** failure: image file '%s' already exists\n",iname);
2594 free(iname);
2595 return NULL;
2596 }
2597
2598 if( g_opts.debug > 2 ) fprintf(stderr,"+d made image filename '%s'\n",iname);
2599
2600 return iname;
2601 }
2602
2603
2604 /*----------------------------------------------------------------------*/
2605 /*! create and set new filenames, based on prefix and image type
2606
2607 \param nim pointer to nifti_image in which to set filenames
2608 \param prefix (required) prefix for output filenames
2609 \param check check for previous existence of filename
2610 (existence is an error condition)
2611 \param set_byte_order flag to set nim->byteorder here
2612 (this is probably a logical place to do so)
2613
2614 \return 0 on successful update
2615
2616 \warning this will free() any existing names and create new ones
2617
2618 \sa nifti_makeimgname, nifti_makehdrname, nifti_type_and_names_match
2619 *//*--------------------------------------------------------------------*/
nifti_set_filenames(nifti_image * nim,const char * prefix,int check,int set_byte_order)2620 int nifti_set_filenames( nifti_image * nim, const char * prefix, int check,
2621 int set_byte_order )
2622 {
2623 int comp = nifti_is_gzfile(prefix);
2624
2625 if( !nim || !prefix ){
2626 fprintf(stderr,"** nifti_set_filenames, bad params %p, %p\n",
2627 (void *)nim,prefix);
2628 return -1;
2629 }
2630
2631 if( g_opts.debug > 1 )
2632 fprintf(stderr,"+d modifying output filenames using prefix %s\n", prefix);
2633
2634 if( nim->fname ) free(nim->fname);
2635 if( nim->iname ) free(nim->iname);
2636 nim->fname = nifti_makehdrname(prefix, nim->nifti_type, check, comp);
2637 nim->iname = nifti_makeimgname(prefix, nim->nifti_type, check, comp);
2638 if( !nim->fname || !nim->iname ){
2639 LNI_FERR("nifti_set_filename","failed to set prefix for",prefix);
2640 return -1;
2641 }
2642
2643 if( set_byte_order ) nim->byteorder = nifti_short_order() ;
2644
2645 if( nifti_set_type_from_names(nim) < 0 )
2646 return -1;
2647
2648 if( g_opts.debug > 2 )
2649 fprintf(stderr,"+d have new filenames %s and %s\n",nim->fname,nim->iname);
2650
2651 return 0;
2652 }
2653
2654
2655 /*--------------------------------------------------------------------------*/
2656 /*! check whether nifti_type matches fname and iname for the nifti_image
2657
2658 - if type 0 or 2, expect .hdr/.img pair
2659 - if type 1, expect .nii (and names must match)
2660
2661 \param nim given nifti_image
2662 \param show_warn if set, print a warning message for any mis-match
2663
2664 \return
2665 - 1 if the values seem to match
2666 - 0 if there is a mis-match
2667 - -1 if there is not sufficient information to create file(s)
2668
2669 \sa NIFTI_FTYPE_* codes in nifti1_io.h
2670 \sa nifti_set_type_from_names, is_valid_nifti_type
2671 *//*------------------------------------------------------------------------*/
nifti_type_and_names_match(nifti_image * nim,int show_warn)2672 int nifti_type_and_names_match( nifti_image * nim, int show_warn )
2673 {
2674 char func[] = "nifti_type_and_names_match";
2675 char * ext_h, * ext_i; /* header and image filename extensions */
2676 int errs = 0; /* error counter */
2677
2678 /* sanity checks */
2679 if( !nim ){
2680 if( show_warn ) fprintf(stderr,"** %s: missing nifti_image\n", func);
2681 return -1;
2682 }
2683 if( !nim->fname ){
2684 if( show_warn ) fprintf(stderr,"** %s: missing header filename\n", func);
2685 errs++;
2686 }
2687 if( !nim->iname ){
2688 if( show_warn ) fprintf(stderr,"** %s: missing image filename\n", func);
2689 errs++;
2690 }
2691 if( !is_valid_nifti_type(nim->nifti_type) ){
2692 if( show_warn )
2693 fprintf(stderr,"** %s: bad nifti_type %d\n", func, nim->nifti_type);
2694 errs++;
2695 }
2696
2697 if( errs ) return -1; /* then do not proceed */
2698
2699 /* get pointers to extensions */
2700 ext_h = nifti_find_file_extension( nim->fname );
2701 ext_i = nifti_find_file_extension( nim->iname );
2702
2703 /* check for filename extensions */
2704 if( !ext_h ){
2705 if( show_warn )
2706 fprintf(stderr,"-d missing NIFTI extension in header filename, %s\n",
2707 nim->fname);
2708 errs++;
2709 }
2710 if( !ext_i ){
2711 if( show_warn )
2712 fprintf(stderr,"-d missing NIFTI extension in image filename, %s\n",
2713 nim->iname);
2714 errs++;
2715 }
2716
2717 if( errs ) return 0; /* do not proceed, but this is just a mis-match */
2718
2719 /* general tests */
2720 if( nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ){ /* .nii */
2721 if( strncmp(ext_h,".nii",4) ) {
2722 if( show_warn )
2723 fprintf(stderr,
2724 "-d NIFTI_FTYPE 1, but no .nii extension in header filename, %s\n",
2725 nim->fname);
2726 errs++;
2727 }
2728 if( strncmp(ext_i,".nii",4) ) {
2729 if( show_warn )
2730 fprintf(stderr,
2731 "-d NIFTI_FTYPE 1, but no .nii extension in image filename, %s\n",
2732 nim->iname);
2733 errs++;
2734 }
2735 if( strcmp(nim->fname, nim->iname) != 0 ){
2736 if( show_warn )
2737 fprintf(stderr,
2738 "-d NIFTI_FTYPE 1, but header and image filenames differ: %s, %s\n",
2739 nim->fname, nim->iname);
2740 errs++;
2741 }
2742 }
2743 else if( (nim->nifti_type == NIFTI_FTYPE_NIFTI1_2) || /* .hdr/.img */
2744 (nim->nifti_type == NIFTI_FTYPE_ANALYZE) )
2745 {
2746 if( strncmp(ext_h,".hdr",4) != 0 ){
2747 if( show_warn )
2748 fprintf(stderr,"-d no '.hdr' extension, but NIFTI type is %d, %s\n",
2749 nim->nifti_type, nim->fname);
2750 errs++;
2751 }
2752 if( strncmp(ext_i,".img",4) != 0 ){
2753 if( show_warn )
2754 fprintf(stderr,"-d no '.img' extension, but NIFTI type is %d, %s\n",
2755 nim->nifti_type, nim->iname);
2756 errs++;
2757 }
2758 }
2759 /* ignore any other nifti_type */
2760
2761 return 1;
2762 }
2763
2764
2765 /*--------------------------------------------------------------------------*/
2766 /*! check whether the given type is on the "approved" list
2767
2768 The code is valid if it is non-negative, and does not exceed
2769 NIFTI_MAX_FTYPE.
2770
2771 \return 1 if nifti_type is valid, 0 otherwise
2772 \sa NIFTI_FTYPE_* codes in nifti1_io.h
2773 *//*------------------------------------------------------------------------*/
is_valid_nifti_type(int nifti_type)2774 int is_valid_nifti_type( int nifti_type )
2775 {
2776 if( nifti_type >= NIFTI_FTYPE_ANALYZE && /* smallest type, 0 */
2777 nifti_type <= NIFTI_MAX_FTYPE )
2778 return 1;
2779 return 0;
2780 }
2781
2782
2783 /*--------------------------------------------------------------------------*/
2784 /*! set the nifti_type field based on fname and iname
2785
2786 Note that nifti_type is changed only when it does not match
2787 the filenames.
2788
2789 \return 0 on success, -1 on error
2790
2791 \sa is_valid_nifti_type, nifti_type_and_names_match
2792 *//*------------------------------------------------------------------------*/
nifti_set_type_from_names(nifti_image * nim)2793 int nifti_set_type_from_names( nifti_image * nim )
2794 {
2795 /* error checking first */
2796 if( !nim ){ fprintf(stderr,"** NSTFN: no nifti_image\n"); return -1; }
2797
2798 if( !nim->fname || !nim->iname ){
2799 fprintf(stderr,"** NSTFN: missing filename(s) fname @ %p, iname @ %p\n",
2800 nim->fname, nim->iname);
2801 return -1;
2802 }
2803
2804 if( ! nifti_validfilename ( nim->fname ) ||
2805 ! nifti_validfilename ( nim->iname ) ||
2806 ! nifti_find_file_extension( nim->fname ) ||
2807 ! nifti_find_file_extension( nim->iname )
2808 ) {
2809 fprintf(stderr,"** NSTFN: invalid filename(s) fname='%s', iname='%s'\n",
2810 nim->fname, nim->iname);
2811 return -1;
2812 }
2813
2814 if( g_opts.debug > 2 )
2815 fprintf(stderr,"-d verify nifti_type from filenames: %d",nim->nifti_type);
2816
2817 /* not too picky here, do what must be done, and then verify */
2818 if( strcmp(nim->fname, nim->iname) == 0 ) /* one file, type 1 */
2819 nim->nifti_type = NIFTI_FTYPE_NIFTI1_1;
2820 else if( nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ) /* cannot be type 1 */
2821 nim->nifti_type = NIFTI_FTYPE_NIFTI1_2;
2822
2823 if( g_opts.debug > 2 ) fprintf(stderr," -> %d\n",nim->nifti_type);
2824
2825 if( g_opts.debug > 1 ) /* warn user about anything strange */
2826 nifti_type_and_names_match(nim, 1);
2827
2828 if( is_valid_nifti_type(nim->nifti_type) ) return 0; /* success! */
2829
2830 fprintf(stderr,"** NSTFN: bad nifti_type %d, for '%s' and '%s'\n",
2831 nim->nifti_type, nim->fname, nim->iname);
2832
2833 return -1;
2834 }
2835
2836
2837 /*--------------------------------------------------------------------------*/
2838 /*! Determine if this is a NIFTI-formatted file.
2839
2840 <pre>
2841 \return 0 if file looks like ANALYZE 7.5 [checks sizeof_hdr field == 348]
2842 1 if file marked as NIFTI (header+data in 1 file)
2843 2 if file marked as NIFTI (header+data in 2 files)
2844 -1 if it can't tell, file doesn't exist, etc.
2845 </pre>
2846 *//*------------------------------------------------------------------------*/
is_nifti_file(const char * hname)2847 int is_nifti_file( const char *hname )
2848 {
2849 struct nifti_1_header nhdr ;
2850 znzFile fp ;
2851 int ii ;
2852 char *tmpname;
2853
2854 /* bad input name? */
2855
2856 if( !nifti_validfilename(hname) ) return -1 ;
2857
2858 /* open file */
2859
2860 tmpname = nifti_findhdrname(hname);
2861 if( tmpname == NULL ){
2862 if( g_opts.debug > 0 )
2863 fprintf(stderr,"** no header file found for '%s'\n",hname);
2864 return -1;
2865 }
2866 fp = znzopen( tmpname , "rb" , nifti_is_gzfile(tmpname) ) ;
2867 free(tmpname);
2868 if (znz_isnull(fp)) return -1 ; /* bad open? */
2869
2870 /* read header, close file */
2871
2872 ii = znzread( &nhdr , 1 , sizeof(nhdr) , fp ) ;
2873 znzclose( fp ) ;
2874 if( ii < (int) sizeof(nhdr) ) return -1 ; /* bad read? */
2875
2876 /* check for NIFTI-ness */
2877
2878 if( NIFTI_VERSION(nhdr) != 0 ){
2879 return ( NIFTI_ONEFILE(nhdr) ) ? 1 : 2 ;
2880 }
2881
2882 /* check for ANALYZE-ness (sizeof_hdr field == 348) */
2883
2884 ii = nhdr.sizeof_hdr ;
2885 if( ii == (int)sizeof(nhdr) ) return 0 ; /* matches */
2886
2887 /* try byte-swapping header */
2888
2889 swap_4(ii) ;
2890 if( ii == (int)sizeof(nhdr) ) return 0 ; /* matches */
2891
2892 return -1 ; /* not good */
2893 }
2894
print_hex_vals(const char * data,int nbytes,FILE * fp)2895 static int print_hex_vals( const char * data, int nbytes, FILE * fp )
2896 {
2897 int c;
2898
2899 if ( !data || nbytes < 1 || !fp ) return -1;
2900
2901 fputs("0x", fp);
2902 for ( c = 0; c < nbytes; c++ )
2903 fprintf(fp, " %x", data[c]);
2904
2905 return 0;
2906 }
2907
2908 /*----------------------------------------------------------------------*/
2909 /*! display the contents of the nifti_1_header (send to stdout)
2910
2911 \param info if non-NULL, print this character string
2912 \param hp pointer to nifti_1_header
2913 *//*--------------------------------------------------------------------*/
disp_nifti_1_header(const char * info,const nifti_1_header * hp)2914 int disp_nifti_1_header( const char * info, const nifti_1_header * hp )
2915 {
2916 int c;
2917
2918 fputs( "-------------------------------------------------------\n", stdout );
2919 if ( info ) fputs( info, stdout );
2920 if ( !hp ){ fputs(" ** no nifti_1_header to display!\n",stdout); return 1; }
2921
2922 fprintf(stdout," nifti_1_header :\n"
2923 " sizeof_hdr = %d\n"
2924 " data_type[10] = ", hp->sizeof_hdr);
2925 print_hex_vals(hp->data_type, 10, stdout);
2926 fprintf(stdout, "\n"
2927 " db_name[18] = ");
2928 print_hex_vals(hp->db_name, 18, stdout);
2929 fprintf(stdout, "\n"
2930 " extents = %d\n"
2931 " session_error = %d\n"
2932 " regular = 0x%x\n"
2933 " dim_info = 0x%x\n",
2934 hp->extents, hp->session_error, hp->regular, hp->dim_info );
2935 fprintf(stdout, " dim[8] =");
2936 for ( c = 0; c < 8; c++ ) fprintf(stdout," %d", hp->dim[c]);
2937 fprintf(stdout, "\n"
2938 " intent_p1 = %f\n"
2939 " intent_p2 = %f\n"
2940 " intent_p3 = %f\n"
2941 " intent_code = %d\n"
2942 " datatype = %d\n"
2943 " bitpix = %d\n"
2944 " slice_start = %d\n"
2945 " pixdim[8] =",
2946 hp->intent_p1, hp->intent_p2, hp->intent_p3, hp->intent_code,
2947 hp->datatype, hp->bitpix, hp->slice_start);
2948 /* break pixdim over 2 lines */
2949 for ( c = 0; c < 4; c++ ) fprintf(stdout," %f", hp->pixdim[c]);
2950 fprintf(stdout, "\n ");
2951 for ( c = 4; c < 8; c++ ) fprintf(stdout," %f", hp->pixdim[c]);
2952 fprintf(stdout, "\n"
2953 " vox_offset = %f\n"
2954 " scl_slope = %f\n"
2955 " scl_inter = %f\n"
2956 " slice_end = %d\n"
2957 " slice_code = %d\n"
2958 " xyzt_units = 0x%x\n"
2959 " cal_max = %f\n"
2960 " cal_min = %f\n"
2961 " slice_duration = %f\n"
2962 " toffset = %f\n"
2963 " glmax = %d\n"
2964 " glmin = %d\n",
2965 hp->vox_offset, hp->scl_slope, hp->scl_inter, hp->slice_end,
2966 hp->slice_code, hp->xyzt_units, hp->cal_max, hp->cal_min,
2967 hp->slice_duration, hp->toffset, hp->glmax, hp->glmin);
2968 fprintf(stdout,
2969 " descrip = '%.80s'\n"
2970 " aux_file = '%.24s'\n"
2971 " qform_code = %d\n"
2972 " sform_code = %d\n"
2973 " quatern_b = %f\n"
2974 " quatern_c = %f\n"
2975 " quatern_d = %f\n"
2976 " qoffset_x = %f\n"
2977 " qoffset_y = %f\n"
2978 " qoffset_z = %f\n"
2979 " srow_x[4] = %f, %f, %f, %f\n"
2980 " srow_y[4] = %f, %f, %f, %f\n"
2981 " srow_z[4] = %f, %f, %f, %f\n"
2982 " intent_name = '%-.16s'\n"
2983 " magic = '%-.4s'\n",
2984 hp->descrip, hp->aux_file, hp->qform_code, hp->sform_code,
2985 hp->quatern_b, hp->quatern_c, hp->quatern_d,
2986 hp->qoffset_x, hp->qoffset_y, hp->qoffset_z,
2987 hp->srow_x[0], hp->srow_x[1], hp->srow_x[2], hp->srow_x[3],
2988 hp->srow_y[0], hp->srow_y[1], hp->srow_y[2], hp->srow_y[3],
2989 hp->srow_z[0], hp->srow_z[1], hp->srow_z[2], hp->srow_z[3],
2990 hp->intent_name, hp->magic);
2991 fputs( "-------------------------------------------------------\n", stdout );
2992 fflush(stdout);
2993
2994 return 0;
2995 }
2996
2997
2998 #undef ERREX
2999 #define ERREX(msg) \
3000 do{ fprintf(stderr,"** ERROR: nifti_convert_nhdr2nim: %s\n", (msg) ) ; \
3001 return NULL ; } while(0)
3002
3003 /*----------------------------------------------------------------------*/
3004 /*! convert a nifti_1_header into a nift1_image
3005
3006 \return an allocated nifti_image, or NULL on failure
3007 *//*--------------------------------------------------------------------*/
nifti_convert_nhdr2nim(struct nifti_1_header nhdr,const char * fname)3008 nifti_image* nifti_convert_nhdr2nim(struct nifti_1_header nhdr,
3009 const char * fname)
3010 {
3011 int ii , doswap , ioff ;
3012 int is_nifti , is_onefile ;
3013 nifti_image *nim;
3014
3015 nim = (nifti_image *) calloc( 1 , sizeof(nifti_image) ) ;
3016 if( !nim ) ERREX("failed to allocate nifti image");
3017
3018 /* be explicit with pointers */
3019 nim->fname = NULL;
3020 nim->iname = NULL;
3021 nim->data = NULL;
3022
3023 /**- check if we must swap bytes */
3024
3025 doswap = need_nhdr_swap(nhdr.dim[0], nhdr.sizeof_hdr); /* swap data flag */
3026
3027 if( doswap < 0 ){
3028 if( doswap == -1 ) ERREX("bad dim[0]") ;
3029 ERREX("bad sizeof_hdr") ; /* else */
3030 }
3031
3032 /**- determine if this is a NIFTI-1 compliant header */
3033
3034 is_nifti = NIFTI_VERSION(nhdr) ;
3035 if( doswap ) {
3036 if ( g_opts.debug > 3 ) disp_nifti_1_header("-d ni1 pre-swap: ", &nhdr);
3037 swap_nifti_header( &nhdr , is_nifti ) ;
3038 }
3039
3040 if ( g_opts.debug > 2 ) disp_nifti_1_header("-d nhdr2nim : ", &nhdr);
3041
3042 if( nhdr.datatype == DT_BINARY ||
3043 nhdr.datatype == DT_UNKNOWN ) ERREX("bad datatype") ;
3044
3045 if( nhdr.dim[1] <= 0 ) ERREX("bad dim[1]") ;
3046
3047 /* fix bad dim[] values in the defined dimension range */
3048 for( ii=2 ; ii <= nhdr.dim[0] ; ii++ )
3049 if( nhdr.dim[ii] <= 0 ) nhdr.dim[ii] = 1 ;
3050
3051 /* fix any remaining bad dim[] values, so garbage does not propagate */
3052 /* (only values 0 or 1 seem rational, otherwise set to arbirary 1) */
3053 for( ii=nhdr.dim[0]+1 ; ii <= 7 ; ii++ )
3054 if( nhdr.dim[ii] != 1 && nhdr.dim[ii] != 0) nhdr.dim[ii] = 1 ;
3055
3056 #if 0 /* rely on dim[0], do not attempt to modify it 16 Nov 2005 [rickr] */
3057
3058 /**- get number of dimensions (ignoring dim[0] now) */
3059 for( ii=7 ; ii >= 2 ; ii-- ) /* loop backwards until we */
3060 if( nhdr.dim[ii] > 1 ) break ; /* find a dim bigger than 1 */
3061 ndim = ii ;
3062 #endif
3063
3064 /**- set bad grid spacings to 1.0 */
3065
3066 for( ii=1 ; ii <= nhdr.dim[0] ; ii++ ){
3067 if( nhdr.pixdim[ii] == 0.0 ||
3068 !IS_GOOD_FLOAT(nhdr.pixdim[ii]) ) nhdr.pixdim[ii] = 1.0 ;
3069 }
3070
3071 is_onefile = is_nifti && NIFTI_ONEFILE(nhdr) ;
3072
3073 if( is_nifti ) nim->nifti_type = (is_onefile) ? NIFTI_FTYPE_NIFTI1_1
3074 : NIFTI_FTYPE_NIFTI1_2 ;
3075 else nim->nifti_type = NIFTI_FTYPE_ANALYZE ;
3076
3077 ii = nifti_short_order() ;
3078 if( doswap ) nim->byteorder = REVERSE_ORDER(ii) ;
3079 else nim->byteorder = ii ;
3080
3081
3082 /**- set dimensions of data array */
3083
3084 nim->ndim = nim->dim[0] = nhdr.dim[0];
3085 nim->nx = nim->dim[1] = nhdr.dim[1];
3086 nim->ny = nim->dim[2] = nhdr.dim[2];
3087 nim->nz = nim->dim[3] = nhdr.dim[3];
3088 nim->nt = nim->dim[4] = nhdr.dim[4];
3089 nim->nu = nim->dim[5] = nhdr.dim[5];
3090 nim->nv = nim->dim[6] = nhdr.dim[6];
3091 nim->nw = nim->dim[7] = nhdr.dim[7];
3092
3093 for( ii=1, nim->nvox=1; ii <= nhdr.dim[0]; ii++ )
3094 nim->nvox *= nhdr.dim[ii];
3095
3096 /**- set the type of data in voxels and how many bytes per voxel */
3097
3098 nim->datatype = nhdr.datatype ;
3099
3100 nifti_datatype_sizes( nim->datatype , &(nim->nbyper) , &(nim->swapsize) ) ;
3101 if( nim->nbyper == 0 ){ free(nim); ERREX("bad datatype"); }
3102
3103 /**- set the grid spacings */
3104
3105 nim->dx = nim->pixdim[1] = nhdr.pixdim[1] ;
3106 nim->dy = nim->pixdim[2] = nhdr.pixdim[2] ;
3107 nim->dz = nim->pixdim[3] = nhdr.pixdim[3] ;
3108 nim->dt = nim->pixdim[4] = nhdr.pixdim[4] ;
3109 nim->du = nim->pixdim[5] = nhdr.pixdim[5] ;
3110 nim->dv = nim->pixdim[6] = nhdr.pixdim[6] ;
3111 nim->dw = nim->pixdim[7] = nhdr.pixdim[7] ;
3112
3113 /**- compute qto_xyz transformation from pixel indexes (i,j,k) to (x,y,z) */
3114
3115 if( !is_nifti || nhdr.qform_code <= 0 ){
3116 /**- if not nifti or qform_code <= 0, use grid spacing for qto_xyz */
3117
3118 nim->qto_xyz.m[0][0] = nim->dx ; /* grid spacings */
3119 nim->qto_xyz.m[1][1] = nim->dy ; /* along diagonal */
3120 nim->qto_xyz.m[2][2] = nim->dz ;
3121
3122 /* off diagonal is zero */
3123
3124 nim->qto_xyz.m[0][1]=nim->qto_xyz.m[0][2]=nim->qto_xyz.m[0][3] = 0.0;
3125 nim->qto_xyz.m[1][0]=nim->qto_xyz.m[1][2]=nim->qto_xyz.m[1][3] = 0.0;
3126 nim->qto_xyz.m[2][0]=nim->qto_xyz.m[2][1]=nim->qto_xyz.m[2][3] = 0.0;
3127
3128 /* last row is always [ 0 0 0 1 ] */
3129
3130 nim->qto_xyz.m[3][0]=nim->qto_xyz.m[3][1]=nim->qto_xyz.m[3][2] = 0.0;
3131 nim->qto_xyz.m[3][3]= 1.0 ;
3132
3133 nim->qform_code = NIFTI_XFORM_UNKNOWN ;
3134
3135 if( g_opts.debug > 1 ) fprintf(stderr,"-d no qform provided\n");
3136 } else {
3137 /**- else NIFTI: use the quaternion-specified transformation */
3138
3139 nim->quatern_b = FIXED_FLOAT( nhdr.quatern_b ) ;
3140 nim->quatern_c = FIXED_FLOAT( nhdr.quatern_c ) ;
3141 nim->quatern_d = FIXED_FLOAT( nhdr.quatern_d ) ;
3142
3143 nim->qoffset_x = FIXED_FLOAT(nhdr.qoffset_x) ;
3144 nim->qoffset_y = FIXED_FLOAT(nhdr.qoffset_y) ;
3145 nim->qoffset_z = FIXED_FLOAT(nhdr.qoffset_z) ;
3146
3147 nim->qfac = (nhdr.pixdim[0] < 0.0) ? -1.0 : 1.0 ; /* left-handedness? */
3148
3149 nim->qto_xyz = nifti_quatern_to_mat44(
3150 nim->quatern_b, nim->quatern_c, nim->quatern_d,
3151 nim->qoffset_x, nim->qoffset_y, nim->qoffset_z,
3152 nim->dx , nim->dy , nim->dz ,
3153 nim->qfac ) ;
3154
3155 nim->qform_code = nhdr.qform_code ;
3156
3157 if( g_opts.debug > 1 )
3158 nifti_disp_matrix_orient("-d qform orientations:\n", nim->qto_xyz);
3159 }
3160
3161 /**- load inverse transformation (x,y,z) -> (i,j,k) */
3162
3163 nim->qto_ijk = nifti_mat44_inverse( nim->qto_xyz ) ;
3164
3165 /**- load sto_xyz affine transformation, if present */
3166
3167 if( !is_nifti || nhdr.sform_code <= 0 ){
3168 /**- if not nifti or sform_code <= 0, then no sto transformation */
3169
3170 nim->sform_code = NIFTI_XFORM_UNKNOWN ;
3171
3172 if( g_opts.debug > 1 ) fprintf(stderr,"-d no sform provided\n");
3173
3174 } else {
3175 /**- else set the sto transformation from srow_*[] */
3176
3177 nim->sto_xyz.m[0][0] = nhdr.srow_x[0] ;
3178 nim->sto_xyz.m[0][1] = nhdr.srow_x[1] ;
3179 nim->sto_xyz.m[0][2] = nhdr.srow_x[2] ;
3180 nim->sto_xyz.m[0][3] = nhdr.srow_x[3] ;
3181
3182 nim->sto_xyz.m[1][0] = nhdr.srow_y[0] ;
3183 nim->sto_xyz.m[1][1] = nhdr.srow_y[1] ;
3184 nim->sto_xyz.m[1][2] = nhdr.srow_y[2] ;
3185 nim->sto_xyz.m[1][3] = nhdr.srow_y[3] ;
3186
3187 nim->sto_xyz.m[2][0] = nhdr.srow_z[0] ;
3188 nim->sto_xyz.m[2][1] = nhdr.srow_z[1] ;
3189 nim->sto_xyz.m[2][2] = nhdr.srow_z[2] ;
3190 nim->sto_xyz.m[2][3] = nhdr.srow_z[3] ;
3191
3192 /* last row is always [ 0 0 0 1 ] */
3193
3194 nim->sto_xyz.m[3][0]=nim->sto_xyz.m[3][1]=nim->sto_xyz.m[3][2] = 0.0;
3195 nim->sto_xyz.m[3][3]= 1.0 ;
3196
3197 nim->sto_ijk = nifti_mat44_inverse( nim->sto_xyz ) ;
3198
3199 nim->sform_code = nhdr.sform_code ;
3200
3201 if( g_opts.debug > 1 )
3202 nifti_disp_matrix_orient("-d sform orientations:\n", nim->sto_xyz);
3203 }
3204
3205 /**- set miscellaneous NIFTI stuff */
3206
3207 if( is_nifti ){
3208 nim->scl_slope = FIXED_FLOAT( nhdr.scl_slope ) ;
3209 nim->scl_inter = FIXED_FLOAT( nhdr.scl_inter ) ;
3210
3211 nim->intent_code = nhdr.intent_code ;
3212
3213 nim->intent_p1 = FIXED_FLOAT( nhdr.intent_p1 ) ;
3214 nim->intent_p2 = FIXED_FLOAT( nhdr.intent_p2 ) ;
3215 nim->intent_p3 = FIXED_FLOAT( nhdr.intent_p3 ) ;
3216
3217 nim->toffset = FIXED_FLOAT( nhdr.toffset ) ;
3218
3219 memcpy(nim->intent_name,nhdr.intent_name,15); nim->intent_name[15] = '\0';
3220
3221 nim->xyz_units = XYZT_TO_SPACE(nhdr.xyzt_units) ;
3222 nim->time_units = XYZT_TO_TIME (nhdr.xyzt_units) ;
3223
3224 nim->freq_dim = DIM_INFO_TO_FREQ_DIM ( nhdr.dim_info ) ;
3225 nim->phase_dim = DIM_INFO_TO_PHASE_DIM( nhdr.dim_info ) ;
3226 nim->slice_dim = DIM_INFO_TO_SLICE_DIM( nhdr.dim_info ) ;
3227
3228 nim->slice_code = nhdr.slice_code ;
3229 nim->slice_start = nhdr.slice_start ;
3230 nim->slice_end = nhdr.slice_end ;
3231 nim->slice_duration = FIXED_FLOAT(nhdr.slice_duration) ;
3232 }
3233
3234 /**- set Miscellaneous ANALYZE stuff */
3235
3236 nim->cal_min = FIXED_FLOAT(nhdr.cal_min) ;
3237 nim->cal_max = FIXED_FLOAT(nhdr.cal_max) ;
3238
3239 memcpy(nim->descrip ,nhdr.descrip ,79) ; nim->descrip [79] = '\0' ;
3240 memcpy(nim->aux_file,nhdr.aux_file,23) ; nim->aux_file[23] = '\0' ;
3241
3242 /**- set ioff from vox_offset (but at least sizeof(header)) */
3243
3244 is_onefile = is_nifti && NIFTI_ONEFILE(nhdr) ;
3245
3246 if( is_onefile ){
3247 ioff = (int)nhdr.vox_offset ;
3248 if( ioff < (int) sizeof(nhdr) ) ioff = (int) sizeof(nhdr) ;
3249 } else {
3250 ioff = (int)nhdr.vox_offset ;
3251 }
3252 nim->iname_offset = ioff ;
3253
3254
3255 /**- deal with file names if set */
3256 if (fname!=NULL) {
3257 nifti_set_filenames(nim,fname,0,0);
3258 if (nim->iname==NULL) { ERREX("bad filename"); }
3259 } else {
3260 nim->fname = NULL;
3261 nim->iname = NULL;
3262 }
3263
3264 /* clear extension fields */
3265 nim->num_ext = 0;
3266 nim->ext_list = NULL;
3267
3268 return nim;
3269 }
3270
3271 #undef ERREX
3272 #define ERREX(msg) \
3273 do{ fprintf(stderr,"** ERROR: nifti_image_open(%s): %s\n", \
3274 (hname != NULL) ? hname : "(null)" , (msg) ) ; \
3275 return fptr ; } while(0)
3276
3277 /***************************************************************
3278 * nifti_image_open
3279 ***************************************************************/
3280 /*! znzFile nifti_image_open( char *hname, char *opts , nifti_image **nim)
3281 \brief Read in NIFTI-1 or ANALYZE-7.5 file (pair) header information into a nifti_image struct.
3282
3283 - The image data is not read from disk (it may be read later using
3284 nifti_image_load(), for example).
3285 - The image data will be stored in whatever data format the
3286 input data is; no scaling will be applied.
3287 - DT_BINARY data is not supported.
3288 - nifti_image_free() can be used to delete the returned struct,
3289 when you are done with it.
3290
3291 \param hname filename of dataset .hdr or .nii file
3292 \param opts options string for opening the header file
3293 \param nim pointer to pointer to nifti_image struct
3294 (this routine allocates the nifti_image struct)
3295 \return file pointer (gzippable) to the file with the image data,
3296 ready for reading.
3297 <br>NULL if something fails badly.
3298 \sa nifti_image_load, nifti_image_free
3299 */
nifti_image_open(const char * hname,char * opts,nifti_image ** nim)3300 znzFile nifti_image_open(const char * hname, char * opts, nifti_image ** nim)
3301 {
3302 znzFile fptr=NULL;
3303 /* open the hdr and reading it in, but do not load the data */
3304 *nim = nifti_image_read(hname,0);
3305 /* open the image file, ready for reading (compressed works for all reads) */
3306 if( ((*nim) == NULL) || ((*nim)->iname == NULL) ||
3307 ((*nim)->nbyper <= 0) || ((*nim)->nvox <= 0) )
3308 ERREX("bad header info") ;
3309
3310 /* open image data file */
3311 fptr = znzopen( (*nim)->iname, opts, nifti_is_gzfile((*nim)->iname) );
3312 if( znz_isnull(fptr) ) ERREX("Can't open data file") ;
3313
3314 return fptr;
3315 }
3316
3317
3318 /*----------------------------------------------------------------------*/
3319 /*! return an allocated and filled nifti_1_header struct
3320
3321 Read the binary header from disk, and swap bytes if necessary.
3322
3323 \return an allocated nifti_1_header struct, or NULL on failure
3324
3325 \param hname name of file containing header
3326 \param swapped if not NULL, return whether header bytes were swapped
3327 \param check flag to check for invalid nifti_1_header
3328
3329 \warning ASCII header type is not supported
3330
3331 \sa nifti_image_read, nifti_image_free, nifti_image_read_bricks
3332 *//*--------------------------------------------------------------------*/
nifti_read_header(const char * hname,int * swapped,int check)3333 nifti_1_header * nifti_read_header(const char * hname, int * swapped, int check)
3334 {
3335 nifti_1_header nhdr, * hptr;
3336 znzFile fp;
3337 int bytes, lswap;
3338 char * hfile;
3339 char fname[] = { "nifti_read_header" };
3340
3341 /* determine file name to use for header */
3342 hfile = nifti_findhdrname(hname);
3343 if( hfile == NULL ){
3344 if( g_opts.debug > 0 )
3345 LNI_FERR(fname,"failed to find header file for", hname);
3346 return NULL;
3347 } else if( g_opts.debug > 1 )
3348 fprintf(stderr,"-d %s: found header filename '%s'\n",fname,hfile);
3349
3350 fp = znzopen( hfile, "rb", nifti_is_gzfile(hfile) );
3351 if( znz_isnull(fp) ){
3352 if( g_opts.debug > 0 ) LNI_FERR(fname,"failed to open header file",hfile);
3353 free(hfile);
3354 return NULL;
3355 }
3356
3357 free(hfile); /* done with filename */
3358
3359 if( has_ascii_header(fp) == 1 ){
3360 znzclose( fp );
3361 if( g_opts.debug > 0 )
3362 LNI_FERR(fname,"ASCII header type not supported",hname);
3363 return NULL;
3364 }
3365
3366 /* read the binary header */
3367 bytes = znzread( &nhdr, 1, sizeof(nhdr), fp );
3368 znzclose( fp ); /* we are done with the file now */
3369
3370 if( bytes < (int)sizeof(nhdr) ){
3371 if( g_opts.debug > 0 ){
3372 LNI_FERR(fname,"bad binary header read for file", hname);
3373 fprintf(stderr," - read %d of %d bytes\n",bytes, (int)sizeof(nhdr));
3374 }
3375 return NULL;
3376 }
3377
3378 /* now just decide on byte swapping */
3379 lswap = need_nhdr_swap(nhdr.dim[0], nhdr.sizeof_hdr); /* swap data flag */
3380 if( check && lswap < 0 ){
3381 LNI_FERR(fname,"bad nifti_1_header for file", hname);
3382 return NULL;
3383 }
3384
3385 if( lswap ) {
3386 if ( g_opts.debug > 3 ) disp_nifti_1_header("-d nhdr pre-swap: ", &nhdr);
3387 swap_nifti_header( &nhdr , NIFTI_VERSION(nhdr) ) ;
3388 }
3389
3390 if ( g_opts.debug > 2 ) disp_nifti_1_header("-d nhdr post-swap: ", &nhdr);
3391
3392 if ( check && ! nifti_hdr_looks_good(&nhdr) ){
3393 LNI_FERR(fname,"nifti_1_header looks bad for file", hname);
3394 return NULL;
3395 }
3396
3397 /* all looks good, so allocate memory for and return the header */
3398 hptr = (nifti_1_header *)malloc(sizeof(nifti_1_header));
3399 if( ! hptr ){
3400 fprintf(stderr,"** nifti_read_hdr: failed to alloc nifti_1_header\n");
3401 return NULL;
3402 }
3403
3404 if( swapped ) *swapped = lswap; /* only if they care <sniff!> */
3405
3406 memcpy(hptr, &nhdr, sizeof(nifti_1_header));
3407
3408 return hptr;
3409 }
3410
3411
3412 /*----------------------------------------------------------------------*/
3413 /*! decide if this nifti_1_header structure looks reasonable
3414
3415 Check dim[0], dim[1], sizeof_hdr, and datatype.
3416 Check magic string for "n+1".
3417 Maybe more tests will follow.
3418
3419 \return 1 if the header seems valid, 0 otherwise
3420
3421 \sa nifti_nim_is_valid, valid_nifti_extensions
3422 *//*--------------------------------------------------------------------*/
nifti_hdr_looks_good(const nifti_1_header * hdr)3423 int nifti_hdr_looks_good(const nifti_1_header * hdr)
3424 {
3425 int nbyper, swapsize;
3426 int c, errs = 0;
3427
3428 /* check dim[0] and sizeof_hdr */
3429 if( need_nhdr_swap(hdr->dim[0], hdr->sizeof_hdr) < 0 ){
3430 if( g_opts.debug > 0 )
3431 fprintf(stderr,"** bad nhdr fields: dim0, sizeof_hdr = %d, %d\n",
3432 hdr->dim[0], hdr->sizeof_hdr);
3433 errs++;
3434 }
3435
3436 /* check the valid dimension sizes (maybe dim[0] is bad) */
3437 for( c = 1; c <= hdr->dim[0] && c <= 7; c++ )
3438 if( hdr->dim[c] <= 0 ){
3439 if( g_opts.debug > 0 )
3440 fprintf(stderr,"** bad nhdr field: dim[%d] = %d\n",c,hdr->dim[c]);
3441 errs++;
3442 }
3443
3444 /* check the magic string */
3445 if( (hdr->magic[0] != 'n') ||
3446 (hdr->magic[1] != 'i' && hdr->magic[1] != '+') ||
3447 (hdr->magic[2] != '1') ||
3448 (hdr->magic[3] != '\0') )
3449 {
3450 if( g_opts.debug > 0 )
3451 fprintf(stderr,
3452 "** bad nhdr field: magic = '%.4s', should be \"n+1\" or \"ni1\"\n"
3453 " (in hex) magic = 0x%02x%02x%02x%02x\n"
3454 " should be = 0x6e2b3100 or 0x6e693100\n",
3455 hdr->magic, hdr->magic[0], hdr->magic[1],
3456 hdr->magic[2], hdr->magic[3]);
3457 errs++;
3458 }
3459
3460 /* check the datatype */
3461 if( hdr->datatype == DT_BINARY || hdr->datatype == DT_UNKNOWN ){
3462 if( g_opts.debug > 0 )
3463 fprintf(stderr,"** bad nhdr field: datatype = %d\n",hdr->datatype);
3464 errs++;
3465 }
3466
3467 nifti_datatype_sizes(hdr->datatype, &nbyper, &swapsize);
3468 if( nbyper == 0 ){
3469 if( g_opts.debug > 0 )
3470 fprintf(stderr,"** bad nhdr field: datatype = %d\n",hdr->datatype);
3471 errs++;
3472 }
3473
3474 if( errs ) return 0; /* problems */
3475
3476 if( g_opts.debug > 2 ) fprintf(stderr,"-d nifti header looks good\n");
3477
3478 return 1; /* looks good */
3479 }
3480
3481
3482 /*----------------------------------------------------------------------
3483 * check whether byte swapping is needed
3484 *
3485 * dim[0] should be in [0,7], and sizeof_hdr should be accurate
3486 *
3487 * \returns > 0 : needs swap
3488 * 0 : does not need swap
3489 * < 0 : error condition
3490 *----------------------------------------------------------------------*/
need_nhdr_swap(short dim0,int hdrsize)3491 static int need_nhdr_swap( short dim0, int hdrsize )
3492 {
3493 short d0 = dim0; /* so we won't have to swap them on the stack */
3494 int hsize = hdrsize;
3495
3496 if( d0 != 0 ){ /* then use it for the check */
3497 if( d0 > 0 && d0 <= 7 ) return 0;
3498
3499 nifti_swap_2bytes(1, &d0); /* swap? */
3500 if( d0 > 0 && d0 <= 7 ) return 1;
3501
3502 if( g_opts.debug > 1 ){
3503 fprintf(stderr,"** bad swapped d0 = %d, unswapped = ", d0);
3504 nifti_swap_2bytes(1, &d0); /* swap? */
3505 fprintf(stderr,"%d\n", d0);
3506 }
3507
3508 return -1; /* bad, naughty d0 */
3509 }
3510
3511 /* dim[0] == 0 should not happen, but could, so try hdrsize */
3512 if( hsize == sizeof(nifti_1_header) ) return 0;
3513
3514 nifti_swap_4bytes(1, &hsize); /* swap? */
3515 if( hsize == sizeof(nifti_1_header) ) return 1;
3516
3517 if( g_opts.debug > 1 ){
3518 fprintf(stderr,"** bad swapped hsize = %d, unswapped = ", hsize);
3519 nifti_swap_4bytes(1, &hsize); /* swap? */
3520 fprintf(stderr,"%d\n", hsize);
3521 }
3522
3523 return -2; /* bad, naughty hsize */
3524 }
3525
3526
3527 /* use macro LNI_FILE_ERROR instead of ERREX()
3528 #undef ERREX
3529 #define ERREX(msg) \
3530 do{ fprintf(stderr,"** ERROR: nifti_image_read(%s): %s\n", \
3531 (hname != NULL) ? hname : "(null)" , (msg) ) ; \
3532 return NULL ; } while(0)
3533 */
3534
3535
3536 /***************************************************************
3537 * nifti_image_read
3538 ***************************************************************/
3539 /*! \brief Read a nifti header and optionally the data, creating a nifti_image.
3540
3541 - The data buffer will be byteswapped if necessary.
3542 - The data buffer will not be scaled.
3543 - The data buffer is allocated with calloc().
3544
3545 \param hname filename of the nifti dataset
3546 \param read_data Flag, true=read data blob, false=don't read blob.
3547 \return A pointer to the nifti_image data structure.
3548
3549 \sa nifti_image_free, nifti_free_extensions, nifti_image_read_bricks
3550 */
nifti_image_read(const char * hname,int read_data)3551 nifti_image *nifti_image_read( const char *hname , int read_data )
3552 {
3553 struct nifti_1_header nhdr ;
3554 nifti_image *nim ;
3555 znzFile fp ;
3556 int rv, ii , filesize, remaining;
3557 char fname[] = { "nifti_image_read" };
3558 char *hfile=NULL;
3559
3560 if( g_opts.debug > 1 ){
3561 fprintf(stderr,"-d image_read from '%s', read_data = %d",hname,read_data);
3562 #ifdef HAVE_ZLIB
3563 fprintf(stderr,", HAVE_ZLIB = 1\n");
3564 #else
3565 fprintf(stderr,", HAVE_ZLIB = 0\n");
3566 #endif
3567 }
3568
3569 /**- determine filename to use for header */
3570 hfile = nifti_findhdrname(hname);
3571 if( hfile == NULL ){
3572 if(g_opts.debug > 0)
3573 LNI_FERR(fname,"failed to find header file for", hname);
3574 return NULL; /* check return */
3575 } else if( g_opts.debug > 1 )
3576 fprintf(stderr,"-d %s: found header filename '%s'\n",fname,hfile);
3577
3578 if( nifti_is_gzfile(hfile) ) filesize = -1; /* unknown */
3579 else filesize = nifti_get_filesize(hfile);
3580
3581 fp = znzopen(hfile, "rb", nifti_is_gzfile(hfile));
3582 if( znz_isnull(fp) ){
3583 if( g_opts.debug > 0 ) LNI_FERR(fname,"failed to open header file",hfile);
3584 free(hfile);
3585 return NULL;
3586 }
3587
3588 rv = has_ascii_header( fp );
3589 if( rv < 0 ){
3590 if( g_opts.debug > 0 ) LNI_FERR(fname,"short header read",hfile);
3591 znzclose( fp );
3592 free(hfile);
3593 return NULL;
3594 }
3595 else if ( rv == 1 ) /* process special file type */
3596 return nifti_read_ascii_image( fp, hfile, filesize, read_data );
3597
3598 /* else, just process normally */
3599
3600 /**- read binary header */
3601
3602 ii = znzread( &nhdr , 1 , sizeof(nhdr) , fp ) ; /* read the thing */
3603
3604 /* keep file open so we can check for exts. after nifti_convert_nhdr2nim() */
3605
3606 if( ii < sizeof(nhdr) ){
3607 if( g_opts.debug > 0 ){
3608 LNI_FERR(fname,"bad binary header read for file", hfile);
3609 fprintf(stderr," - read %d of %d bytes\n",ii, (int)sizeof(nhdr));
3610 }
3611 znzclose(fp) ;
3612 free(hfile);
3613 return NULL;
3614 }
3615
3616 /* create output image struct and set it up */
3617
3618 /**- convert all nhdr fields to nifti_image fields */
3619 nim = nifti_convert_nhdr2nim(nhdr,hfile);
3620
3621 if( nim == NULL ){
3622 znzclose( fp ) ; /* close the file */
3623 if( g_opts.debug > 0 )
3624 LNI_FERR(fname,"cannot create nifti image from header",hfile);
3625 free(hfile); /* had to save this for debug message */
3626 return NULL;
3627 }
3628
3629 if( g_opts.debug > 3 ){
3630 fprintf(stderr,"+d nifti_image_read(), have nifti image:\n");
3631 if( g_opts.debug > 2 ) nifti_image_infodump(nim);
3632 }
3633
3634 /**- check for extensions (any errors here means no extensions) */
3635 if( NIFTI_ONEFILE(nhdr) ) remaining = nim->iname_offset - sizeof(nhdr);
3636 else remaining = filesize - sizeof(nhdr);
3637
3638 (void)nifti_read_extensions(nim, fp, remaining);
3639
3640 znzclose( fp ) ; /* close the file */
3641 free(hfile);
3642
3643 /**- read the data if desired, then bug out */
3644 if( read_data ){
3645 if( nifti_image_load( nim ) < 0 ){
3646 nifti_image_free(nim); /* take ball, go home. */
3647 return NULL;
3648 }
3649 }
3650 else nim->data = NULL ;
3651
3652 return nim ;
3653 }
3654
3655
3656 /*----------------------------------------------------------------------
3657 * has_ascii_header - see if the NIFTI header is an ASCII format
3658 *
3659 * If the file starts with the ASCII string "<nifti_image", then
3660 * process the dataset as a type-3 .nia file.
3661 *
3662 * return: -1 on error, 1 if true, or 0 if false
3663 *
3664 * NOTE: this is NOT part of the NIFTI-1 standard
3665 *----------------------------------------------------------------------*/
has_ascii_header(znzFile fp)3666 static int has_ascii_header( znzFile fp )
3667 {
3668 char buf[16];
3669 int nread;
3670
3671 if( znz_isnull(fp) ) return 0;
3672
3673 nread = znzread( buf, 1, 12, fp );
3674 buf[12] = '\0';
3675
3676 if( nread < 12 ) return -1;
3677
3678 znzrewind(fp); /* move back to the beginning, and check */
3679
3680 if( strcmp(buf, "<nifti_image") == 0 ) return 1;
3681
3682 return 0;
3683 }
3684
3685
3686 /*----------------------------------------------------------------------*/
3687 /*! nifti_read_ascii_image - process as a type-3 .nia image file
3688
3689 return NULL on failure
3690
3691 NOTE: this is NOT part of the NIFTI-1 standard
3692 *//*--------------------------------------------------------------------*/
nifti_read_ascii_image(znzFile fp,char * fname,int flen,int read_data)3693 nifti_image * nifti_read_ascii_image(znzFile fp, char *fname, int flen,
3694 int read_data)
3695 {
3696 nifti_image * nim;
3697 int slen, txt_size, remain, rv = 0;
3698 char * sbuf, lfunc[25] = { "nifti_read_ascii_image" };
3699
3700 if( nifti_is_gzfile(fname) ){
3701 LNI_FERR(lfunc, "compressed file with negative offset", fname);
3702 free(fname); znzclose(fp); return NULL;
3703 }
3704 slen = flen; /* slen will be our buffer length */
3705
3706 if( g_opts.debug > 1 )
3707 fprintf(stderr,"-d %s: have ASCII NIFTI file of size %d\n",fname,slen);
3708
3709 if( slen > 65530 ) slen = 65530 ;
3710 sbuf = (char *)calloc(sizeof(char),slen+1) ;
3711 if( !sbuf ){
3712 fprintf(stderr,"** %s: failed to alloc %d bytes for sbuf",lfunc,65530);
3713 free(fname); znzclose(fp); return NULL;
3714 }
3715 znzread( sbuf , 1 , slen , fp ) ;
3716 nim = nifti_image_from_ascii( sbuf, &txt_size ) ; free( sbuf ) ;
3717 if( nim == NULL ){
3718 LNI_FERR(lfunc,"failed nifti_image_from_ascii()",fname);
3719 free(fname); znzclose(fp); return NULL;
3720 }
3721 nim->nifti_type = NIFTI_FTYPE_ASCII ;
3722
3723 /* compute remaining space for extensions */
3724 remain = flen - txt_size - (int)nifti_get_volsize(nim);
3725 if( remain > 4 ){
3726 /* read extensions (reposition file pointer, first) */
3727 znzseek(fp, txt_size, SEEK_SET);
3728 (void) nifti_read_extensions(nim, fp, remain);
3729 }
3730
3731 free(fname);
3732 znzclose( fp ) ;
3733
3734 nim->iname_offset = -1 ; /* check from the end of the file */
3735
3736 if( read_data ) rv = nifti_image_load( nim ) ;
3737 else nim->data = NULL ;
3738
3739 /* check for nifti_image_load() failure, maybe bail out */
3740 if( read_data && rv != 0 ){
3741 if( g_opts.debug > 1 )
3742 fprintf(stderr,"-d failed image_load, free nifti image struct\n");
3743 free(nim);
3744 return NULL;
3745 }
3746
3747 return nim ;
3748 }
3749
3750
3751 /*----------------------------------------------------------------------
3752 * Read the extensions into the nifti_image struct 08 Dec 2004 [rickr]
3753 *
3754 * This function is called just after the header struct is read in, and
3755 * it is assumed the file pointer has not moved. The value in remain
3756 * is assumed to be accurate, reflecting the bytes of space for potential
3757 * extensions.
3758 *
3759 * return the number of extensions read in, or < 0 on error
3760 *----------------------------------------------------------------------*/
nifti_read_extensions(nifti_image * nim,znzFile fp,int remain)3761 static int nifti_read_extensions( nifti_image *nim, znzFile fp, int remain )
3762 {
3763 nifti1_extender extdr; /* defines extension existence */
3764 nifti1_extension extn; /* single extension to process */
3765 nifti1_extension * Elist; /* list of processed extensions */
3766 int posn, count;
3767
3768 if( !nim || znz_isnull(fp) ) {
3769 if( g_opts.debug > 0 )
3770 fprintf(stderr,"** nifti_read_extensions: bad inputs (%p,%p)\n",
3771 (void *)nim, (void *)fp);
3772 return -1;
3773 }
3774
3775 posn = znztell(fp);
3776
3777 if( (posn != sizeof(nifti_1_header)) &&
3778 (nim->nifti_type != NIFTI_FTYPE_ASCII) )
3779 fprintf(stderr,"** WARNING: posn not header size (%d, %d)\n",
3780 posn, (int)sizeof(nifti_1_header));
3781
3782 if( g_opts.debug > 2 )
3783 fprintf(stderr,"-d nre: posn = %d, offset = %d, type = %d, remain = %d\n",
3784 posn, nim->iname_offset, nim->nifti_type, remain);
3785
3786 if( remain < 16 ){
3787 if( g_opts.debug > 2 ){
3788 if( g_opts.skip_blank_ext )
3789 fprintf(stderr,"-d no extender in '%s' is okay, as "
3790 "skip_blank_ext is set\n",nim->fname);
3791 else
3792 fprintf(stderr,"-d no space for extensions\n");
3793 }
3794 return 0;
3795 }
3796
3797 count = znzread( extdr.extension, 1, 4, fp ); /* get extender */
3798
3799 if( count < 4 ){
3800 if( g_opts.debug > 1 )
3801 fprintf(stderr,"-d file '%s' is too short for an extender\n",
3802 nim->fname);
3803 return 0;
3804 }
3805
3806 if( extdr.extension[0] != 1 ){
3807 if( g_opts.debug > 2 )
3808 fprintf(stderr,"-d extender[0] (%d) shows no extensions for '%s'\n",
3809 extdr.extension[0], nim->fname);
3810 return 0;
3811 }
3812
3813 remain -= 4;
3814 if( g_opts.debug > 2 )
3815 fprintf(stderr,"-d found valid 4-byte extender, remain = %d\n", remain);
3816
3817 /* so we expect extensions, but have no idea of how many there may be */
3818
3819 count = 0;
3820 Elist = NULL;
3821 while (nifti_read_next_extension(&extn, nim, remain, fp) > 0)
3822 {
3823 if( nifti_add_exten_to_list(&extn, &Elist, count+1) < 0 ){
3824 if( g_opts.debug > 0 )
3825 fprintf(stderr,"** failed adding ext %d to list\n", count);
3826 return -1;
3827 }
3828
3829 /* we have a new extension */
3830 if( g_opts.debug > 1 ){
3831 fprintf(stderr,"+d found extension #%d, code = 0x%x, size = %d\n",
3832 count, extn.ecode, extn.esize);
3833 if( extn.ecode == NIFTI_ECODE_AFNI && g_opts.debug > 2 ) /* ~XML */
3834 fprintf(stderr," AFNI extension: %.*s\n",
3835 extn.esize-8,extn.edata);
3836 else if( extn.ecode == NIFTI_ECODE_COMMENT && g_opts.debug > 2 )
3837 fprintf(stderr," COMMENT extension: %.*s\n", /* TEXT */
3838 extn.esize-8,extn.edata);
3839 }
3840 remain -= extn.esize;
3841 count++;
3842 }
3843
3844 if( g_opts.debug > 2 ) fprintf(stderr,"+d found %d extension(s)\n", count);
3845
3846 nim->num_ext = count;
3847 nim->ext_list = Elist;
3848
3849 return count;
3850 }
3851
3852
3853 /*----------------------------------------------------------------------*/
3854 /*! nifti_add_extension - add an extension, with a copy of the data
3855
3856 Add an extension to the nim->ext_list array.
3857 Fill this extension with a copy of the data, noting the
3858 length and extension code.
3859
3860 \param nim - nifti_image to add extension to
3861 \param data - raw extension data
3862 \param length - length of raw extension data
3863 \param ecode - extension code
3864
3865 \sa extension codes NIFTI_ECODE_* in nifti1_io.h
3866 \sa nifti_free_extensions, valid_nifti_extensions, nifti_copy_extensions
3867
3868 \return 0 on success, -1 on error (and free the entire list)
3869 *//*--------------------------------------------------------------------*/
nifti_add_extension(nifti_image * nim,const char * data,int len,int ecode)3870 int nifti_add_extension(nifti_image *nim, const char * data, int len, int ecode)
3871 {
3872 nifti1_extension ext;
3873
3874 /* error are printed in functions */
3875 if( nifti_fill_extension(&ext, data, len, ecode) ) return -1;
3876 if( nifti_add_exten_to_list(&ext, &nim->ext_list, nim->num_ext+1)) return -1;
3877
3878 nim->num_ext++; /* success, so increment */
3879
3880 return 0;
3881 }
3882
3883
3884 /*----------------------------------------------------------------------*/
3885 /* nifti_add_exten_to_list - add a new nifti1_extension to the list
3886
3887 We will append via "malloc, copy and free", because on an error,
3888 the old data pointers must all be released (sorry realloc(), only
3889 quality dolphins get to become part of St@rk!st brand tunafish).
3890
3891 return 0 on success, -1 on error (and free the entire list)
3892 *//*--------------------------------------------------------------------*/
nifti_add_exten_to_list(nifti1_extension * new_ext,nifti1_extension ** list,int new_length)3893 static int nifti_add_exten_to_list( nifti1_extension * new_ext,
3894 nifti1_extension ** list, int new_length )
3895 {
3896 nifti1_extension * tmplist;
3897 int count;
3898
3899 tmplist = *list;
3900 *list = (nifti1_extension *)malloc(new_length * sizeof(nifti1_extension));
3901
3902 /* check for failure first */
3903 if( ! *list ){
3904 fprintf(stderr,"** failed to alloc %d extension structs (%d bytes)\n",
3905 new_length, new_length*(int)sizeof(nifti1_extension));
3906 if( !tmplist ) return -1; /* no old list to lose */
3907
3908 for ( count = 0; count < new_length-1; count++ )
3909 if( tmplist[count].edata ) free(tmplist[count].edata);
3910 free(tmplist);
3911 return -1;
3912 }
3913
3914 /* we have memory, so copy the old and insert the new */
3915 memcpy(*list, tmplist, (new_length-1)*sizeof(nifti1_extension));
3916
3917 /* for some reason, I just don't like struct copy... */
3918 (*list)[new_length-1].esize = new_ext->esize;
3919 (*list)[new_length-1].ecode = new_ext->ecode;
3920 (*list)[new_length-1].edata = new_ext->edata;
3921
3922 if( g_opts.debug > 2 )
3923 fprintf(stderr,"+d allocated and appended extension #%d to list\n",
3924 new_length);
3925
3926 return 0;
3927 }
3928
3929
3930 /*----------------------------------------------------------------------*/
3931 /* nifti_fill_extension - given data and length, fill an extension struct
3932
3933 Allocate memory for data, copy data, set the size and code.
3934
3935 return 0 on success, -1 on error (and free the entire list)
3936 *//*--------------------------------------------------------------------*/
nifti_fill_extension(nifti1_extension * ext,const char * data,int len,int ecode)3937 static int nifti_fill_extension( nifti1_extension *ext, const char * data,
3938 int len, int ecode)
3939 {
3940 int esize;
3941
3942 if( !ext || !data || len < 0 ){
3943 fprintf(stderr,"** fill_ext: bad params (%p,%p,%d)\n",
3944 (void *)ext, data, len);
3945 return -1;
3946 } else if( ! nifti_is_valid_ecode(ecode) ){
3947 fprintf(stderr,"** fill_ext: invalid ecode %d\n", ecode);
3948 return -1;
3949 }
3950
3951 /* compute esize, first : len+8, and take ceiling up to a mult of 16 */
3952 esize = len+8;
3953 if( esize & 0xf ) esize = (esize + 0xf) & ~0xf;
3954 ext->esize = esize;
3955
3956 /* allocate esize-8 (maybe more than len), using calloc for fill */
3957 ext->edata = (char *)calloc(esize-8, sizeof(char));
3958 if( !ext->edata ){
3959 fprintf(stderr,"** NFE: failed to alloc %d bytes for extension\n",len);
3960 return -1;
3961 }
3962
3963 memcpy(ext->edata, data, len); /* copy the data, using len */
3964 ext->ecode = ecode; /* set the ecode */
3965
3966 if( g_opts.debug > 2 )
3967 fprintf(stderr,"+d alloc %d bytes for ext len %d, ecode %d, esize %d\n",
3968 esize-8, len, ecode, esize);
3969
3970 return 0;
3971 }
3972
3973
3974 /*----------------------------------------------------------------------
3975 * nifti_read_next_extension - read a single extension from the file
3976 *
3977 * return (>= 0 is okay):
3978 *
3979 * success : esize
3980 * no extension : 0
3981 * error : -1
3982 *----------------------------------------------------------------------*/
nifti_read_next_extension(nifti1_extension * nex,nifti_image * nim,int remain,znzFile fp)3983 static int nifti_read_next_extension( nifti1_extension * nex, nifti_image *nim,
3984 int remain, znzFile fp )
3985 {
3986 int swap = nim->byteorder != nifti_short_order();
3987 int count, size, code;
3988
3989 /* first clear nex */
3990 nex->esize = nex->ecode = 0;
3991 nex->edata = NULL;
3992
3993 if( remain < 16 ){
3994 if( g_opts.debug > 2 )
3995 fprintf(stderr,"-d only %d bytes remain, so no extension\n", remain);
3996 return 0;
3997 }
3998
3999 /* must start with 4-byte size and code */
4000 count = znzread( &size, 4, 1, fp );
4001 if( count == 1 ) count += znzread( &code, 4, 1, fp );
4002
4003 if( count != 2 ){
4004 if( g_opts.debug > 2 )
4005 fprintf(stderr,"-d current extension read failed\n");
4006 znzseek(fp, -4*count, SEEK_CUR); /* back up past any read */
4007 return 0; /* no extension, no error condition */
4008 }
4009
4010 if( swap ){
4011 if( g_opts.debug > 2 )
4012 fprintf(stderr,"-d pre-swap exts: code %d, size %d\n", code, size);
4013
4014 nifti_swap_4bytes(1, &size);
4015 nifti_swap_4bytes(1, &code);
4016 }
4017
4018 if( g_opts.debug > 2 )
4019 fprintf(stderr,"-d potential extension: code %d, size %d\n", code, size);
4020
4021 if( !nifti_check_extension(nim, size, code, remain) ){
4022 if( znzseek(fp, -8, SEEK_CUR) < 0 ){ /* back up past any read */
4023 fprintf(stderr,"** failure to back out of extension read!\n");
4024 return -1;
4025 }
4026 return 0;
4027 }
4028
4029 /* now get the actual data */
4030 nex->esize = size;
4031 nex->ecode = code;
4032
4033 size -= 8; /* subtract space for size and code in extension */
4034 nex->edata = (char *)malloc(size * sizeof(char));
4035 if( !nex->edata ){
4036 fprintf(stderr,"** failed to allocate %d bytes for extension\n",size);
4037 return -1;
4038 }
4039
4040 count = znzread(nex->edata, 1, size, fp);
4041 if( count < size ){
4042 if( g_opts.debug > 0 )
4043 fprintf(stderr,"-d read only %d (of %d) bytes for extension\n",
4044 count, size);
4045 free(nex->edata);
4046 nex->edata = NULL;
4047 return -1;
4048 }
4049
4050 /* success! */
4051 if( g_opts.debug > 2 )
4052 fprintf(stderr,"+d successfully read extension, code %d, size %d\n",
4053 nex->ecode, nex->esize);
4054
4055 return nex->esize;
4056 }
4057
4058
4059 /*----------------------------------------------------------------------*/
4060 /*! for each extension, check code, size and data pointer
4061 *//*--------------------------------------------------------------------*/
valid_nifti_extensions(const nifti_image * nim)4062 int valid_nifti_extensions(const nifti_image * nim)
4063 {
4064 nifti1_extension * ext;
4065 int c, errs;
4066
4067 if( nim->num_ext <= 0 || nim->ext_list == NULL ){
4068 if( g_opts.debug > 2 ) fprintf(stderr,"-d empty extension list\n");
4069 return 0;
4070 }
4071
4072 /* for each extension, check code, size and data pointer */
4073 ext = nim->ext_list;
4074 errs = 0;
4075 for ( c = 0; c < nim->num_ext; c++ ){
4076 if( ! nifti_is_valid_ecode(ext->ecode) ) {
4077 if( g_opts.debug > 1 )
4078 fprintf(stderr,"-d ext %d, invalid code %d\n", c, ext->ecode);
4079 errs++;
4080 }
4081
4082 if( ext->esize <= 0 ){
4083 if( g_opts.debug > 1 )
4084 fprintf(stderr,"-d ext %d, bad size = %d\n", c, ext->esize);
4085 errs++;
4086 } else if( ext->esize & 0xf ){
4087 if( g_opts.debug > 1 )
4088 fprintf(stderr,"-d ext %d, size %d not multiple of 16\n",
4089 c, ext->esize);
4090 errs++;
4091 }
4092
4093 if( ext->edata == NULL ){
4094 if( g_opts.debug > 1 ) fprintf(stderr,"-d ext %d, missing data\n", c);
4095 errs++;
4096 }
4097
4098 ext++;
4099 }
4100
4101 if( errs > 0 ){
4102 if( g_opts.debug > 0 )
4103 fprintf(stderr,"-d had %d extension errors, none will be written\n",
4104 errs);
4105 return 0;
4106 }
4107
4108 /* if we're here, we're good */
4109 return 1;
4110 }
4111
4112
4113 /*----------------------------------------------------------------------*/
4114 /*! check whether the extension code is valid
4115
4116 \return 1 if valid, 0 otherwise
4117 *//*--------------------------------------------------------------------*/
nifti_is_valid_ecode(int ecode)4118 int nifti_is_valid_ecode( int ecode )
4119 {
4120 if( ecode < NIFTI_ECODE_IGNORE || /* minimum code number (0) */
4121 ecode > NIFTI_MAX_ECODE || /* maximum code number */
4122 ecode & 1 ) /* cannot be odd */
4123 return 0;
4124
4125 return 1;
4126 }
4127
4128
4129 /*----------------------------------------------------------------------
4130 * check for valid size and code, as well as can be done
4131 *----------------------------------------------------------------------*/
nifti_check_extension(nifti_image * nim,int size,int code,int rem)4132 static int nifti_check_extension(nifti_image *nim, int size, int code, int rem)
4133 {
4134 /* check for bad code before bad size */
4135 if( ! nifti_is_valid_ecode(code) ) {
4136 if( g_opts.debug > 2 )
4137 fprintf(stderr,"-d invalid extension code %d\n",code);
4138 return 0;
4139 }
4140
4141 if( size < 16 ){
4142 if( g_opts.debug > 2 )
4143 fprintf(stderr,"-d ext size %d, no extension\n",size);
4144 return 0;
4145 }
4146
4147 if( size > rem ){
4148 if( g_opts.debug > 2 )
4149 fprintf(stderr,"-d ext size %d, space %d, no extension\n", size, rem);
4150 return 0;
4151 }
4152
4153 if( size & 0xf ){
4154 if( g_opts.debug > 2 )
4155 fprintf(stderr,"-d nifti extension size %d not multiple of 16\n",size);
4156 return 0;
4157 }
4158
4159 if( nim->nifti_type == NIFTI_FTYPE_ASCII && size > LNI_MAX_NIA_EXT_LEN ){
4160 if( g_opts.debug > 2 )
4161 fprintf(stderr,"-d NVE, bad nifti_type 3 size %d\n", size);
4162 return 0;
4163 }
4164
4165 return 1;
4166 }
4167
4168
4169 /*----------------------------------------------------------------------
4170 * nifti_image_load_prep - prepare to read data
4171 *
4172 * Check nifti_image fields, open the file and seek to the appropriate
4173 * offset for reading.
4174 *
4175 * return NULL on failure
4176 *----------------------------------------------------------------------*/
nifti_image_load_prep(nifti_image * nim)4177 static znzFile nifti_image_load_prep( nifti_image *nim )
4178 {
4179 /* set up data space, open data file and seek, then call nifti_read_buffer */
4180 size_t ntot , ii , ioff;
4181 znzFile fp;
4182 char fname[] = { "nifti_image_load_prep" };
4183
4184 /**- perform sanity checks */
4185 if( nim == NULL || nim->iname == NULL ||
4186 nim->nbyper <= 0 || nim->nvox <= 0 )
4187 {
4188 if ( g_opts.debug > 0 ){
4189 if( !nim ) fprintf(stderr,"** ERROR: N_image_load: no nifti image\n");
4190 else fprintf(stderr,"** ERROR: N_image_load: bad params (%p,%d,%d)\n",
4191 nim->iname, nim->nbyper, nim->nvox);
4192 }
4193 return NULL;
4194 }
4195
4196 ntot = nifti_get_volsize(nim) ; /* total bytes to read */
4197
4198 /**- open image data file */
4199
4200 fp = znzopen(nim->iname, "rb", nifti_is_gzfile(nim->iname));
4201 if( znz_isnull(fp) ){
4202 if( g_opts.debug > 0 ) LNI_FERR(fname,"cannot open data file",nim->iname);
4203 return NULL;
4204 }
4205
4206 /**- get image offset: a negative offset means to figure from end of file */
4207 if( nim->iname_offset < 0 ){
4208 if( nifti_is_gzfile(nim->iname) ){
4209 if( g_opts.debug > 0 )
4210 LNI_FERR(fname,"negative offset for compressed file",nim->iname);
4211 znzclose(fp);
4212 return NULL;
4213 }
4214 ii = nifti_get_filesize( nim->iname ) ;
4215 if( ii <= 0 ){
4216 if( g_opts.debug > 0 ) LNI_FERR(fname,"empty data file",nim->iname);
4217 znzclose(fp);
4218 return NULL;
4219 }
4220 ioff = (ii > ntot) ? ii-ntot : 0 ;
4221 } else { /* non-negative offset */
4222 ioff = nim->iname_offset ; /* means use it directly */
4223 }
4224
4225 /**- seek to the appropriate read position */
4226 if( znzseek(fp , ioff , SEEK_SET) < 0 ){
4227 fprintf(stderr,"** could not seek to offset %d in file '%s'\n",
4228 (int)ioff, nim->iname);
4229 znzclose(fp);
4230 return NULL;
4231 }
4232
4233 /**- and return the File pointer */
4234 return fp;
4235 }
4236
4237
4238 /*----------------------------------------------------------------------
4239 * nifti_image_load
4240 *----------------------------------------------------------------------*/
4241 /*! \fn int nifti_image_load( nifti_image *nim )
4242 \brief Load the image blob into a previously initialized nifti_image.
4243
4244 - If not yet set, the data buffer is allocated with calloc().
4245 - The data buffer will be byteswapped if necessary.
4246 - The data buffer will not be scaled.
4247
4248 This function is used to read the image from disk. It should be used
4249 after a function such as nifti_image_read(), so that the nifti_image
4250 structure is already initialized.
4251
4252 \param nim pointer to a nifti_image (previously initialized)
4253 \return 0 on success, -1 on failure
4254 \sa nifti_image_read, nifti_image_free, nifti_image_unload
4255 */
nifti_image_load(nifti_image * nim)4256 int nifti_image_load( nifti_image *nim )
4257 {
4258 /* set up data space, open data file and seek, then call nifti_read_buffer */
4259 size_t ntot , ii ;
4260 znzFile fp ;
4261
4262 /**- open the file and position the FILE pointer */
4263 fp = nifti_image_load_prep( nim );
4264
4265 if( fp == NULL ){
4266 if( g_opts.debug > 0 )
4267 fprintf(stderr,"** nifti_image_load, failed load_prep\n");
4268 return -1;
4269 }
4270
4271 ntot = nifti_get_volsize(nim);
4272
4273 /**- if the data pointer is not yet set, get memory space for the image */
4274
4275 if( nim->data == NULL )
4276 {
4277 nim->data = (void *)calloc(1,ntot) ; /* create image memory */
4278 if( nim->data == NULL ){
4279 if( g_opts.debug > 0 )
4280 fprintf(stderr,"** failed to alloc %d bytes for image data\n",
4281 (int)ntot);
4282 znzclose(fp);
4283 return -1;
4284 }
4285 }
4286
4287 /**- now that everything is set up, do the reading */
4288 ii = nifti_read_buffer(fp,nim->data,ntot,nim);
4289 if( ii < ntot ){
4290 znzclose(fp) ;
4291 free(nim->data) ;
4292 nim->data = NULL ;
4293 return -1 ; /* errors were printed in nifti_read_buffer() */
4294 }
4295
4296 /**- close the file */
4297 znzclose( fp ) ;
4298
4299 return 0 ;
4300 }
4301
4302
4303 /* 30 Nov 2004 [rickr]
4304 #undef ERREX
4305 #define ERREX(msg) \
4306 do{ fprintf(stderr,"** ERROR: nifti_read_buffer: %s\n",(msg)) ; \
4307 return 0; } while(0)
4308 */
4309
4310 /*----------------------------------------------------------------------*/
4311 /*! read ntot bytes of data from an open file and byte swaps if necessary
4312
4313 note that nifti_image is required for information on datatype, bsize
4314 (for any needed byte swapping), etc.
4315
4316 This function does not allocate memory, so dataptr must be valid.
4317 *//*--------------------------------------------------------------------*/
nifti_read_buffer(znzFile fp,void * dataptr,size_t ntot,nifti_image * nim)4318 size_t nifti_read_buffer(znzFile fp, void* dataptr, size_t ntot,
4319 nifti_image *nim)
4320 {
4321 size_t ii;
4322
4323 if( dataptr == NULL ){
4324 if( g_opts.debug > 0 )
4325 fprintf(stderr,"** ERROR: nifti_read_buffer: NULL dataptr\n");
4326 return -1;
4327 }
4328
4329 ii = znzread( dataptr , 1 , ntot , fp ) ; /* data input */
4330
4331 /* if read was short, fail */
4332 if( ii < ntot ){
4333 if( g_opts.debug > 0 )
4334 fprintf(stderr,"++ WARNING: nifti_read_buffer(%s):\n"
4335 " data bytes needed = %u\n"
4336 " data bytes input = %u\n"
4337 " number missing = %u (set to 0)\n",
4338 nim->iname , (unsigned int)ntot ,
4339 (unsigned int)ii , (unsigned int)(ntot-ii) ) ;
4340 /* memset( (char *)(dataptr)+ii , 0 , ntot-ii ) ; now failure [rickr] */
4341 return -1 ;
4342 }
4343
4344 /* byte swap array if needed */
4345
4346 if( nim->swapsize > 1 && nim->byteorder != nifti_short_order() )
4347 nifti_swap_Nbytes( ntot / nim->swapsize , nim->swapsize , dataptr ) ;
4348
4349 #ifdef isfinite
4350 {
4351 /* check input float arrays for goodness, and fix bad floats */
4352 int fix_count = 0 ;
4353
4354 switch( nim->datatype ){
4355
4356 case NIFTI_TYPE_FLOAT32:
4357 case NIFTI_TYPE_COMPLEX64:{
4358 register float *far = (float *)dataptr ; register int jj,nj ;
4359 nj = ntot / sizeof(float) ;
4360 for( jj=0 ; jj < nj ; jj++ ) /* count fixes 30 Nov 2004 [rickr] */
4361 if( !IS_GOOD_FLOAT(far[jj]) ){
4362 far[jj] = 0 ;
4363 fix_count++ ;
4364 }
4365 }
4366 break ;
4367
4368 case NIFTI_TYPE_FLOAT64:
4369 case NIFTI_TYPE_COMPLEX128:{
4370 register double *far = (double *)dataptr ; register int jj,nj ;
4371 nj = ntot / sizeof(double) ;
4372 for( jj=0 ; jj < nj ; jj++ ) /* count fixes 30 Nov 2004 [rickr] */
4373 if( !IS_GOOD_FLOAT(far[jj]) ){
4374 far[jj] = 0 ;
4375 fix_count++ ;
4376 }
4377 }
4378 break ;
4379
4380 }
4381
4382 if( g_opts.debug > 1 )
4383 fprintf(stderr,"+d in image, %d bad floats were set to 0\n", fix_count);
4384 }
4385 #endif
4386
4387 return ii;
4388 }
4389
4390 /*--------------------------------------------------------------------------*/
4391 /*! Unload the data in a nifti_image struct, but keep the metadata.
4392 *//*------------------------------------------------------------------------*/
nifti_image_unload(nifti_image * nim)4393 void nifti_image_unload( nifti_image *nim )
4394 {
4395 if( nim != NULL && nim->data != NULL ){
4396 free(nim->data) ; nim->data = NULL ;
4397 }
4398 return ;
4399 }
4400
4401 /*--------------------------------------------------------------------------*/
4402 /*! free 'everything' about a nifti_image struct (including the passed struct)
4403
4404 free (only fields which are not NULL):
4405 - fname and iname
4406 - data
4407 - any ext_list[i].edata
4408 - ext_list
4409 - nim
4410 *//*------------------------------------------------------------------------*/
nifti_image_free(nifti_image * nim)4411 void nifti_image_free( nifti_image *nim )
4412 {
4413 if( nim == NULL ) return ;
4414 if( nim->fname != NULL ) free(nim->fname) ;
4415 if( nim->iname != NULL ) free(nim->iname) ;
4416 if( nim->data != NULL ) free(nim->data ) ;
4417 (void)nifti_free_extensions( nim ) ;
4418 free(nim) ; return ;
4419 }
4420
4421
4422 /*--------------------------------------------------------------------------*/
4423 /*! free the nifti extensions
4424
4425 - If any edata pointer is set in the extension list, free() it.
4426 - Free ext_list, if it is set.
4427 - Clear num_ext and ext_list from nim.
4428
4429 \return 0 on success, -1 on error
4430
4431 \sa nifti_add_extension, nifti_copy_extensions
4432 *//*------------------------------------------------------------------------*/
nifti_free_extensions(nifti_image * nim)4433 int nifti_free_extensions( nifti_image *nim )
4434 {
4435 int c ;
4436 if( nim == NULL ) return -1;
4437 if( nim->num_ext > 0 && nim->ext_list ){
4438 for( c = 0; c < nim->num_ext; c++ )
4439 if ( nim->ext_list[c].edata ) free(nim->ext_list[c].edata);
4440 free(nim->ext_list);
4441 }
4442 /* or if it is inconsistent, warn the user (if we are not in quiet mode) */
4443 else if ( (nim->num_ext > 0 || nim->ext_list != NULL) && (g_opts.debug > 0) )
4444 fprintf(stderr,"** warning: nifti extension num/ptr mismatch (%d,%p)\n",
4445 nim->num_ext, (void *)nim->ext_list);
4446
4447 if( g_opts.debug > 2 )
4448 fprintf(stderr,"+d free'd %d extension(s)\n", nim->num_ext);
4449
4450 nim->num_ext = 0;
4451 nim->ext_list = NULL;
4452
4453 return 0;
4454 }
4455
4456
4457 /*--------------------------------------------------------------------------*/
4458 /*! Print to stdout some info about a nifti_image struct.
4459 *//*------------------------------------------------------------------------*/
nifti_image_infodump(const nifti_image * nim)4460 void nifti_image_infodump( const nifti_image *nim )
4461 {
4462 char *str = nifti_image_to_ascii( nim ) ;
4463 /* stdout -> stderr 2 Dec 2004 [rickr] */
4464 if( str != NULL ){ fputs(str,stderr) ; free(str) ; }
4465 return ;
4466 }
4467
4468
4469 /*--------------------------------------------------------------------------
4470 * nifti_write_buffer just check for a null znzFile and call znzwrite
4471 *--------------------------------------------------------------------------*/
4472 /*! \fn size_t nifti_write_buffer(znzFile fp, void *buffer, size_t numbytes)
4473 \brief write numbytes of buffer to file, fp
4474
4475 \param fp File pointer (from znzopen) to gzippable nifti datafile
4476 \param buffer data buffer to be written
4477 \param numbytes number of bytes in buffer to write
4478 \return number of bytes successfully written
4479 */
nifti_write_buffer(znzFile fp,const void * buffer,size_t numbytes)4480 size_t nifti_write_buffer(znzFile fp, const void *buffer, size_t numbytes)
4481 {
4482 /* Write all the image data at once (no swapping here) */
4483 size_t ss;
4484 if (znz_isnull(fp)){
4485 fprintf(stderr,"** ERROR: nifti_write_buffer: null file pointer\n");
4486 return 0;
4487 }
4488 ss = znzwrite( (void*)buffer , 1 , numbytes , fp ) ;
4489 return ss;
4490 }
4491
4492
4493 /*----------------------------------------------------------------------*/
4494 /*! write the nifti_image data to file (from nim->data or from NBL)
4495
4496 If NBL is not NULL, write the data from that structure. Otherwise,
4497 write it out from nim->data. No swapping is done here.
4498
4499 \param fp : File pointer
4500 \param nim : nifti_image corresponding to the data
4501 \param NBL : optional source of write data (if NULL use nim->data)
4502
4503 \return 0 on success, -1 on failure
4504
4505 Note: the nifti_image byte_order is set as that of the current CPU.
4506 This is because such a conversion was made to the data upon
4507 reading, while byte_order was not set (so the programs would
4508 know what format the data was on disk). Effectively, since
4509 byte_order should match what is on disk, it should bet set to
4510 that of the current CPU whenever new filenames are assigned.
4511 *//*--------------------------------------------------------------------*/
nifti_write_all_data(znzFile fp,nifti_image * nim,const nifti_brick_list * NBL)4512 int nifti_write_all_data(znzFile fp, nifti_image * nim,
4513 const nifti_brick_list * NBL)
4514 {
4515 size_t ss;
4516 int bnum;
4517
4518 if( !NBL ){ /* just write one buffer and get out of here */
4519 if( nim->data == NULL ){
4520 fprintf(stderr,"** NWAD: no image data to write\n");
4521 return -1;
4522 }
4523
4524 ss = nifti_write_buffer(fp,nim->data,nim->nbyper * nim->nvox);
4525 if (ss < (size_t)(nim->nbyper * nim->nvox)){
4526 fprintf(stderr,
4527 "** ERROR: NWAD: wrote only %d of %d bytes to file\n",
4528 (int)ss, nim->nbyper * nim->nvox);
4529 return -1;
4530 }
4531
4532 if( g_opts.debug > 1 )
4533 fprintf(stderr,"+d wrote single image of %d bytes\n",(int)ss);
4534 } else {
4535 if( ! NBL->bricks || NBL->nbricks <= 0 || NBL->bsize <= 0 ){
4536 fprintf(stderr,"** NWAD: no brick data to write (%p,%d,%d)\n",
4537 (void *)NBL->bricks, NBL->nbricks, NBL->bsize);
4538 return -1;
4539 }
4540
4541 for( bnum = 0; bnum < NBL->nbricks; bnum++ ){
4542 ss = nifti_write_buffer(fp, NBL->bricks[bnum], NBL->bsize);
4543 if( ss < (size_t)NBL->bsize ){
4544 fprintf(stderr,
4545 "** NWAD ERROR: wrote %d of %d bytes of brick %d of %d to file",
4546 (int)ss, NBL->bsize, bnum+1, NBL->nbricks);
4547 return -1;
4548 }
4549 }
4550 if( g_opts.debug > 1 )
4551 fprintf(stderr,"+d wrote image of %d brick(s), each of %d bytes\n",
4552 NBL->nbricks, NBL->bsize);
4553 }
4554
4555 /* mark as being in this CPU byte order */
4556 nim->byteorder = nifti_short_order() ;
4557
4558 return 0;
4559 }
4560
4561 /* return number of extensions written, or -1 on error */
nifti_write_extensions(znzFile fp,nifti_image * nim)4562 static int nifti_write_extensions(znzFile fp, nifti_image *nim)
4563 {
4564 nifti1_extension * list;
4565 char extdr[4] = { 0, 0, 0, 0 };
4566 int c, size, ok = 1;
4567
4568 if( znz_isnull(fp) || !nim || nim->num_ext < 0 ){
4569 if( g_opts.debug > 0 )
4570 fprintf(stderr,"** nifti_write_extensions, bad params\n");
4571 return -1;
4572 }
4573
4574 /* if no extensions and user requests it, skip extender */
4575 if( g_opts.skip_blank_ext && (nim->num_ext == 0 || ! nim->ext_list ) ){
4576 if( g_opts.debug > 1 )
4577 fprintf(stderr,"-d no exts and skip_blank_ext set, "
4578 "so skipping 4-byte extender\n");
4579 return 0;
4580 }
4581
4582 /* if invalid extension list, clear num_ext */
4583 if( ! valid_nifti_extensions(nim) ) nim->num_ext = 0;
4584
4585 /* write out extender block */
4586 if( nim->num_ext > 0 ) extdr[0] = 1;
4587 if( nifti_write_buffer(fp, extdr, 4) != 4 ){
4588 fprintf(stderr,"** failed to write extender\n");
4589 return -1;
4590 }
4591
4592 list = nim->ext_list;
4593 for ( c = 0; c < nim->num_ext; c++ ){
4594 size = nifti_write_buffer(fp, &list->esize, sizeof(int));
4595 ok = (size == (int)sizeof(int));
4596 if( ok ){
4597 size = nifti_write_buffer(fp, &list->ecode, sizeof(int));
4598 ok = (size == (int)sizeof(int));
4599 }
4600 if( ok ){
4601 size = nifti_write_buffer(fp, list->edata, list->esize - 8);
4602 ok = (size == list->esize - 8);
4603 }
4604
4605 if( !ok ){
4606 fprintf(stderr,"** failed while writing extension #%d\n",c);
4607 return -1;
4608 }
4609
4610 list++;
4611 }
4612
4613 if( g_opts.debug > 1 )
4614 fprintf(stderr,"+d wrote out %d extension(s)\n", nim->num_ext);
4615
4616 return nim->num_ext;
4617 }
4618
4619
4620 /*----------------------------------------------------------------------*/
4621 /*! basic initialization of a nifti_image struct (to a 1x1x1 image)
4622 *//*--------------------------------------------------------------------*/
nifti_simple_init_nim(void)4623 nifti_image* nifti_simple_init_nim(void)
4624 {
4625 nifti_image *nim;
4626 struct nifti_1_header nhdr;
4627 int nbyper, swapsize;
4628
4629 memset(&nhdr,0,sizeof(nhdr)) ; /* zero out header, to be safe */
4630
4631 nhdr.sizeof_hdr = sizeof(nhdr) ;
4632 nhdr.regular = 'r' ; /* for some stupid reason */
4633
4634 nhdr.dim[0] = 3 ;
4635 nhdr.dim[1] = 1 ; nhdr.dim[2] = 1 ; nhdr.dim[3] = 1 ;
4636 nhdr.dim[4] = 0 ;
4637
4638 nhdr.pixdim[0] = 0.0 ;
4639 nhdr.pixdim[1] = 1.0 ; nhdr.pixdim[2] = 1.0 ;
4640 nhdr.pixdim[3] = 1.0 ;
4641
4642 nhdr.datatype = DT_FLOAT32 ;
4643 nifti_datatype_sizes( nhdr.datatype , &nbyper, &swapsize );
4644 nhdr.bitpix = 8 * nbyper ;
4645
4646 nim = nifti_convert_nhdr2nim(nhdr,NULL);
4647 nim->fname = NULL;
4648 nim->iname = NULL;
4649 return nim;
4650 }
4651
4652
4653 /*----------------------------------------------------------------------*/
4654 /*! convert a nifti_image structure to a nifti_1_header struct
4655
4656 No allocation is done, this should be used via structure copy.
4657 As in:
4658 <pre>
4659 nifti_1_header my_header;
4660 my_header = nifti_convert_nim2nhdr(my_nim_pointer);
4661 </pre>
4662 *//*--------------------------------------------------------------------*/
nifti_convert_nim2nhdr(const nifti_image * nim)4663 struct nifti_1_header nifti_convert_nim2nhdr(const nifti_image * nim)
4664 {
4665 struct nifti_1_header nhdr;
4666
4667 memset(&nhdr,0,sizeof(nhdr)) ; /* zero out header, to be safe */
4668
4669
4670 /**- load the ANALYZE-7.5 generic parts of the header struct */
4671
4672 nhdr.sizeof_hdr = sizeof(nhdr) ;
4673 nhdr.regular = 'r' ; /* for some stupid reason */
4674
4675 nhdr.dim[0] = nim->ndim ;
4676 nhdr.dim[1] = nim->nx ; nhdr.dim[2] = nim->ny ; nhdr.dim[3] = nim->nz ;
4677 nhdr.dim[4] = nim->nt ; nhdr.dim[5] = nim->nu ; nhdr.dim[6] = nim->nv ;
4678 nhdr.dim[7] = nim->nw ;
4679
4680 nhdr.pixdim[0] = 0.0 ;
4681 nhdr.pixdim[1] = nim->dx ; nhdr.pixdim[2] = nim->dy ;
4682 nhdr.pixdim[3] = nim->dz ; nhdr.pixdim[4] = nim->dt ;
4683 nhdr.pixdim[5] = nim->du ; nhdr.pixdim[6] = nim->dv ;
4684 nhdr.pixdim[7] = nim->dw ;
4685
4686 nhdr.datatype = nim->datatype ;
4687 nhdr.bitpix = 8 * nim->nbyper ;
4688
4689 if( nim->cal_max > nim->cal_min ){
4690 nhdr.cal_max = nim->cal_max ;
4691 nhdr.cal_min = nim->cal_min ;
4692 }
4693
4694 if( nim->scl_slope != 0.0 ){
4695 nhdr.scl_slope = nim->scl_slope ;
4696 nhdr.scl_inter = nim->scl_inter ;
4697 }
4698
4699 if( nim->descrip[0] != '\0' ){
4700 memcpy(nhdr.descrip ,nim->descrip ,79) ; nhdr.descrip[79] = '\0' ;
4701 }
4702 if( nim->aux_file[0] != '\0' ){
4703 memcpy(nhdr.aux_file ,nim->aux_file ,23) ; nhdr.aux_file[23] = '\0' ;
4704 }
4705
4706 /**- Load NIFTI specific stuff into the header */
4707
4708 if( nim->nifti_type > NIFTI_FTYPE_ANALYZE ){ /* then not ANALYZE */
4709
4710 if( nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ) strcpy(nhdr.magic,"n+1") ;
4711 else strcpy(nhdr.magic,"ni1") ;
4712
4713 nhdr.pixdim[1] = fabs(nhdr.pixdim[1]) ; nhdr.pixdim[2] = fabs(nhdr.pixdim[2]) ;
4714 nhdr.pixdim[3] = fabs(nhdr.pixdim[3]) ; nhdr.pixdim[4] = fabs(nhdr.pixdim[4]) ;
4715 nhdr.pixdim[5] = fabs(nhdr.pixdim[5]) ; nhdr.pixdim[6] = fabs(nhdr.pixdim[6]) ;
4716 nhdr.pixdim[7] = fabs(nhdr.pixdim[7]) ;
4717
4718 nhdr.intent_code = nim->intent_code ;
4719 nhdr.intent_p1 = nim->intent_p1 ;
4720 nhdr.intent_p2 = nim->intent_p2 ;
4721 nhdr.intent_p3 = nim->intent_p3 ;
4722 if( nim->intent_name[0] != '\0' ){
4723 memcpy(nhdr.intent_name,nim->intent_name,15) ;
4724 nhdr.intent_name[15] = '\0' ;
4725 }
4726
4727 nhdr.vox_offset = (float) nim->iname_offset ;
4728 nhdr.xyzt_units = SPACE_TIME_TO_XYZT( nim->xyz_units, nim->time_units ) ;
4729 nhdr.toffset = nim->toffset ;
4730
4731 if( nim->qform_code > 0 ){
4732 nhdr.qform_code = nim->qform_code ;
4733 nhdr.quatern_b = nim->quatern_b ;
4734 nhdr.quatern_c = nim->quatern_c ;
4735 nhdr.quatern_d = nim->quatern_d ;
4736 nhdr.qoffset_x = nim->qoffset_x ;
4737 nhdr.qoffset_y = nim->qoffset_y ;
4738 nhdr.qoffset_z = nim->qoffset_z ;
4739 nhdr.pixdim[0] = (nim->qfac >= 0.0) ? 1.0 : -1.0 ;
4740 }
4741
4742 if( nim->sform_code > 0 ){
4743 nhdr.sform_code = nim->sform_code ;
4744 nhdr.srow_x[0] = nim->sto_xyz.m[0][0] ;
4745 nhdr.srow_x[1] = nim->sto_xyz.m[0][1] ;
4746 nhdr.srow_x[2] = nim->sto_xyz.m[0][2] ;
4747 nhdr.srow_x[3] = nim->sto_xyz.m[0][3] ;
4748 nhdr.srow_y[0] = nim->sto_xyz.m[1][0] ;
4749 nhdr.srow_y[1] = nim->sto_xyz.m[1][1] ;
4750 nhdr.srow_y[2] = nim->sto_xyz.m[1][2] ;
4751 nhdr.srow_y[3] = nim->sto_xyz.m[1][3] ;
4752 nhdr.srow_z[0] = nim->sto_xyz.m[2][0] ;
4753 nhdr.srow_z[1] = nim->sto_xyz.m[2][1] ;
4754 nhdr.srow_z[2] = nim->sto_xyz.m[2][2] ;
4755 nhdr.srow_z[3] = nim->sto_xyz.m[2][3] ;
4756 }
4757
4758 nhdr.dim_info = FPS_INTO_DIM_INFO( nim->freq_dim ,
4759 nim->phase_dim , nim->slice_dim ) ;
4760 nhdr.slice_code = nim->slice_code ;
4761 nhdr.slice_start = nim->slice_start ;
4762 nhdr.slice_end = nim->slice_end ;
4763 nhdr.slice_duration = nim->slice_duration ;
4764 }
4765
4766 return nhdr;
4767 }
4768
4769
4770 /*----------------------------------------------------------------------*/
4771 /*! \fn int nifti_copy_extensions(nifti_image * nim_dest, nifti_image * nim_src)
4772 \brief copy the nifti1_extension list from src to dest
4773
4774 Duplicate the list of nifti1_extensions. The dest structure must
4775 be clear of extensions.
4776 \return 0 on success, -1 on failure
4777
4778 \sa nifti_add_extension, nifti_free_extensions
4779 */
nifti_copy_extensions(nifti_image * nim_dest,const nifti_image * nim_src)4780 int nifti_copy_extensions(nifti_image * nim_dest, const nifti_image * nim_src)
4781 {
4782 char * data;
4783 size_t bytes;
4784 int c, size, old_size;
4785
4786 if( nim_dest->num_ext > 0 || nim_dest->ext_list != NULL ){
4787 fprintf(stderr,"** will not copy extensions over existing ones\n");
4788 return -1;
4789 }
4790
4791 if( g_opts.debug > 1 )
4792 fprintf(stderr,"+d duplicating %d extension(s)\n", nim_src->num_ext);
4793
4794 if( nim_src->num_ext <= 0 ) return 0;
4795
4796 bytes = nim_src->num_ext * sizeof(nifti1_extension); /* I'm lazy */
4797 nim_dest->ext_list = (nifti1_extension *)malloc(bytes);
4798 if( !nim_dest->ext_list ){
4799 fprintf(stderr,"** failed to allocate %d nifti1_extension structs\n",
4800 nim_src->num_ext);
4801 return -1;
4802 }
4803
4804 /* copy the extension data */
4805 nim_dest->num_ext = 0;
4806 for( c = 0; c < nim_src->num_ext; c++ ){
4807 size = old_size = nim_src->ext_list[c].esize;
4808 if( size & 0xf ) size = (size + 0xf) & ~0xf; /* make multiple of 16 */
4809 if( g_opts.debug > 2 )
4810 fprintf(stderr,"+d dup'ing ext #%d of size %d (from size %d)\n",
4811 c, size, old_size);
4812 data = (char *)calloc(size,sizeof(char)); /* calloc, maybe size > old */
4813 if( !data ){
4814 fprintf(stderr,"** failed to alloc %d bytes for extention\n", size);
4815 if( c == 0 ) { free(nim_dest->ext_list); nim_dest->ext_list = NULL; }
4816 /* otherwise, keep what we have (a.o.t. deleting them all) */
4817 return -1;
4818 }
4819 /* finally, fill the new structure */
4820 nim_dest->ext_list[c].esize = size;
4821 nim_dest->ext_list[c].ecode = nim_src->ext_list[c].ecode;
4822 nim_dest->ext_list[c].edata = data;
4823 memcpy(data, nim_src->ext_list[c].edata, old_size);
4824
4825 nim_dest->num_ext++;
4826 }
4827
4828 return 0;
4829 }
4830
4831
4832 /*----------------------------------------------------------------------*/
4833 /*! compute the total size of all extensions
4834
4835 \return the total of all esize fields
4836
4837 Note that each esize includes 4 bytes for ecode, 4 bytes for esize,
4838 and the bytes used for the data. Each esize also needs to be a
4839 multiple of 16, so it may be greater than the sum of its 3 parts.
4840 *//*--------------------------------------------------------------------*/
nifti_extension_size(nifti_image * nim)4841 int nifti_extension_size(nifti_image *nim)
4842 {
4843 int c, size = 0;
4844
4845 if( !nim || nim->num_ext <= 0 ) return 0;
4846
4847 if( g_opts.debug > 2 ) fprintf(stderr,"-d ext sizes:");
4848
4849 for ( c = 0; c < nim->num_ext; c++ ){
4850 size += nim->ext_list[c].esize;
4851 if( g_opts.debug > 2 ) fprintf(stderr," %d",nim->ext_list[c].esize);
4852 }
4853
4854 if( g_opts.debug > 2 ) fprintf(stderr," (total = %d)\n",size);
4855
4856 return size;
4857 }
4858
4859
4860 /*----------------------------------------------------------------------*/
4861 /*! set the nifti_image iname_offset field, based on nifti_type
4862
4863 - if writing to 2 files, set offset to 0
4864 - if writing to a single NIFTI-1 file, set the offset to
4865 352 + total extension size, then align to 16-byte boundary
4866 - if writing an ASCII header, set offset to -1
4867 *//*--------------------------------------------------------------------*/
nifti_set_iname_offset(nifti_image * nim)4868 void nifti_set_iname_offset(nifti_image *nim)
4869 {
4870 int offset;
4871
4872 switch( nim->nifti_type ){
4873
4874 default: /* writing into 2 files */
4875 /* we only write files with 0 offset in the 2 file format */
4876 nim->iname_offset = 0 ;
4877 break ;
4878
4879 /* NIFTI-1 single binary file - always update */
4880 case NIFTI_FTYPE_NIFTI1_1:
4881 offset = nifti_extension_size(nim)+sizeof(struct nifti_1_header)+4;
4882 /* be sure offset is aligned to a 16 byte boundary */
4883 if ( ( offset % 16 ) != 0 ) offset = ((offset + 0xf) & ~0xf);
4884 if( nim->iname_offset != offset ){
4885 if( g_opts.debug > 1 )
4886 fprintf(stderr,"+d changing offset from %d to %d\n",
4887 nim->iname_offset, offset);
4888 nim->iname_offset = offset;
4889 }
4890 break ;
4891
4892 /* non-standard case: NIFTI-1 ASCII header + binary data (single file) */
4893 case NIFTI_FTYPE_ASCII:
4894 nim->iname_offset = -1 ; /* compute offset from filesize */
4895 break ;
4896 }
4897 }
4898
4899
4900 /*----------------------------------------------------------------------*/
4901 /*! write the nifti_image dataset to disk, optionally including data
4902
4903 This is just a front-end for nifti_image_write_hdr_img2.
4904
4905 \param nim nifti_image to write to disk
4906 \param write_data write options (see nifti_image_write_hdr_img2)
4907 \param opts file open options ("wb" from nifti_image_write)
4908
4909 \sa nifti_image_write, nifti_image_write_hdr_img2, nifti_image_free,
4910 nifti_set_filenames
4911 *//*--------------------------------------------------------------------*/
nifti_image_write_hdr_img(nifti_image * nim,int write_data,const char * opts)4912 znzFile nifti_image_write_hdr_img( nifti_image *nim , int write_data ,
4913 const char* opts )
4914 {
4915 return nifti_image_write_hdr_img2(nim,write_data,opts,NULL,NULL);
4916 }
4917
4918
4919 #undef ERREX
4920 #define ERREX(msg) \
4921 do{ fprintf(stderr,"** ERROR: nifti_image_write_hdr_img: %s\n",(msg)) ; \
4922 return fp ; } while(0)
4923
4924
4925 /* ----------------------------------------------------------------------*/
4926 /*! This writes the header (and optionally the image data) to file
4927 *
4928 * If the image data file is left open it returns a valid znzFile handle.
4929 * It also uses imgfile as the open image file is not null, and modifies
4930 * it inside.
4931 *
4932 * \param nim nifti_image to write to disk
4933 * \param write_opts flags whether to write data and/or close file (see below)
4934 * \param opts file-open options, probably "wb" from nifti_image_write()
4935 * \param imgfile optional open znzFile struct, for writing image data
4936 (may be NULL)
4937 * \param NBL optional nifti_brick_list, containing the image data
4938 (may be NULL)
4939 *
4940 * Values for write_opts mode are based on two binary flags
4941 * ( 0/1 for no-write/write data, and 0/2 for close/leave-open files ) :
4942 * - 0 = do not write data and close (do not open data file)
4943 * - 1 = write data and close
4944 * - 2 = do not write data and leave data file open
4945 * - 3 = write data and leave data file open
4946 *
4947 * \sa nifti_image_write, nifti_image_write_hdr_img, nifti_image_free,
4948 * nifti_set_filenames
4949 *//*---------------------------------------------------------------------*/
nifti_image_write_hdr_img2(nifti_image * nim,int write_opts,const char * opts,znzFile imgfile,const nifti_brick_list * NBL)4950 znzFile nifti_image_write_hdr_img2(nifti_image *nim, int write_opts,
4951 const char * opts, znzFile imgfile, const nifti_brick_list * NBL)
4952 {
4953 struct nifti_1_header nhdr ;
4954 znzFile fp=NULL;
4955 size_t ss ;
4956 int write_data, leave_open;
4957 char func[] = { "nifti_image_write_hdr_img2" };
4958
4959 write_data = write_opts & 1; /* just separate the bits now */
4960 leave_open = write_opts & 2;
4961
4962 if( ! nim ) ERREX("NULL input") ;
4963 if( ! nifti_validfilename(nim->fname) ) ERREX("bad fname input") ;
4964 if( write_data && ! nim->data && ! NBL ) ERREX("no image data") ;
4965
4966 nifti_set_iname_offset(nim);
4967
4968 if( g_opts.debug > 1 ){
4969 fprintf(stderr,"-d writing nifti file '%s'...\n", nim->fname);
4970 if( g_opts.debug > 2 )
4971 fprintf(stderr,"-d nifti type %d, offset %d\n",
4972 nim->nifti_type, nim->iname_offset);
4973 }
4974
4975 if( nim->nifti_type == NIFTI_FTYPE_ASCII ) /* non-standard case */
4976 return nifti_write_ascii_image(nim,NBL,opts,write_data,leave_open);
4977
4978 nhdr = nifti_convert_nim2nhdr(nim); /* create the nifti1_header struct */
4979
4980 /* if writing to 2 files, make sure iname is set and different from fname */
4981 if( nim->nifti_type != NIFTI_FTYPE_NIFTI1_1 ){
4982 if( nim->iname && strcmp(nim->iname,nim->fname) == 0 ){
4983 free(nim->iname) ; nim->iname = NULL ;
4984 }
4985 if( nim->iname == NULL ){ /* then make a new one */
4986 nim->iname = nifti_makeimgname(nim->fname,nim->nifti_type,0,0);
4987 if( nim->iname == NULL ) return NULL;
4988 }
4989 }
4990
4991 /* if we have an imgfile and will write the header there, use it */
4992 if( ! znz_isnull(imgfile) && nim->nifti_type == NIFTI_FTYPE_NIFTI1_1 ){
4993 if( g_opts.debug > 2 ) fprintf(stderr,"+d using passed file for hdr\n");
4994 fp = imgfile;
4995 }
4996 else {
4997 if( g_opts.debug > 2 )
4998 fprintf(stderr,"+d opening output file '%s'\n",nim->fname);
4999 fp = znzopen( nim->fname , opts , nifti_is_gzfile(nim->fname) ) ;
5000 if( znz_isnull(fp) ){
5001 LNI_FERR(func,"cannot open output file",nim->fname);
5002 return fp;
5003 }
5004 }
5005
5006 /* write the header and extensions */
5007
5008 ss = znzwrite(&nhdr , 1 , sizeof(nhdr) , fp); /* write header */
5009 if( ss < sizeof(nhdr) ){
5010 LNI_FERR(func,"bad header write to output file",nim->fname);
5011 znzclose(fp); return fp;
5012 }
5013
5014 /* partial file exists, and errors have been printed, so ignore return */
5015 if( nim->nifti_type != NIFTI_FTYPE_ANALYZE )
5016 (void)nifti_write_extensions(fp,nim);
5017
5018 /* if the header is all we want, we are done */
5019 if( ! write_data && ! leave_open ){
5020 if( g_opts.debug > 2 ) fprintf(stderr,"-d header is all we want: done\n");
5021 znzclose(fp); return(fp);
5022 }
5023
5024 if( nim->nifti_type != NIFTI_FTYPE_NIFTI1_1 ){ /* get a new file pointer */
5025 znzclose(fp); /* first, close header file */
5026 if( ! znz_isnull(imgfile) ){
5027 if(g_opts.debug > 2) fprintf(stderr,"+d using passed file for img\n");
5028 fp = imgfile;
5029 }
5030 else {
5031 if( g_opts.debug > 2 )
5032 fprintf(stderr,"+d opening img file '%s'\n", nim->iname);
5033 fp = znzopen( nim->iname , opts , nifti_is_gzfile(nim->iname) ) ;
5034 if( znz_isnull(fp) ) ERREX("cannot open image file") ;
5035 }
5036 }
5037
5038 znzseek(fp, nim->iname_offset, SEEK_SET); /* in any case, seek to offset */
5039
5040 if( write_data ) nifti_write_all_data(fp,nim,NBL);
5041 if( ! leave_open ) znzclose(fp);
5042
5043 return fp;
5044 }
5045
5046
5047 /*----------------------------------------------------------------------*/
5048 /*! write a nifti_image to disk in ASCII format
5049 *//*--------------------------------------------------------------------*/
nifti_write_ascii_image(nifti_image * nim,const nifti_brick_list * NBL,const char * opts,int write_data,int leave_open)5050 znzFile nifti_write_ascii_image(nifti_image *nim, const nifti_brick_list * NBL,
5051 const char *opts, int write_data, int leave_open)
5052 {
5053 znzFile fp;
5054 char * hstr;
5055
5056 hstr = nifti_image_to_ascii( nim ) ; /* get header in ASCII form */
5057 if( ! hstr ){ fprintf(stderr,"** failed image_to_ascii()\n"); return NULL; }
5058
5059 fp = znzopen( nim->fname , opts , nifti_is_gzfile(nim->fname) ) ;
5060 if( znz_isnull(fp) ){
5061 free(hstr);
5062 fprintf(stderr,"** failed to open '%s' for ascii write\n",nim->fname);
5063 return fp;
5064 }
5065
5066 znzputs(hstr,fp); /* header */
5067 nifti_write_extensions(fp,nim); /* extensions */
5068
5069 if ( write_data ) { nifti_write_all_data(fp,nim,NBL); } /* data */
5070 if ( ! leave_open ) { znzclose(fp); }
5071 free(hstr);
5072 return fp; /* returned but may be closed */
5073 }
5074
5075
5076 /*--------------------------------------------------------------------------*/
5077 /*! Write a nifti_image to disk.
5078
5079 Since data is properly byte-swapped upon reading, it is assumed
5080 to be in the byte-order of the current CPU at write time. Thus,
5081 nim->byte_order should match that of the current CPU. Note that
5082 the nifti_set_filenames() function takes the flag, set_byte_order.
5083
5084 The following fields of nim affect how the output appears:
5085 - nifti_type = 0 ==> ANALYZE-7.5 format file pair will be written
5086 - nifti_type = 1 ==> NIFTI-1 format single file will be written
5087 (data offset will be 352+extensions)
5088 - nifti_type = 2 ==> NIFTI_1 format file pair will be written
5089 - nifti_type = 3 ==> NIFTI_1 ASCII single file will be written
5090 - fname is the name of the output file (header or header+data)
5091 - if a file pair is being written, iname is the name of the data file
5092 - existing files WILL be overwritten with extreme prejudice
5093 - if qform_code > 0, the quatern_*, qoffset_*, and qfac fields determine
5094 the qform output, NOT the qto_xyz matrix; if you want to compute these
5095 fields from the qto_xyz matrix, you can use the utility function
5096 nifti_mat44_to_quatern()
5097
5098 \sa nifti_image_write_bricks, nifti_image_free, nifti_set_filenames,
5099 nifti_image_write_hdr_img
5100 *//*------------------------------------------------------------------------*/
nifti_image_write(nifti_image * nim)5101 void nifti_image_write( nifti_image *nim )
5102 {
5103 znzFile fp = nifti_image_write_hdr_img(nim,1,"wb");
5104 if( fp ){
5105 if( g_opts.debug > 2 ) fprintf(stderr,"-d niw: done with znzFile\n");
5106 free(fp);
5107 }
5108 if( g_opts.debug > 1 ) fprintf(stderr,"-d nifti_image_write: done\n");
5109 }
5110
5111
5112 /*----------------------------------------------------------------------*/
5113 /*! similar to nifti_image_write, but data is in NBL struct, not nim->data
5114
5115 \sa nifti_image_write, nifti_image_free, nifti_set_filenames, nifti_free_NBL
5116 *//*--------------------------------------------------------------------*/
nifti_image_write_bricks(nifti_image * nim,const nifti_brick_list * NBL)5117 void nifti_image_write_bricks( nifti_image *nim, const nifti_brick_list * NBL )
5118 {
5119 znzFile fp = nifti_image_write_hdr_img2(nim,1,"wb",NULL,NBL);
5120 if( fp ){
5121 if( g_opts.debug > 2 ) fprintf(stderr,"-d niwb: done with znzFile\n");
5122 free(fp);
5123 }
5124 if( g_opts.debug > 1 ) fprintf(stderr,"-d niwb: done writing bricks\n");
5125 }
5126
5127
5128 /*----------------------------------------------------------------------*/
5129 /*! copy the nifti_image structure, without data
5130
5131 Duplicate the structure, including fname, iname and extensions.
5132 Leave the data pointer as NULL.
5133 *//*--------------------------------------------------------------------*/
nifti_copy_nim_info(const nifti_image * src)5134 nifti_image * nifti_copy_nim_info(const nifti_image * src)
5135 {
5136 nifti_image *dest;
5137 dest = (nifti_image *)calloc(1,sizeof(nifti_image));
5138 if( !dest ){
5139 fprintf(stderr,"** NCNI: failed to alloc nifti_image\n");
5140 return NULL;
5141 }
5142 memcpy(dest, src, sizeof(nifti_image));
5143 if( src->fname ) dest->fname = nifti_strdup(src->fname);
5144 if( src->iname ) dest->iname = nifti_strdup(src->iname);
5145 dest->num_ext = 0;
5146 dest->ext_list = NULL;
5147 /* errors will be printed in NCE(), continue in either case */
5148 (void)nifti_copy_extensions(dest, src);
5149
5150 dest->data = NULL;
5151
5152 return dest;
5153 }
5154
5155
5156 /*------------------------------------------------------------------------*/
5157 /* Un-escape a C string in place -- that is, convert XML escape sequences
5158 back into their characters. (This can be done in place since the
5159 replacement is always smaller than the input.) Escapes recognized are:
5160 - < -> <
5161 - > -> >
5162 - " -> "
5163 - ' -> '
5164 - & -> &
5165 Also replace CR LF pair (Microsoft), or CR alone (Macintosh) with
5166 LF (Unix), per the XML standard.
5167 Return value is number of replacements made (if you care).
5168 --------------------------------------------------------------------------*/
5169
5170 #undef CR
5171 #undef LF
5172 #define CR 0x0D
5173 #define LF 0x0A
5174
unescape_string(char * str)5175 static int unescape_string( char *str )
5176 {
5177 int ii,jj , nn,ll ;
5178
5179 if( str == NULL ) return 0 ; /* no string? */
5180 ll = strlen(str) ; if( ll == 0 ) return 0 ;
5181
5182 /* scan for escapes: &something; */
5183
5184 for( ii=jj=nn=0 ; ii<ll ; ii++,jj++ ){ /* scan at ii; results go in at jj */
5185
5186 if( str[ii] == '&' ){ /* start of escape? */
5187
5188 if( ii+3 < ll && /* < */
5189 str[ii+1] == 'l' &&
5190 str[ii+2] == 't' &&
5191 str[ii+3] == ';' ){ str[jj] = '<' ; ii += 3 ; nn++ ; }
5192
5193 else if( ii+3 < ll && /* > */
5194 str[ii+1] == 'g' &&
5195 str[ii+2] == 't' &&
5196 str[ii+3] == ';' ){ str[jj] = '>' ; ii += 3 ; nn++ ; }
5197
5198 else if( ii+5 < ll && /* " */
5199 str[ii+1] == 'q' &&
5200 str[ii+2] == 'u' &&
5201 str[ii+3] == 'o' &&
5202 str[ii+4] == 't' &&
5203 str[ii+5] == ';' ){ str[jj] = '"' ; ii += 5 ; nn++ ; }
5204
5205 else if( ii+5 < ll && /* ' */
5206 str[ii+1] == 'a' &&
5207 str[ii+2] == 'p' &&
5208 str[ii+3] == 'o' &&
5209 str[ii+4] == 's' &&
5210 str[ii+5] == ';' ){ str[jj] = '\'' ; ii += 5 ; nn++ ; }
5211
5212 else if( ii+4 < ll && /* & */
5213 str[ii+1] == 'a' &&
5214 str[ii+2] == 'm' &&
5215 str[ii+3] == 'p' &&
5216 str[ii+4] == ';' ){ str[jj] = '&' ; ii += 4 ; nn++ ; }
5217
5218 /* although the comments above don't mention it,
5219 we also look for XML style numeric escapes
5220 of the forms   (decimal) and ý (hex) */
5221
5222 else if( ii+3 < ll &&
5223 str[ii+1] == '#' &&
5224 isdigit(str[ii+2]) ){ /* &#dec; */
5225
5226 unsigned int val='?' ; int kk=ii+3 ;
5227 while( kk < ll && kk != ';' ) kk++ ;
5228 sscanf( str+ii+2 , "%u" , &val ) ;
5229 str[jj] = (char) val ; ii = kk ; nn++ ;
5230 }
5231
5232 else if( ii+4 < ll &&
5233 str[ii+1] == '#' &&
5234 str[ii+2] == 'x' &&
5235 isxdigit(str[ii+3]) ){ /* &#hex; */
5236
5237 unsigned int val='?' ; int kk=ii+4 ;
5238 while( kk < ll && kk != ';' ) kk++ ;
5239 sscanf( str+ii+3 , "%x" , &val ) ;
5240 str[jj] = (char) val ; ii = kk ; nn++ ;
5241 }
5242
5243 /* didn't start a recognized escape, so just copy as normal */
5244
5245 else if( jj < ii ){ str[jj] = str[ii] ; }
5246
5247 } else if( str[ii] == CR ) { /* is a carriage return */
5248
5249 if( str[ii+1] == LF ){ str[jj] = LF ; ii++ ; nn++ ; } /* CR LF */
5250 else { str[jj] = LF ; ; nn++ ; } /* CR only */
5251
5252 } else { /* is a normal character, just copy to output */
5253
5254 if( jj < ii ){ str[jj] = str[ii] ; }
5255 }
5256
5257 /* at this point, ii=index of last character used up in scan
5258 jj=index of last character written to (jj <= ii) */
5259 }
5260
5261 if( jj < ll ) str[jj] = '\0' ; /* end string properly */
5262
5263 return nn ;
5264 }
5265
5266 /*------------------------------------------------------------------------*/
5267 /* Quotize (and escapize) one string, returning a new string.
5268 Approximately speaking, this is the inverse of unescape_string().
5269 The result should be free()-ed when you are done with it.
5270 --------------------------------------------------------------------------*/
5271
escapize_string(const char * str)5272 static char *escapize_string( const char * str )
5273 {
5274 int ii,jj , lstr,lout ;
5275 char *out ;
5276
5277 if( str == NULL || (lstr=strlen(str)) == 0 ){ /* 0 length */
5278 out = nifti_strdup("''") ; return out ; /* string?? */
5279 }
5280
5281 lout = 4 ; /* initialize length of output */
5282 for( ii=0 ; ii < lstr ; ii++ ){ /* count characters for output */
5283 switch( str[ii] ){
5284 case '&': lout += 5 ; break ; /* replace '&' with "&" */
5285
5286 case '<':
5287 case '>': lout += 4 ; break ; /* replace '<' with "<" */
5288
5289 case '"' :
5290 case '\'': lout += 6 ; break ; /* replace '"' with """ */
5291
5292 case CR:
5293 case LF: lout += 6 ; break ; /* replace CR with "
"
5294 LF with "
" */
5295
5296 default: lout++ ; break ; /* copy all other chars */
5297 }
5298 }
5299 out = (char *)calloc(1,lout) ; /* allocate output string */
5300 if( !out ){
5301 fprintf(stderr,"** escapize_string: failed to alloc %d bytes\n",lout);
5302 return NULL;
5303 }
5304 out[0] = '\'' ; /* opening quote mark */
5305 for( ii=0,jj=1 ; ii < lstr ; ii++ ){
5306 switch( str[ii] ){
5307 default: out[jj++] = str[ii] ; break ; /* normal characters */
5308
5309 case '&': memcpy(out+jj,"&",5) ; jj+=5 ; break ;
5310
5311 case '<': memcpy(out+jj,"<",4) ; jj+=4 ; break ;
5312 case '>': memcpy(out+jj,">",4) ; jj+=4 ; break ;
5313
5314 case '"' : memcpy(out+jj,""",6) ; jj+=6 ; break ;
5315
5316 case '\'': memcpy(out+jj,"'",6) ; jj+=6 ; break ;
5317
5318 case CR: memcpy(out+jj,"
",6) ; jj+=6 ; break ;
5319 case LF: memcpy(out+jj,"
",6) ; jj+=6 ; break ;
5320 }
5321 }
5322 out[jj++] = '\'' ; /* closing quote mark */
5323 out[jj] = '\0' ; /* terminate the string */
5324 return out ;
5325 }
5326
5327 /*---------------------------------------------------------------------------*/
5328 /*! Dump the information in a NIFTI image header to an XML-ish ASCII string
5329 that can later be converted back into a NIFTI header in
5330 nifti_image_from_ascii().
5331
5332 The resulting string can be free()-ed when you are done with it.
5333 *//*-------------------------------------------------------------------------*/
nifti_image_to_ascii(const nifti_image * nim)5334 char *nifti_image_to_ascii( const nifti_image *nim )
5335 {
5336 char *buf , *ebuf ; int nbuf ;
5337
5338 if( nim == NULL ) return NULL ; /* stupid caller */
5339
5340 buf = (char *)calloc(1,65534); nbuf = 0; /* longer than needed, to be safe */
5341 if( !buf ){
5342 fprintf(stderr,"** NITA: failed to alloc %d bytes\n",65534);
5343 return NULL;
5344 }
5345
5346 sprintf( buf , "<nifti_image\n" ) ; /* XML-ish opener */
5347
5348 sprintf( buf+strlen(buf) , " nifti_type = '%s'\n" ,
5349 (nim->nifti_type == NIFTI_FTYPE_NIFTI1_1) ? "NIFTI-1+"
5350 :(nim->nifti_type == NIFTI_FTYPE_NIFTI1_2) ? "NIFTI-1"
5351 :(nim->nifti_type == NIFTI_FTYPE_ASCII ) ? "NIFTI-1A"
5352 : "ANALYZE-7.5" ) ;
5353
5354 /** Strings that we don't control (filenames, etc.) that might
5355 contain "weird" characters (like quotes) are "escaped":
5356 - A few special characters are replaced by XML-style escapes, using
5357 the function escapize_string().
5358 - On input, function unescape_string() reverses this process.
5359 - The result is that the NIFTI ASCII-format header is XML-compliant. */
5360
5361 ebuf = escapize_string(nim->fname) ;
5362 sprintf( buf+strlen(buf) , " header_filename = %s\n",ebuf); free(ebuf);
5363
5364 ebuf = escapize_string(nim->iname) ;
5365 sprintf( buf+strlen(buf) , " image_filename = %s\n", ebuf); free(ebuf);
5366
5367 sprintf( buf+strlen(buf) , " image_offset = '%d'\n" , nim->iname_offset );
5368
5369 sprintf( buf+strlen(buf), " ndim = '%d'\n", nim->ndim);
5370 sprintf( buf+strlen(buf), " nx = '%d'\n", nim->nx );
5371 if( nim->ndim > 1 ) sprintf( buf+strlen(buf), " ny = '%d'\n", nim->ny );
5372 if( nim->ndim > 2 ) sprintf( buf+strlen(buf), " nz = '%d'\n", nim->nz );
5373 if( nim->ndim > 3 ) sprintf( buf+strlen(buf), " nt = '%d'\n", nim->nt );
5374 if( nim->ndim > 4 ) sprintf( buf+strlen(buf), " nu = '%d'\n", nim->nu );
5375 if( nim->ndim > 5 ) sprintf( buf+strlen(buf), " nv = '%d'\n", nim->nv );
5376 if( nim->ndim > 6 ) sprintf( buf+strlen(buf), " nw = '%d'\n", nim->nw );
5377 sprintf( buf+strlen(buf), " dx = '%g'\n", nim->dx );
5378 if( nim->ndim > 1 ) sprintf( buf+strlen(buf), " dy = '%g'\n", nim->dy );
5379 if( nim->ndim > 2 ) sprintf( buf+strlen(buf), " dz = '%g'\n", nim->dz );
5380 if( nim->ndim > 3 ) sprintf( buf+strlen(buf), " dt = '%g'\n", nim->dt );
5381 if( nim->ndim > 4 ) sprintf( buf+strlen(buf), " du = '%g'\n", nim->du );
5382 if( nim->ndim > 5 ) sprintf( buf+strlen(buf), " dv = '%g'\n", nim->dv );
5383 if( nim->ndim > 6 ) sprintf( buf+strlen(buf), " dw = '%g'\n", nim->dw );
5384
5385 sprintf( buf+strlen(buf) , " datatype = '%d'\n" , nim->datatype ) ;
5386 sprintf( buf+strlen(buf) , " datatype_name = '%s'\n" ,
5387 nifti_datatype_string(nim->datatype) ) ;
5388
5389 sprintf( buf+strlen(buf) , " nvox = '%d'\n" , nim->nvox ) ;
5390 sprintf( buf+strlen(buf) , " nbyper = '%d'\n" , nim->nbyper ) ;
5391
5392 sprintf( buf+strlen(buf) , " byteorder = '%s'\n" ,
5393 (nim->byteorder==MSB_FIRST) ? "MSB_FIRST" : "LSB_FIRST" ) ;
5394
5395 if( nim->cal_min < nim->cal_max ){
5396 sprintf( buf+strlen(buf) , " cal_min = '%g'\n", nim->cal_min ) ;
5397 sprintf( buf+strlen(buf) , " cal_max = '%g'\n", nim->cal_max ) ;
5398 }
5399
5400 if( nim->scl_slope != 0.0 ){
5401 sprintf( buf+strlen(buf) , " scl_slope = '%g'\n" , nim->scl_slope ) ;
5402 sprintf( buf+strlen(buf) , " scl_inter = '%g'\n" , nim->scl_inter ) ;
5403 }
5404
5405 if( nim->intent_code > 0 ){
5406 sprintf( buf+strlen(buf) , " intent_code = '%d'\n", nim->intent_code ) ;
5407 sprintf( buf+strlen(buf) , " intent_code_name = '%s'\n" ,
5408 nifti_intent_string(nim->intent_code) ) ;
5409 sprintf( buf+strlen(buf) , " intent_p1 = '%g'\n" , nim->intent_p1 ) ;
5410 sprintf( buf+strlen(buf) , " intent_p2 = '%g'\n" , nim->intent_p2 ) ;
5411 sprintf( buf+strlen(buf) , " intent_p3 = '%g'\n" , nim->intent_p3 ) ;
5412
5413 if( nim->intent_name[0] != '\0' ){
5414 ebuf = escapize_string(nim->intent_name) ;
5415 sprintf( buf+strlen(buf) , " intent_name = %s\n",ebuf) ;
5416 free(ebuf) ;
5417 }
5418 }
5419
5420 if( nim->toffset != 0.0 )
5421 sprintf( buf+strlen(buf) , " toffset = '%g'\n",nim->toffset ) ;
5422
5423 if( nim->xyz_units > 0 )
5424 sprintf( buf+strlen(buf) ,
5425 " xyz_units = '%d'\n"
5426 " xyz_units_name = '%s'\n" ,
5427 nim->xyz_units , nifti_units_string(nim->xyz_units) ) ;
5428
5429 if( nim->time_units > 0 )
5430 sprintf( buf+strlen(buf) ,
5431 " time_units = '%d'\n"
5432 " time_units_name = '%s'\n" ,
5433 nim->time_units , nifti_units_string(nim->time_units) ) ;
5434
5435 if( nim->freq_dim > 0 )
5436 sprintf( buf+strlen(buf) , " freq_dim = '%d'\n",nim->freq_dim ) ;
5437 if( nim->phase_dim > 0 )
5438 sprintf( buf+strlen(buf) , " phase_dim = '%d'\n",nim->phase_dim ) ;
5439 if( nim->slice_dim > 0 )
5440 sprintf( buf+strlen(buf) , " slice_dim = '%d'\n",nim->slice_dim ) ;
5441 if( nim->slice_code > 0 )
5442 sprintf( buf+strlen(buf) ,
5443 " slice_code = '%d'\n"
5444 " slice_code_name = '%s'\n" ,
5445 nim->slice_code , nifti_slice_string(nim->slice_code) ) ;
5446 if( nim->slice_start >= 0 && nim->slice_end > nim->slice_start )
5447 sprintf( buf+strlen(buf) ,
5448 " slice_start = '%d'\n"
5449 " slice_end = '%d'\n" , nim->slice_start , nim->slice_end ) ;
5450 if( nim->slice_duration != 0.0 )
5451 sprintf( buf+strlen(buf) , " slice_duration = '%g'\n",
5452 nim->slice_duration ) ;
5453
5454 if( nim->descrip[0] != '\0' ){
5455 ebuf = escapize_string(nim->descrip) ;
5456 sprintf( buf+strlen(buf) , " descrip = %s\n",ebuf) ;
5457 free(ebuf) ;
5458 }
5459
5460 if( nim->aux_file[0] != '\0' ){
5461 ebuf = escapize_string(nim->aux_file) ;
5462 sprintf( buf+strlen(buf) , " aux_file = %s\n",ebuf) ;
5463 free(ebuf) ;
5464 }
5465
5466 if( nim->qform_code > 0 ){
5467 int i,j,k ;
5468
5469 sprintf( buf+strlen(buf) ,
5470 " qform_code = '%d'\n"
5471 " qform_code_name = '%s'\n"
5472 " qto_xyz_matrix = '%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g'\n" ,
5473 nim->qform_code , nifti_xform_string(nim->qform_code) ,
5474 nim->qto_xyz.m[0][0] , nim->qto_xyz.m[0][1] ,
5475 nim->qto_xyz.m[0][2] , nim->qto_xyz.m[0][3] ,
5476 nim->qto_xyz.m[1][0] , nim->qto_xyz.m[1][1] ,
5477 nim->qto_xyz.m[1][2] , nim->qto_xyz.m[1][3] ,
5478 nim->qto_xyz.m[2][0] , nim->qto_xyz.m[2][1] ,
5479 nim->qto_xyz.m[2][2] , nim->qto_xyz.m[2][3] ,
5480 nim->qto_xyz.m[3][0] , nim->qto_xyz.m[3][1] ,
5481 nim->qto_xyz.m[3][2] , nim->qto_xyz.m[3][3] ) ;
5482
5483 sprintf( buf+strlen(buf) ,
5484 " qto_ijk_matrix = '%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g'\n" ,
5485 nim->qto_ijk.m[0][0] , nim->qto_ijk.m[0][1] ,
5486 nim->qto_ijk.m[0][2] , nim->qto_ijk.m[0][3] ,
5487 nim->qto_ijk.m[1][0] , nim->qto_ijk.m[1][1] ,
5488 nim->qto_ijk.m[1][2] , nim->qto_ijk.m[1][3] ,
5489 nim->qto_ijk.m[2][0] , nim->qto_ijk.m[2][1] ,
5490 nim->qto_ijk.m[2][2] , nim->qto_ijk.m[2][3] ,
5491 nim->qto_ijk.m[3][0] , nim->qto_ijk.m[3][1] ,
5492 nim->qto_ijk.m[3][2] , nim->qto_ijk.m[3][3] ) ;
5493
5494 sprintf( buf+strlen(buf) ,
5495 " quatern_b = '%g'\n"
5496 " quatern_c = '%g'\n"
5497 " quatern_d = '%g'\n"
5498 " qoffset_x = '%g'\n"
5499 " qoffset_y = '%g'\n"
5500 " qoffset_z = '%g'\n"
5501 " qfac = '%g'\n" ,
5502 nim->quatern_b , nim->quatern_c , nim->quatern_d ,
5503 nim->qoffset_x , nim->qoffset_y , nim->qoffset_z , nim->qfac ) ;
5504
5505 nifti_mat44_to_orientation( nim->qto_xyz , &i,&j,&k ) ;
5506 if( i > 0 && j > 0 && k > 0 )
5507 sprintf( buf+strlen(buf) ,
5508 " qform_i_orientation = '%s'\n"
5509 " qform_j_orientation = '%s'\n"
5510 " qform_k_orientation = '%s'\n" ,
5511 nifti_orientation_string(i) ,
5512 nifti_orientation_string(j) ,
5513 nifti_orientation_string(k) ) ;
5514 }
5515
5516 if( nim->sform_code > 0 ){
5517 int i,j,k ;
5518
5519 sprintf( buf+strlen(buf) ,
5520 " sform_code = '%d'\n"
5521 " sform_code_name = '%s'\n"
5522 " sto_xyz_matrix = '%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g'\n" ,
5523 nim->sform_code , nifti_xform_string(nim->sform_code) ,
5524 nim->sto_xyz.m[0][0] , nim->sto_xyz.m[0][1] ,
5525 nim->sto_xyz.m[0][2] , nim->sto_xyz.m[0][3] ,
5526 nim->sto_xyz.m[1][0] , nim->sto_xyz.m[1][1] ,
5527 nim->sto_xyz.m[1][2] , nim->sto_xyz.m[1][3] ,
5528 nim->sto_xyz.m[2][0] , nim->sto_xyz.m[2][1] ,
5529 nim->sto_xyz.m[2][2] , nim->sto_xyz.m[2][3] ,
5530 nim->sto_xyz.m[3][0] , nim->sto_xyz.m[3][1] ,
5531 nim->sto_xyz.m[3][2] , nim->sto_xyz.m[3][3] ) ;
5532
5533 sprintf( buf+strlen(buf) ,
5534 " sto_ijk matrix = '%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g'\n" ,
5535 nim->sto_ijk.m[0][0] , nim->sto_ijk.m[0][1] ,
5536 nim->sto_ijk.m[0][2] , nim->sto_ijk.m[0][3] ,
5537 nim->sto_ijk.m[1][0] , nim->sto_ijk.m[1][1] ,
5538 nim->sto_ijk.m[1][2] , nim->sto_ijk.m[1][3] ,
5539 nim->sto_ijk.m[2][0] , nim->sto_ijk.m[2][1] ,
5540 nim->sto_ijk.m[2][2] , nim->sto_ijk.m[2][3] ,
5541 nim->sto_ijk.m[3][0] , nim->sto_ijk.m[3][1] ,
5542 nim->sto_ijk.m[3][2] , nim->sto_ijk.m[3][3] ) ;
5543
5544 nifti_mat44_to_orientation( nim->sto_xyz , &i,&j,&k ) ;
5545 if( i > 0 && j > 0 && k > 0 )
5546 sprintf( buf+strlen(buf) ,
5547 " sform_i_orientation = '%s'\n"
5548 " sform_j_orientation = '%s'\n"
5549 " sform_k_orientation = '%s'\n" ,
5550 nifti_orientation_string(i) ,
5551 nifti_orientation_string(j) ,
5552 nifti_orientation_string(k) ) ;
5553 }
5554
5555 sprintf( buf+strlen(buf) , " num_ext = '%d'\n", nim->num_ext ) ;
5556
5557 sprintf( buf+strlen(buf) , "/>\n" ) ; /* XML-ish closer */
5558
5559 nbuf = strlen(buf) ;
5560 buf = (char *)realloc((void *)buf, nbuf+1); /* cut back to proper length */
5561 if( !buf ) fprintf(stderr,"** NITA: failed to realloc %d bytes\n",nbuf+1);
5562 return buf ;
5563 }
5564
5565 /*---------------------------------------------------------------------------*/
5566
5567 /*----------------------------------------------------------------------*/
5568 /*! get the byte order for this CPU
5569
5570 - LSB_FIRST means least significant byte, first (little endian)
5571 - MSB_FIRST means most significant byte, first (big endian)
5572 *//*--------------------------------------------------------------------*/
nifti_short_order(void)5573 int nifti_short_order(void) /* determine this CPU's byte order */
5574 {
5575 union { unsigned char bb[2] ;
5576 short ss ; } fred ;
5577
5578 fred.bb[0] = 1 ; fred.bb[1] = 0 ;
5579
5580 return (fred.ss == 1) ? LSB_FIRST : MSB_FIRST ;
5581 }
5582
5583 /*---------------------------------------------------------------------------*/
5584
5585 #undef QQNUM
5586 #undef QNUM
5587 #undef QSTR
5588
5589 /* macro to check lhs string against "n1"; if it matches,
5590 interpret rhs string as a number, and put it into nim->"n2" */
5591
5592 #define QQNUM(n1,n2) if( strcmp(lhs,#n1)==0 ) nim->n2=strtod(rhs,NULL)
5593
5594 /* same, but where "n1" == "n2" */
5595
5596 #define QNUM(nam) QQNUM(nam,nam)
5597
5598 /* macro to check lhs string against "nam"; if it matches,
5599 put rhs string into nim->"nam" string, with max length = "ml" */
5600
5601 #define QSTR(nam,ml) if( strcmp(lhs,#nam) == 0 ) \
5602 strncpy(nim->nam,rhs,ml), nim->nam[ml]='\0'
5603
5604 /*---------------------------------------------------------------------------*/
5605 /*! Take an XML-ish ASCII string and create a NIFTI image header to match.
5606
5607 NULL is returned if enough information isn't present in the input string.
5608 - The image data can later be loaded with nifti_image_load().
5609 - The struct returned here can be liberated with nifti_image_free().
5610 - Not a lot of error checking is done here to make sure that the
5611 input values are reasonable!
5612 *//*-------------------------------------------------------------------------*/
nifti_image_from_ascii(const char * str,int * bytes_read)5613 nifti_image *nifti_image_from_ascii( const char *str, int * bytes_read )
5614 {
5615 char lhs[1024] , rhs[1024] ;
5616 int ii , spos, nn , slen ;
5617 nifti_image *nim ; /* will be output */
5618
5619 if( str == NULL || *str == '\0' ) return NULL ; /* bad input!? */
5620
5621 /* scan for opening string */
5622
5623 spos = 0 ; slen = strlen(str) ;
5624 ii = sscanf( str+spos , "%1023s%n" , lhs , &nn ) ; spos += nn ;
5625 if( ii == 0 || strcmp(lhs,"<nifti_image") != 0 ) return NULL ;
5626
5627 /* create empty image struct */
5628
5629 nim = (nifti_image *) calloc( 1 , sizeof(nifti_image) ) ;
5630 if( !nim ){
5631 fprintf(stderr,"** NIFA: failed to alloc nifti_image\n");
5632 return NULL;
5633 }
5634
5635 nim->nx = nim->ny = nim->nz = nim->nt
5636 = nim->nu = nim->nv = nim->nw = 1 ;
5637 nim->dx = nim->dy = nim->dz = nim->dt
5638 = nim->du = nim->dv = nim->dw = nim->qfac = 1.0 ;
5639
5640 nim->byteorder = nifti_short_order() ;
5641
5642 /* starting at str[spos], scan for "equations" of the form
5643 lhs = 'rhs'
5644 and assign rhs values into the struct component named by lhs */
5645
5646 while(1){
5647
5648 while( isspace(str[spos]) ) spos++ ; /* skip whitespace */
5649 if( str[spos] == '\0' ) break ; /* end of string? */
5650
5651 /* get lhs string */
5652
5653 ii = sscanf( str+spos , "%1023s%n" , lhs , &nn ) ; spos += nn ;
5654 if( ii == 0 || strcmp(lhs,"/>") == 0 ) break ; /* end of input? */
5655
5656 /* skip whitespace and the '=' marker */
5657
5658 while( isspace(str[spos]) || str[spos] == '=' ) spos++ ;
5659 if( str[spos] == '\0' ) break ; /* end of string? */
5660
5661 /* if next character is a quote ', copy everything up to next '
5662 otherwise, copy everything up to next nonblank */
5663
5664 if( str[spos] == '\'' ){
5665 ii = spos+1 ;
5666 while( str[ii] != '\0' && str[ii] != '\'' ) ii++ ;
5667 nn = ii-spos-1 ; if( nn > 1023 ) nn = 1023 ;
5668 memcpy(rhs,str+spos+1,nn) ; rhs[nn] = '\0' ;
5669 spos = (str[ii] == '\'') ? ii+1 : ii ;
5670 } else {
5671 ii = sscanf( str+spos , "%1023s%n" , rhs , &nn ) ; spos += nn ;
5672 if( ii == 0 ) break ; /* nothing found? */
5673 }
5674 unescape_string(rhs) ; /* remove any XML escape sequences */
5675
5676 /* Now can do the assignment, based on lhs string.
5677 Start with special cases that don't fit the QNUM/QSTR macros. */
5678
5679 if( strcmp(lhs,"nifti_type") == 0 ){
5680 if( strcmp(rhs,"ANALYZE-7.5") == 0 )
5681 nim->nifti_type = NIFTI_FTYPE_ANALYZE ;
5682 else if( strcmp(rhs,"NIFTI-1+") == 0 )
5683 nim->nifti_type = NIFTI_FTYPE_NIFTI1_1 ;
5684 else if( strcmp(rhs,"NIFTI-1") == 0 )
5685 nim->nifti_type = NIFTI_FTYPE_NIFTI1_2 ;
5686 else if( strcmp(rhs,"NIFTI-1A") == 0 )
5687 nim->nifti_type = NIFTI_FTYPE_ASCII ;
5688 }
5689 else if( strcmp(lhs,"header_filename") == 0 ){
5690 nim->fname = nifti_strdup(rhs) ;
5691 }
5692 else if( strcmp(lhs,"image_filename") == 0 ){
5693 nim->iname = nifti_strdup(rhs) ;
5694 }
5695 else if( strcmp(lhs,"sto_xyz_matrix") == 0 ){
5696 sscanf( rhs , "%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f" ,
5697 &(nim->sto_xyz.m[0][0]) , &(nim->sto_xyz.m[0][1]) ,
5698 &(nim->sto_xyz.m[0][2]) , &(nim->sto_xyz.m[0][3]) ,
5699 &(nim->sto_xyz.m[1][0]) , &(nim->sto_xyz.m[1][1]) ,
5700 &(nim->sto_xyz.m[1][2]) , &(nim->sto_xyz.m[1][3]) ,
5701 &(nim->sto_xyz.m[2][0]) , &(nim->sto_xyz.m[2][1]) ,
5702 &(nim->sto_xyz.m[2][2]) , &(nim->sto_xyz.m[2][3]) ,
5703 &(nim->sto_xyz.m[3][0]) , &(nim->sto_xyz.m[3][1]) ,
5704 &(nim->sto_xyz.m[3][2]) , &(nim->sto_xyz.m[3][3]) ) ;
5705 }
5706 else if( strcmp(lhs,"byteorder") == 0 ){
5707 if( strcmp(rhs,"MSB_FIRST") == 0 ) nim->byteorder = MSB_FIRST ;
5708 if( strcmp(rhs,"LSB_FIRST") == 0 ) nim->byteorder = LSB_FIRST ;
5709 }
5710 else QQNUM(image_offset,iname_offset) ;
5711 else QNUM(datatype) ;
5712 else QNUM(ndim) ;
5713 else QNUM(nx) ;
5714 else QNUM(ny) ;
5715 else QNUM(nz) ;
5716 else QNUM(nt) ;
5717 else QNUM(nu) ;
5718 else QNUM(nv) ;
5719 else QNUM(nw) ;
5720 else QNUM(dx) ;
5721 else QNUM(dy) ;
5722 else QNUM(dz) ;
5723 else QNUM(dt) ;
5724 else QNUM(du) ;
5725 else QNUM(dv) ;
5726 else QNUM(dw) ;
5727 else QNUM(cal_min) ;
5728 else QNUM(cal_max) ;
5729 else QNUM(scl_slope) ;
5730 else QNUM(scl_inter) ;
5731 else QNUM(intent_code) ;
5732 else QNUM(intent_p1) ;
5733 else QNUM(intent_p2) ;
5734 else QNUM(intent_p3) ;
5735 else QSTR(intent_name,15) ;
5736 else QNUM(toffset) ;
5737 else QNUM(xyz_units) ;
5738 else QNUM(time_units) ;
5739 else QSTR(descrip,79) ;
5740 else QSTR(aux_file,23) ;
5741 else QNUM(qform_code) ;
5742 else QNUM(quatern_b) ;
5743 else QNUM(quatern_c) ;
5744 else QNUM(quatern_d) ;
5745 else QNUM(qoffset_x) ;
5746 else QNUM(qoffset_y) ;
5747 else QNUM(qoffset_z) ;
5748 else QNUM(qfac) ;
5749 else QNUM(sform_code) ;
5750 else QNUM(freq_dim) ;
5751 else QNUM(phase_dim) ;
5752 else QNUM(slice_dim) ;
5753 else QNUM(slice_code) ;
5754 else QNUM(slice_start) ;
5755 else QNUM(slice_end) ;
5756 else QNUM(slice_duration) ;
5757 else QNUM(num_ext) ;
5758
5759 } /* end of while loop */
5760
5761 if( bytes_read ) *bytes_read = spos+1; /* "process" last '\n' */
5762
5763 /* do miscellaneous checking and cleanup */
5764
5765 if( nim->ndim <= 0 ){ nifti_image_free(nim); return NULL; } /* bad! */
5766
5767 nifti_datatype_sizes( nim->datatype, &(nim->nbyper), &(nim->swapsize) );
5768 if( nim->nbyper == 0 ){ nifti_image_free(nim); return NULL; } /* bad! */
5769
5770 nim->dim[0] = nim->ndim ;
5771 nim->dim[1] = nim->nx ; nim->pixdim[1] = nim->dx ;
5772 nim->dim[2] = nim->ny ; nim->pixdim[2] = nim->dy ;
5773 nim->dim[3] = nim->nz ; nim->pixdim[3] = nim->dz ;
5774 nim->dim[4] = nim->nt ; nim->pixdim[4] = nim->dt ;
5775 nim->dim[5] = nim->nu ; nim->pixdim[5] = nim->du ;
5776 nim->dim[6] = nim->nv ; nim->pixdim[6] = nim->dv ;
5777 nim->dim[7] = nim->nw ; nim->pixdim[7] = nim->dw ;
5778
5779 nim->nvox = nim->nx * nim->ny * nim->nz
5780 * nim->nt * nim->nu * nim->nv * nim->nw ;
5781
5782 if( nim->qform_code > 0 )
5783 nim->qto_xyz = nifti_quatern_to_mat44(
5784 nim->quatern_b, nim->quatern_c, nim->quatern_d,
5785 nim->qoffset_x, nim->qoffset_y, nim->qoffset_z,
5786 nim->dx , nim->dy , nim->dz ,
5787 nim->qfac ) ;
5788 else
5789 nim->qto_xyz = nifti_quatern_to_mat44(
5790 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
5791 nim->dx , nim->dy , nim->dz , 0.0 ) ;
5792
5793
5794 nim->qto_ijk = nifti_mat44_inverse( nim->qto_xyz ) ;
5795
5796 if( nim->sform_code > 0 )
5797 nim->sto_ijk = nifti_mat44_inverse( nim->sto_xyz ) ;
5798
5799 return nim ;
5800 }
5801
5802
5803 /*---------------------------------------------------------------------------*/
5804 /*! validate the nifti_image
5805
5806 \return 1 if the structure seems valid, otherwise 0
5807
5808 \sa nifti_nim_has_valid_dims, nifti_hdr_looks_good
5809 *//*-------------------------------------------------------------------------*/
nifti_nim_is_valid(nifti_image * nim,int complain)5810 int nifti_nim_is_valid(nifti_image * nim, int complain)
5811 {
5812 int errs = 0;
5813
5814 if( !nim ){
5815 fprintf(stderr,"** is_valid_nim: nim is NULL\n");
5816 return 0;
5817 }
5818
5819 if( g_opts.debug > 2 ) fprintf(stderr,"-d nim_is_valid check...\n");
5820
5821 /**- check that dim[] matches the individual values ndim, nx, ny, ... */
5822 if( ! nifti_nim_has_valid_dims(nim,complain) ){
5823 if( !complain ) return 0;
5824 errs++;
5825 }
5826
5827 /* might check nbyper, pixdim, q/sforms, swapsize, nifti_type, ... */
5828
5829 /**- be explicit in return of 0 or 1 */
5830 if( errs > 0 ) return 0;
5831 else return 1;
5832 }
5833
5834 /*---------------------------------------------------------------------------*/
5835 /*! validate nifti dimensions
5836
5837 \return 1 if valid, 0 if not
5838
5839 \sa nifti_nim_is_valid, nifti_hdr_looks_good
5840
5841 rely on dim[] as the master
5842 *//*-------------------------------------------------------------------------*/
nifti_nim_has_valid_dims(nifti_image * nim,int complain)5843 int nifti_nim_has_valid_dims(nifti_image * nim, int complain)
5844 {
5845 int c, prod, errs = 0;
5846
5847 /**- start with dim[0]: failure here is considered terminal */
5848 if( nim->dim[0] <= 0 || nim->dim[0] > 7 ){
5849 errs++;
5850 if( complain )
5851 fprintf(stderr,"** NVd: dim[0] (%d) out of range [1,7]\n",nim->dim[0]);
5852 return 0;
5853 }
5854
5855 /**- check whether ndim equals dim[0] */
5856 if( nim->ndim != nim->dim[0] ){
5857 errs++;
5858 if( ! complain ) return 0;
5859 fprintf(stderr,"** NVd: ndim != dim[0] (%d,%d)\n",nim->ndim,nim->dim[0]);
5860 }
5861
5862 /**- compare each dim[i] to the proper nx, ny, ... */
5863 if( ( (nim->dim[0] >= 1) && (nim->dim[1] != nim->nx) ) ||
5864 ( (nim->dim[0] >= 2) && (nim->dim[2] != nim->ny) ) ||
5865 ( (nim->dim[0] >= 3) && (nim->dim[3] != nim->nz) ) ||
5866 ( (nim->dim[0] >= 4) && (nim->dim[4] != nim->nt) ) ||
5867 ( (nim->dim[0] >= 5) && (nim->dim[5] != nim->nu) ) ||
5868 ( (nim->dim[0] >= 6) && (nim->dim[6] != nim->nv) ) ||
5869 ( (nim->dim[0] >= 7) && (nim->dim[7] != nim->nw) ) ){
5870 errs++;
5871 if( !complain ) return 0;
5872 fprintf(stderr,"** NVd mismatch: dims = %d,%d,%d,%d,%d,%d,%d\n"
5873 " nxyz... = %d,%d,%d,%d,%d,%d,%d\n",
5874 nim->dim[1], nim->dim[2], nim->dim[3],
5875 nim->dim[4], nim->dim[5], nim->dim[6], nim->dim[7],
5876 nim->nx, nim->ny, nim->nz,
5877 nim->nt, nim->nu, nim->nv, nim->nw );
5878 }
5879
5880 /**- check the dimensions, and that their product matches nvox */
5881 prod = 1;
5882 for( c = 1; c <= nim->dim[0]; c++ ){
5883 if( nim->dim[c] > 0)
5884 prod *= nim->dim[c];
5885 else if( nim->dim[c] <= 0 ){
5886 if( !complain ) return 0;
5887 fprintf(stderr,"** NVd: dim[%d] (=%d) <= 0\n",c, nim->dim[c]);
5888 errs++;
5889 }
5890 }
5891 if( prod != nim->nvox ){
5892 if( ! complain ) return 0;
5893 fprintf(stderr,"** NVd: nvox does not match dimension product (%d, %d)\n",
5894 nim->nvox, prod);
5895 errs++;
5896 }
5897
5898 /**- if debug, warn about any remaining dim that is neither 0, nor 1 */
5899 /* (values in dims above dim[0] are undefined, as reminded by Cinly
5900 Ooi and Alle Meije Wink) 16 Nov 2005 [rickr] */
5901 if( g_opts.debug > 1 )
5902 for( c = nim->dim[0]+1; c <= 7; c++ )
5903 if( nim->dim[c] != 0 && nim->dim[c] != 1 )
5904 fprintf(stderr,"** NVd warning: dim[%d] = %d, but ndim = %d\n",
5905 c, nim->dim[c], nim->dim[0]);
5906
5907 if( g_opts.debug > 2 )
5908 fprintf(stderr,"-d nim_has_valid_dims check, errs = %d\n", errs);
5909
5910 /**- return invalid or valid */
5911 if( errs > 0 ) return 0;
5912 else return 1;
5913 }
5914
5915
5916 /*---------------------------------------------------------------------------*/
5917 /*! read a nifti image, collapsed across dimensions according to dims[8] <pre>
5918
5919 This function may be used to read parts of a nifti dataset, such as
5920 the time series for a single voxel, or perhaps a slice. It is similar
5921 to nifti_image_load(), though the passed 'data' parameter is used for
5922 returning the image, not nim->data.
5923
5924 \param nim given nifti_image struct, corresponding to the data file
5925 \param dims given list of dimensions (see below)
5926 \param data pointer to data pointer (if *data is NULL, data will be
5927 allocated, otherwise not)
5928
5929 Here, dims is an array of 8 ints, similar to nim->dim[8]. While dims[0]
5930 is unused at this point, the other indices specify which dimensions to
5931 collapse (and at which index), and which not to collapse. If dims[i] is
5932 set to -1, then that entire dimension will be read in, from index 0 to
5933 index (nim->dim[i] - 1). If dims[i] >= 0, then only that index will be
5934 read in (so dims[i] must also be < nim->dim[i]).
5935
5936 Example: given nim->dim[8] = { 4, 64, 64, 21, 80, 1, 1, 1 } (4-D dataset)
5937
5938 if dims[8] = { 0, 5, 4, 17, -1, -1, -1, -1 }
5939 -> read time series for voxel i,j,k = 5,4,17
5940
5941 if dims[8] = { 0, -1, -1, -1, 17, -1, -1, -1 }
5942 -> read single volume at time point 17
5943
5944 Example: given nim->dim[8] = { 6, 64, 64, 21, 80, 4, 3, 1 } (6-D dataset)
5945
5946 if dims[8] = { 0, 5, 4, 17, -1, 2, 1, 0 }
5947 -> read time series for the voxel i,j,k = 5,4,17, and dim 5,6 = 2,1
5948
5949 if dims[8] = { 0, 5, 4, -1, -1, 0, 0, 0 }
5950 -> read time series for slice at i,j = 5,4, and dim 5,6,7 = 0,0,0
5951 (note that dims[7] is not relevant, but must be 0 or -1)
5952
5953 If *data is NULL, then *data will be set as a pointer to new memory,
5954 allocated here for the resulting collapsed image data.
5955
5956 e.g. { int dims[8] = { 0, 5, 4, 17, -1, -1, -1, -1 };
5957 void * data = NULL;
5958 ret_val = nifti_read_collapsed_image(nim, dims, &data);
5959 if( ret_val > 0 ){
5960 process_time_series(data);
5961 if( data != NULL ) free(data);
5962 }
5963 }
5964
5965 NOTE: If *data is not NULL, then it will be assumed that it points to
5966 valid memory, sufficient to hold the results. This is done for
5967 speed and possibly repeated calls to this function.
5968
5969 e.g. { int dims[8] = { 0, -1, -1, -1, -1, -1, -1, -1 };
5970 void * data = NULL;
5971 for( zslice = 0; zslice < nzslices; zslice++ ){
5972 dims[3] = zslice;
5973 ret_val = nifti_read_collapsed_image(nim, dims, &data);
5974 if( ret_val > 0 ) process_slice(zslice, data);
5975 }
5976 if( data != NULL ) free(data);
5977 }
5978
5979 \return
5980 - the total number of bytes read, or < 0 on failure
5981 - the read and byte-swapped data, in 'data' </pre>
5982
5983 \sa nifti_image_read, nifti_image_free, nifti_image_read_bricks
5984 nifti_image_load
5985 *//*-------------------------------------------------------------------------*/
nifti_read_collapsed_image(nifti_image * nim,const int dims[8],void ** data)5986 int nifti_read_collapsed_image( nifti_image * nim, const int dims [8],
5987 void ** data )
5988 {
5989 znzFile fp;
5990 int pivots[8], prods[8], nprods; /* sizes are bounded by dims[], so 8 */
5991 int c, bytes;
5992
5993 /** - check pointers for sanity */
5994 if( !nim || !dims || !data ){
5995 fprintf(stderr,"** nifti_RCI: bad params %p, %p, %p\n",
5996 (void *)nim, (void *)dims, (void *)data);
5997 return -1;
5998 }
5999
6000 if( g_opts.debug > 2 ){
6001 fprintf(stderr,"-d read_collapsed_image:\n dims =");
6002 for(c = 0; c < 8; c++) fprintf(stderr," %3d", dims[c]);
6003 fprintf(stderr,"\n nim->dims =");
6004 for(c = 0; c < 8; c++) fprintf(stderr," %3d", nim->dim[c]);
6005 fputc('\n', stderr);
6006 }
6007
6008 /** - verify that dim[] makes sense */
6009 if( ! nifti_nim_is_valid(nim, g_opts.debug > 0) ){
6010 fprintf(stderr,"** invalid nim (file is '%s')\n", nim->fname );
6011 return -1;
6012 }
6013
6014 /** - verify that dims[] makes sense for this dataset */
6015 for( c = 1; c <= nim->dim[0]; c++ ){
6016 if( dims[c] >= nim->dim[c] ){
6017 fprintf(stderr,"** nifti_RCI: dims[%d] >= nim->dim[%d] (%d,%d)\n",
6018 c, c, dims[c], nim->dim[c]);
6019 return -1;
6020 }
6021 }
6022
6023 /** - prepare pivot list - pivots are fixed indices */
6024 if( make_pivot_list(nim, dims, pivots, prods, &nprods) < 0 ) return -1;
6025
6026 bytes = rci_alloc_mem(data, prods, nprods, nim->nbyper);
6027 if( bytes < 0 ) return -1;
6028
6029 /** - open the image file for reading at the appropriate offset */
6030 fp = nifti_image_load_prep( nim );
6031 if( ! fp ){ free(*data); *data = NULL; return -1; } /* failure */
6032
6033 /** - call the recursive reading function, passing nim, the pivot info,
6034 location to store memory, and file pointer and position */
6035 c = rci_read_data(nim, pivots,prods,nprods,dims,
6036 (char *)*data, fp, znztell(fp));
6037 znzclose(fp); /* in any case, close the file */
6038 if( c < 0 ){ free(*data); *data = NULL; return -1; } /* failure */
6039
6040 if( g_opts.debug > 1 )
6041 fprintf(stderr,"+d read %d bytes of collapsed image from %s\n",
6042 bytes, nim->fname);
6043
6044 return bytes;
6045 }
6046
6047
6048 /* read the data from the file pointed to by fp
6049
6050 - this a recursive function, so start with the base case
6051 - data is now (char *) for easy incrementing
6052
6053 return 0 on success, < 0 on failure
6054 */
rci_read_data(nifti_image * nim,int * pivots,int * prods,int nprods,const int dims[],char * data,znzFile fp,int base_offset)6055 static int rci_read_data(nifti_image * nim, int * pivots, int * prods,
6056 int nprods, const int dims[], char * data, znzFile fp, int base_offset)
6057 {
6058 int c, sublen, offset, read_size;
6059
6060 /* bad check first - base_offset may not have been checked */
6061 if( base_offset < 0 || nprods <= 0 ){
6062 fprintf(stderr,"** rci_read_data, bad params, %d,%d\n",
6063 nprods, base_offset);
6064 return -1;
6065 }
6066
6067 /* base case: actually read the data */
6068 if( nprods == 1 ){
6069 int nread, bytes;
6070
6071 /* make sure things look good here */
6072 if( *pivots != 0 ){
6073 fprintf(stderr,"** rciRD: final pivot == %d!\n", *pivots);
6074 return -1;
6075 }
6076
6077 /* so just seek and read (prods[0] * nbyper) bytes from the file */
6078 znzseek(fp, base_offset, SEEK_SET);
6079 bytes = prods[0] * nim->nbyper;
6080 nread = nifti_read_buffer(fp, data, bytes, nim);
6081 if( nread != bytes ){
6082 fprintf(stderr,"** rciRD: read only %d of %d bytes from '%s'\n",
6083 nread, bytes, nim->fname);
6084 return -1;
6085 } else if( g_opts.debug > 3 )
6086 fprintf(stderr,"+d successful read of %d bytes at offset %d\n",
6087 bytes, base_offset);
6088
6089 return 0; /* done with base case - return success */
6090 }
6091
6092 /* not the base case, so do a set of reduced reads */
6093
6094 /* compute size of sub-brick: all dimensions below pivot */
6095 for( c = 1, sublen = 1; c < *pivots; c++ ) sublen *= nim->dim[c];
6096
6097 /* compute number of values to read, i.e. remaining prods */
6098 for( c = 1, read_size = 1; c < nprods; c++ ) read_size *= prods[c];
6099 read_size *= nim->nbyper; /* and multiply by bytes per voxel */
6100
6101 /* now repeatedly compute offsets, and recursively read */
6102 for( c = 0; c < prods[0]; c++ ){
6103 /* offset is (c * sub-block size (including pivot dim)) */
6104 /* + (dims[] index into pivot sub-block) */
6105 /* the unneeded multiplication is to make this more clear */
6106 offset = c * sublen * nim->dim[*pivots] + sublen * dims[*pivots];
6107 offset *= nim->nbyper;
6108
6109 if( g_opts.debug > 3 )
6110 fprintf(stderr,"-d reading %d bytes, foff %d + %d, doff %d\n",
6111 read_size, base_offset, offset, c*read_size);
6112
6113 /* now read the next level down, adding this offset */
6114 if( rci_read_data(nim, pivots+1, prods+1, nprods-1, dims,
6115 data + c * read_size, fp, base_offset + offset) < 0 )
6116 return -1;
6117 }
6118
6119 return 0;
6120 }
6121
6122
6123 /* allocate memory for all collapsed image data
6124
6125 If *data is already set, do not allocate, but still calculate
6126 size for debug report.
6127
6128 return total size on success, and < 0 on failure
6129 */
rci_alloc_mem(void ** data,int prods[8],int nprods,int nbyper)6130 static int rci_alloc_mem(void ** data, int prods[8], int nprods, int nbyper )
6131 {
6132 int size, index;
6133
6134 if( nbyper < 0 || nprods < 1 || nprods > 8 ){
6135 fprintf(stderr,"** rci_am: bad params, %d, %d\n", nbyper, nprods);
6136 return -1;
6137 }
6138
6139 for( index = 0, size = 1; index < nprods; index++ )
6140 size *= prods[index];
6141
6142 size *= nbyper;
6143
6144 if( ! *data ){ /* then allocate what is needed */
6145 if( g_opts.debug > 1 )
6146 fprintf(stderr,"+d alloc %d (= %d x %d) bytes for collapsed image\n",
6147 size, size/nbyper, nbyper);
6148
6149 *data = malloc(size); /* actually allocate the memory */
6150 if( ! *data ){
6151 fprintf(stderr,"** rci_am: failed to alloc %d bytes for data\n", size);
6152 return -1;
6153 }
6154 } else if( g_opts.debug > 1 )
6155 fprintf(stderr,"-d rci_am: *data already set, need %d (%d x %d) bytes\n",
6156 size, size/nbyper, nbyper);
6157
6158 return size;
6159 }
6160
6161
6162 /* prepare a pivot list for reading
6163
6164 The pivot points are the indices into dims where the calling function
6165 wants to collapse a dimension. The last pivot should always be zero
6166 (note that we have space for that in the lists).
6167 */
make_pivot_list(nifti_image * nim,const int dims[],int pivots[],int prods[],int * nprods)6168 static int make_pivot_list(nifti_image * nim, const int dims[], int pivots[],
6169 int prods[], int * nprods )
6170 {
6171 int len, index;
6172
6173 len = 0;
6174 index = nim->dim[0];
6175 while( index > 0 ){
6176 prods[len] = 1;
6177 while( index > 0 && (nim->dim[index] == 1 || dims[index] == -1) ){
6178 prods[len] *= nim->dim[index];
6179 index--;
6180 }
6181 pivots[len] = index;
6182 len++;
6183 index--; /* fine, let it drop out at -1 */
6184 }
6185
6186 /* make sure to include 0 as a pivot (instead of just 1, if it is) */
6187 if( pivots[len-1] != 0 ){
6188 pivots[len] = 0;
6189 prods[len] = 1;
6190 len++;
6191 }
6192
6193 *nprods = len;
6194
6195 if( g_opts.debug > 2 ){
6196 fprintf(stderr,"+d pivot list created, pivots :");
6197 for(index = 0; index < len; index++) fprintf(stderr," %d", pivots[index]);
6198 fprintf(stderr,", prods :");
6199 for(index = 0; index < len; index++) fprintf(stderr," %d", prods[index]);
6200 fputc('\n',stderr);
6201 }
6202
6203 return 0;
6204 }
6205
6206
6207 #undef ISEND
6208 #define ISEND(c) ( (c)==']' || (c)=='}' || (c)=='\0' )
6209
6210 /*---------------------------------------------------------------------*/
6211 /*! Get an integer list in the range 0..(nvals-1), from the
6212 character string str. If we call the output pointer fred,
6213 then fred[0] = number of integers in the list (> 0), and
6214 fred[i] = i-th integer in the list for i=1..fred[0].
6215 If on return, fred == NULL or fred[0] == 0, then something is
6216 wrong, and the caller must deal with that.
6217
6218 Syntax of input string:
6219 - initial '{' or '[' is skipped, if present
6220 - ends when '}' or ']' or end of string is found
6221 - contains entries separated by commas
6222 - entries have one of these forms:
6223 - a single number
6224 - a dollar sign '$', which means nvals-1
6225 - a sequence of consecutive numbers in the form "a..b" or
6226 "a-b", where "a" and "b" are single numbers (or '$')
6227 - a sequence of evenly spaced numbers in the form
6228 "a..b(c)" or "a-b(c)", where "c" encodes the step
6229 - Example: "[2,7..4,3..9(2)]" decodes to the list
6230 2 7 6 5 4 3 5 7 9
6231 - entries should be in the range 0..nvals-1
6232
6233 (borrowed, with permission, from thd_intlist.c)
6234 *//*-------------------------------------------------------------------*/
nifti_get_intlist(int nvals,const char * str)6235 int * nifti_get_intlist( int nvals , const char * str )
6236 {
6237 int *subv = NULL ;
6238 int ii , ipos , nout , slen ;
6239 int ibot,itop,istep , nused ;
6240 char *cpt ;
6241
6242 /* Meaningless input? */
6243 if( nvals < 1 ) return NULL ;
6244
6245 /* No selection list? */
6246 if( str == NULL || str[0] == '\0' ) return NULL ;
6247
6248 /* skip initial '[' or '{' */
6249 subv = (int *) malloc( sizeof(int) * 2 ) ;
6250 subv[0] = nout = 0 ;
6251
6252 ipos = 0 ;
6253 if( str[ipos] == '[' || str[ipos] == '{' ) ipos++ ;
6254
6255 if( g_opts.debug > 1 )
6256 fprintf(stderr,"-d making int_list (vals = %d) from '%s'\n", nvals, str);
6257
6258 /**- for each sub-selector until end of input... */
6259
6260 slen = strlen(str) ;
6261 while( ipos < slen && !ISEND(str[ipos]) ){
6262
6263 while( isspace(str[ipos]) ) ipos++ ; /* skip blanks */
6264 if( ISEND(str[ipos]) ) break ; /* done */
6265
6266 /**- get starting value */
6267
6268 if( str[ipos] == '$' ){ /* special case */
6269 ibot = nvals-1 ; ipos++ ;
6270 } else { /* decode an integer */
6271 ibot = strtol( str+ipos , &cpt , 10 ) ;
6272 if( ibot < 0 ){
6273 fprintf(stderr,"** ERROR: list index %d is out of range 0..%d\n",
6274 ibot,nvals-1) ;
6275 free(subv) ; return NULL ;
6276 }
6277 if( ibot >= nvals ){
6278 fprintf(stderr,"** ERROR: list index %d is out of range 0..%d\n",
6279 ibot,nvals-1) ;
6280 free(subv) ; return NULL ;
6281 }
6282 nused = (cpt-(str+ipos)) ;
6283 if( ibot == 0 && nused == 0 ){
6284 fprintf(stderr,"** ERROR: list syntax error '%s'\n",str+ipos) ;
6285 free(subv) ; return NULL ;
6286 }
6287 ipos += nused ;
6288 }
6289
6290 while( isspace(str[ipos]) ) ipos++ ; /* skip blanks */
6291
6292 /**- if that's it for this sub-selector, add one value to list */
6293
6294 if( str[ipos] == ',' || ISEND(str[ipos]) ){
6295 nout++ ;
6296 subv = (int *) realloc( (char *)subv , sizeof(int) * (nout+1) ) ;
6297 subv[0] = nout ;
6298 subv[nout] = ibot ;
6299 if( ISEND(str[ipos]) ) break ; /* done */
6300 ipos++ ; continue ; /* re-start loop at next sub-selector */
6301 }
6302
6303 /**- otherwise, must have '..' or '-' as next inputs */
6304
6305 if( str[ipos] == '-' ){
6306 ipos++ ;
6307 } else if( str[ipos] == '.' && str[ipos+1] == '.' ){
6308 ipos++ ; ipos++ ;
6309 } else {
6310 fprintf(stderr,"** ERROR: index list syntax is bad: '%s'\n",
6311 str+ipos) ;
6312 free(subv) ; return NULL ;
6313 }
6314
6315 /**- get ending value for loop now */
6316
6317 if( str[ipos] == '$' ){ /* special case */
6318 itop = nvals-1 ; ipos++ ;
6319 } else { /* decode an integer */
6320 itop = strtol( str+ipos , &cpt , 10 ) ;
6321 if( itop < 0 ){
6322 fprintf(stderr,"** ERROR: index %d is out of range 0..%d\n",
6323 itop,nvals-1) ;
6324 free(subv) ; return NULL ;
6325 }
6326 if( itop >= nvals ){
6327 fprintf(stderr,"** ERROR: index %d is out of range 0..%d\n",
6328 itop,nvals-1) ;
6329 free(subv) ; return NULL ;
6330 }
6331 nused = (cpt-(str+ipos)) ;
6332 if( itop == 0 && nused == 0 ){
6333 fprintf(stderr,"** ERROR: index list syntax error '%s'\n",str+ipos) ;
6334 free(subv) ; return NULL ;
6335 }
6336 ipos += nused ;
6337 }
6338
6339 /**- set default loop step */
6340
6341 istep = (ibot <= itop) ? 1 : -1 ;
6342
6343 while( isspace(str[ipos]) ) ipos++ ; /* skip blanks */
6344
6345 /**- check if we have a non-default loop step */
6346
6347 if( str[ipos] == '(' ){ /* decode an integer */
6348 ipos++ ;
6349 istep = strtol( str+ipos , &cpt , 10 ) ;
6350 if( istep == 0 ){
6351 fprintf(stderr,"** ERROR: index loop step is 0!\n") ;
6352 free(subv) ; return NULL ;
6353 }
6354 nused = (cpt-(str+ipos)) ;
6355 ipos += nused ;
6356 if( str[ipos] == ')' ) ipos++ ;
6357 if( (ibot-itop)*istep > 0 ){
6358 fprintf(stderr,"** WARNING: index list '%d..%d(%d)' means nothing\n",
6359 ibot,itop,istep ) ;
6360 }
6361 }
6362
6363 /**- add values to output */
6364
6365 for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
6366 nout++ ;
6367 subv = (int *) realloc( (char *)subv , sizeof(int) * (nout+1) ) ;
6368 subv[0] = nout ;
6369 subv[nout] = ii ;
6370 }
6371
6372 /**- check if we have a comma to skip over */
6373
6374 while( isspace(str[ipos]) ) ipos++ ; /* skip blanks */
6375 if( str[ipos] == ',' ) ipos++ ; /* skip commas */
6376
6377 } /* end of loop through selector string */
6378
6379 if( g_opts.debug > 1 ) {
6380 fprintf(stderr,"+d int_list (vals = %d): ", subv[0]);
6381 for( ii = 1; ii <= subv[0]; ii++ ) fprintf(stderr,"%d ", subv[ii]);
6382 fputc('\n',stderr);
6383 }
6384
6385 if( subv[0] == 0 ){ free(subv); subv = NULL; }
6386 return subv ;
6387 }
6388