1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME : copy_data
3 @DESCRIPTION: File containing routines to copy data when reshaping.
4 @METHOD :
5 @GLOBALS :
6 @CREATED : October 25, 1994 (Peter Neelin)
7 @MODIFIED :
8 * $Log: copy_data.c,v $
9 * Revision 6.8 2008-01-17 02:33:02 rotor
10 * * removed all rcsids
11 * * removed a bunch of ^L's that somehow crept in
12 * * removed old (and outdated) BUGS file
13 *
14 * Revision 6.7 2008/01/13 09:38:54 stever
15 * Avoid compiler warnings about functions and variables that are defined
16 * but not used. Remove some such functions and variables,
17 * conditionalize some, and move static declarations out of header files
18 * into C files.
19 *
20 * Revision 6.6 2008/01/12 19:08:15 stever
21 * Add __attribute__ ((unused)) to all rcsid variables.
22 *
23 * Revision 6.5 2006/05/19 00:35:58 bert
24 * Add config.h to several files that might need it
25 *
26 * Revision 6.4 2004/11/01 22:38:39 bert
27 * Eliminate all references to minc_def.h
28 *
29 * Revision 6.3 2001/04/17 18:40:24 neelin
30 * Modifications to work with NetCDF 3.x
31 * In particular, changed NC_LONG to NC_INT (and corresponding longs to ints).
32 * Changed NC_UNSPECIFIED to NC_NAT.
33 * A few fixes to the configure script.
34 *
35 * Revision 6.2 1999/10/19 14:45:28 neelin
36 * Fixed Log subsitutions for CVS
37 *
38 * Revision 6.1 1998/08/19 13:05:28 neelin
39 * Added code to free minmax buffer.
40 *
41 * Revision 6.0 1997/09/12 13:24:12 neelin
42 * Release of minc version 0.6
43 *
44 * Revision 5.0 1997/08/21 13:25:10 neelin
45 * Release of minc version 0.5
46 *
47 * Revision 4.0 1997/05/07 20:01:44 neelin
48 * Release of minc version 0.4
49 *
50 * Revision 3.1 1995/10/03 13:34:14 neelin
51 * Fixed bug in truncate_input_vectors - was not handling out-of-range
52 * start values properly.
53 *
54 * Revision 3.0 1995/05/15 19:32:36 neelin
55 * Release of minc version 0.3
56 *
57 * Revision 1.5 1995/03/20 13:32:03 neelin
58 * Fixed -normalize option.
59 *
60 * Revision 1.4 1994/12/02 09:08:56 neelin
61 * Moved nd_loop to proglib.
62 *
63 * Revision 1.3 94/11/23 11:46:38 neelin
64 * Handle image-min/max properly when using icv for normalization.
65 *
66 * Revision 1.2 94/11/22 08:45:11 neelin
67 * Fixed handling of normalization for number of image dimensions > 2.
68 * Added appropriate default values of image-max and image-min.
69 *
70 * Revision 1.1 94/11/02 16:21:09 neelin
71 * Initial revision
72 *
73 @COPYRIGHT :
74 Copyright 1994 Peter Neelin, McConnell Brain Imaging Centre,
75 Montreal Neurological Institute, McGill University.
76 Permission to use, copy, modify, and distribute this
77 software and its documentation for any purpose and without
78 fee is hereby granted, provided that the above copyright
79 notice appear in all copies. The author and McGill University
80 make no representations about the suitability of this
81 software for any purpose. It is provided "as is" without
82 express or implied warranty.
83 ---------------------------------------------------------------------------- */
84
85 #if HAVE_CONFIG_H
86 #include "config.h"
87 #endif
88
89 #include <stdlib.h>
90 #include <stdio.h>
91 #include <float.h>
92 #include <limits.h>
93 #include <string.h>
94 #include <math.h>
95 #include <minc.h>
96 #include <nd_loop.h>
97 #include "mincreshape.h"
98
99 #define ROUND( x ) ((long) ((x) + ( ((x) >= 0) ? 0.5 : (-0.5) ) ))
100
101 static void get_num_minmax_values(Reshape_info *reshape_info,
102 long *block_start, long *block_count,
103 long *num_min_values, long *num_max_values);
104 static void handle_normalization(Reshape_info *reshape_info,
105 long *block_start,
106 long *block_count,
107 double *minmax_buffer,
108 double *fillvalue);
109 static void get_block_min_and_max(Reshape_info *reshape_info,
110 long *block_start,
111 long *block_count,
112 double *minmax_buffer,
113 double *minimum,
114 double *maximum);
115 static void truncate_input_vectors(Reshape_info *reshape_info,
116 long *input_start,
117 long *input_count);
118 static void translate_output_to_input(Reshape_info *reshape_info,
119 long *output_start,
120 long *output_count,
121 long *input_start,
122 long *input_count);
123 static void translate_input_to_output(Reshape_info *reshape_info,
124 long *input_start,
125 long *input_count,
126 long *output_start,
127 long *output_count);
128 static void copy_the_chunk(Reshape_info *reshape_info,
129 long chunk_start[],
130 long chunk_count[],
131 void *chunk_data,
132 double fillvalue);
133 static void convert_value_from_double(double dvalue,
134 nc_type datatype, int is_signed,
135 void *ptr);
136
137
138
139 /* ----------------------------- MNI Header -----------------------------------
140 @NAME : copy_data
141 @INPUT : reshape_info - information for reshaping volume
142 @OUTPUT : (none)
143 @RETURNS : (none)
144 @DESCRIPTION: Copies data from one input volume to another, reorganizing
145 it according to the reshaping info.
146 @METHOD :
147 @GLOBALS :
148 @CALLS :
149 @CREATED : October 25, 1994 (Peter Neelin)
150 @MODIFIED :
151 ---------------------------------------------------------------------------- */
copy_data(Reshape_info * reshape_info)152 void copy_data(Reshape_info *reshape_info)
153 {
154 int idim, odim, out_ndims;
155 long block_begin[MAX_VAR_DIMS], block_end[MAX_VAR_DIMS];
156 long block_count[MAX_VAR_DIMS];
157 long block_cur_start[MAX_VAR_DIMS], block_cur_count[MAX_VAR_DIMS];
158 long chunk_begin[MAX_VAR_DIMS], chunk_end[MAX_VAR_DIMS];
159 long chunk_count[MAX_VAR_DIMS];
160 long chunk_cur_start[MAX_VAR_DIMS], chunk_cur_count[MAX_VAR_DIMS];
161 long total_size;
162 long num_min_values, num_max_values, num_values;
163 double fillvalue, *minmax_buffer;
164 void *chunk_data;
165
166 /* Get number of dimensions */
167 out_ndims = reshape_info->output_ndims;
168
169 /* Set up variables for looping through blocks */
170 for (odim=0; odim < out_ndims; odim++) {
171 idim = reshape_info->map_out_to_in[odim];
172 block_begin[odim] = 0;
173 block_end[odim] = ABS(reshape_info->input_count[idim]);
174 if (reshape_info->dim_used_in_block[odim])
175 block_count[odim] = ABS(reshape_info->input_count[idim]);
176 else
177 block_count[odim] = 1;
178 }
179
180 /* Figure out size of chunks and allocate space */
181 total_size = nctypelen(reshape_info->output_datatype);
182 for (odim=0; odim < out_ndims; odim++) {
183 total_size *= reshape_info->chunk_count[odim];
184 }
185 chunk_data = malloc(total_size);
186
187 /* Get enough space for image-min and max values for a block */
188 get_num_minmax_values(reshape_info, NULL, block_count,
189 &num_min_values, &num_max_values);
190 num_values = ((num_min_values > num_max_values) ?
191 num_min_values : num_max_values);
192 if (num_values > 0)
193 minmax_buffer = malloc(num_values * sizeof(double));
194 else
195 minmax_buffer = NULL;
196
197 /* Print log message */
198 if (reshape_info->verbose) {
199 (void) fprintf(stderr, "Copying chunks:");
200 (void) fflush(stderr);
201 }
202
203 /* Loop through blocks */
204
205 nd_begin_looping(block_begin, block_cur_start, out_ndims);
206 while (!nd_end_of_loop(block_cur_start, block_end, out_ndims)) {
207
208 /* Set up count for current block */
209 nd_update_current_count(block_cur_start, block_count, block_end,
210 block_cur_count, out_ndims);
211
212 /* Set up chunk begin, end and count */
213 for (odim=0; odim < out_ndims; odim++) {
214 chunk_begin[odim] = block_cur_start[odim];
215 chunk_end[odim] = chunk_begin[odim] + block_cur_count[odim];
216 chunk_count[odim] = reshape_info->chunk_count[odim];
217 }
218
219 /* Set up icv for normalization, set output image-max/min and
220 calculate pixel fill value to use for current block */
221 handle_normalization(reshape_info, block_cur_start, block_cur_count,
222 minmax_buffer, &fillvalue);
223
224 /* Loop through chunks */
225
226 nd_begin_looping(chunk_begin, chunk_cur_start, out_ndims);
227 while (!nd_end_of_loop(chunk_cur_start, chunk_end, out_ndims)) {
228
229 /* Set up count for current chunk */
230 nd_update_current_count(chunk_cur_start, chunk_count, chunk_end,
231 chunk_cur_count, out_ndims);
232
233 /* Print log message for chunk */
234 if (reshape_info->verbose) {
235 (void) fprintf(stderr, ".");
236 (void) fflush(stderr);
237 }
238
239 /* Copy the chunk */
240 copy_the_chunk(reshape_info,
241 chunk_cur_start, chunk_cur_count, chunk_data,
242 fillvalue);
243
244 /* Increment chunk loop count */
245 nd_increment_loop(chunk_cur_start, chunk_begin, chunk_count,
246 chunk_end, out_ndims);
247
248 }
249
250 /* Increment block loop count */
251 nd_increment_loop(block_cur_start, block_begin, block_count,
252 block_end, out_ndims);
253
254 }
255
256 /* Free the chunk space */
257 free(chunk_data);
258
259 /* Free minmax buffer */
260 if (minmax_buffer != NULL) {
261 free(minmax_buffer);
262 }
263
264 /* Print ending log message */
265 if (reshape_info->verbose) {
266 (void) fprintf(stderr, "Done.\n");
267 (void) fflush(stderr);
268 }
269
270
271 }
272
273 /* ----------------------------- MNI Header -----------------------------------
274 @NAME : get_num_minmax_values
275 @INPUT : reshape_info - information for reshaping volume
276 block_start - start for a block (or NULL)
277 block_count - count for a block
278 @OUTPUT : num_min_values - number of image-min values to be read
279 num_max_values - number of image-max values to be read
280 @RETURNS : (nothing)
281 @DESCRIPTION: Gets the number of image-min and image-max values that
282 correspond to a block. If block_start is NULL, then
283 it is assumed to translate to an input start of [0,0,...].
284 Note that only the true number of values for the specified
285 block is computed (specifying a hyperslab that goes beyond
286 file extents does not give a bigger number of values).
287 @METHOD :
288 @GLOBALS :
289 @CALLS :
290 @CREATED : October 25, 1994 (Peter Neelin)
291 @MODIFIED :
292 ---------------------------------------------------------------------------- */
get_num_minmax_values(Reshape_info * reshape_info,long * block_start,long * block_count,long * num_min_values,long * num_max_values)293 static void get_num_minmax_values(Reshape_info *reshape_info,
294 long *block_start, long *block_count,
295 long *num_min_values, long *num_max_values)
296 {
297 int iloop, idim, ndims;
298 int varid, inimgid;
299 long size;
300 long minmax_count[MAX_VAR_DIMS];
301 long input_block_start[MAX_VAR_DIMS];
302 long input_block_count[MAX_VAR_DIMS];
303 long *num_values;
304
305 /* Check for icv normalization */
306 if (reshape_info->do_icv_normalization) {
307 *num_min_values = 0;
308 *num_max_values = 0;
309 return;
310 }
311
312 /* Translate output block count to input count */
313 translate_output_to_input(reshape_info, block_start, block_count,
314 input_block_start, input_block_count);
315 if (block_start != NULL) {
316 truncate_input_vectors(reshape_info, input_block_start,
317 input_block_count);
318 }
319 inimgid = ncvarid(reshape_info->inmincid, MIimage);
320
321 /* Loop over image-min and image-max */
322
323 for (iloop=0; iloop < 2; iloop++) {
324
325 /* Get varid and pointer to return value */
326 ncopts = 0;
327 switch (iloop) {
328 case 0:
329 varid = ncvarid(reshape_info->inmincid, MIimagemin);
330 num_values = num_min_values; /* Pointer to long */
331 break;
332 case 1:
333 varid = ncvarid(reshape_info->inmincid, MIimagemax);
334 num_values = num_max_values; /* Pointer to long */
335 break;
336 }
337 ncopts = NCOPTS_DEFAULT;
338
339 /* Translate block count to min or max count and work out the
340 total number of values. */
341 size = 0;
342 if (varid != MI_ERROR) {
343 (void) ncvarinq(reshape_info->inmincid, varid, NULL, NULL,
344 &ndims, NULL, NULL);
345 (void) mitranslate_coords(reshape_info->inmincid,
346 inimgid, input_block_count,
347 varid, minmax_count);
348 size = 1;
349 for (idim=0; idim < ndims; idim++) {
350 size *= minmax_count[idim];
351 }
352 }
353 *num_values = size;
354 }
355
356 }
357
358 /* ----------------------------- MNI Header -----------------------------------
359 @NAME : handle_normalization
360 @INPUT : reshape_info - information for reshaping volume
361 block_start - start of current block
362 block_count - count for current block
363 minmax_buffer - buffer space for getting min and max
364 values
365 @OUTPUT : fillvalue - pixel fill value to use for this block
366 @RETURNS : (none)
367 @DESCRIPTION: Sets up icv for normalization to ensure that block is
368 internally normalized. Output image-max and min are set.
369 The appropriate pixel fill value is calculated for this
370 min and max (applies to the whole block).
371 @METHOD :
372 @GLOBALS :
373 @CALLS :
374 @CREATED : October 25, 1994 (Peter Neelin)
375 @MODIFIED :
376 ---------------------------------------------------------------------------- */
handle_normalization(Reshape_info * reshape_info,long * block_start,long * block_count,double * minmax_buffer,double * fillvalue)377 static void handle_normalization(Reshape_info *reshape_info,
378 long *block_start,
379 long *block_count,
380 double *minmax_buffer,
381 double *fillvalue)
382 {
383 int iloop;
384 int inmincid, inimgid, varid, icvid;
385 long minmax_start[MAX_VAR_DIMS];
386 double minimum, maximum, *extreme, valid_min, valid_max, denom;
387 char *varname;
388
389 /* Get input minc id, image id and icv id*/
390 inmincid = reshape_info->inmincid;
391 inimgid = ncvarid(inmincid, MIimage);
392 icvid = reshape_info->icvid;
393
394 /* Get input min and max for block */
395 get_block_min_and_max(reshape_info, block_start, block_count,
396 minmax_buffer, &minimum, &maximum);
397
398 /* Modify the icv if necessary */
399 if (reshape_info->do_block_normalization) {
400 (void) miicv_detach(icvid);
401 (void) miicv_setdbl(icvid, MI_ICV_IMAGE_MIN, minimum);
402 (void) miicv_setdbl(icvid, MI_ICV_IMAGE_MAX, maximum);
403 (void) miicv_setint(icvid, MI_ICV_USER_NORM, TRUE);
404 (void) miicv_setint(icvid, MI_ICV_DO_NORM, TRUE);
405 (void) miicv_attach(icvid, inmincid, inimgid);
406 }
407
408 /* Save the image max and min for the block */
409 for (iloop=0; iloop < 2; iloop++) {
410
411 /* Get varid and pointer to min or max value */
412 switch (iloop) {
413 case 0:
414 varname = MIimagemin;
415 extreme = &minimum;
416 break;
417 case 1:
418 varname = MIimagemax;
419 extreme = &maximum;
420 break;
421 }
422
423 /* Save the value */
424 ncopts = 0;
425 varid = ncvarid(reshape_info->outmincid, varname);
426 ncopts = NCOPTS_DEFAULT;
427 if (varid != MI_ERROR) {
428 (void) mitranslate_coords(reshape_info->outmincid,
429 reshape_info->outimgid, block_start,
430 varid, minmax_start);
431 (void) mivarput1(reshape_info->outmincid, varid,
432 minmax_start, NC_DOUBLE, NULL, extreme);
433 }
434 }
435
436 /* Calculate the pixel fill value */
437 *fillvalue = ((reshape_info->fillvalue == NOFILL) ? 0.0 :
438 reshape_info->fillvalue);
439 if ((reshape_info->output_datatype != NC_FLOAT) &&
440 (reshape_info->output_datatype != NC_DOUBLE) &&
441 (*fillvalue != FILL)) {
442 (void) miicv_inqdbl(icvid, MI_ICV_VALID_MIN, &valid_min);
443 (void) miicv_inqdbl(icvid, MI_ICV_VALID_MAX, &valid_max);
444 denom = maximum - minimum;
445 if (denom == 0.0) {
446 *fillvalue = valid_min;
447 }
448 else {
449 *fillvalue = (*fillvalue - minimum) *
450 (valid_max - valid_min) / denom + valid_min;
451 }
452 }
453
454
455 }
456
457 /* ----------------------------- MNI Header -----------------------------------
458 @NAME : get_block_min_and_max
459 @INPUT : reshape_info - information for reshaping volume
460 block_start - start of current block
461 block_count - count for current block
462 minmax_buffer - buffer space for getting min and max
463 values
464 @OUTPUT : minimum - input minimum for block
465 maximum - input maximum for block
466 @RETURNS : (none)
467 @DESCRIPTION: Gets the min and max for the input file for a given output
468 block.
469 @METHOD :
470 @GLOBALS :
471 @CALLS :
472 @CREATED : October 25, 1994 (Peter Neelin)
473 @MODIFIED :
474 ---------------------------------------------------------------------------- */
get_block_min_and_max(Reshape_info * reshape_info,long * block_start,long * block_count,double * minmax_buffer,double * minimum,double * maximum)475 static void get_block_min_and_max(Reshape_info *reshape_info,
476 long *block_start,
477 long *block_count,
478 double *minmax_buffer,
479 double *minimum,
480 double *maximum)
481 {
482 int iloop;
483 long num_min_values, num_max_values, ivalue;
484 int inmincid, inimgid, varid, icvid;
485 long minmax_start[MAX_VAR_DIMS], minmax_count[MAX_VAR_DIMS];
486 long input_block_start[MAX_VAR_DIMS], input_block_count[MAX_VAR_DIMS];
487 double *extreme;
488 long num_values;
489 char *varname;
490 double sign, default_extreme;
491
492 /* Get input minc id, image id and icv id*/
493 inmincid = reshape_info->inmincid;
494 inimgid = ncvarid(inmincid, MIimage);
495 icvid = reshape_info->icvid;
496
497 /* Is the icv doing the normalization? */
498 if (reshape_info->do_icv_normalization) {
499 (void) miicv_inqdbl(icvid, MI_ICV_NORM_MIN, minimum);
500 (void) miicv_inqdbl(icvid, MI_ICV_NORM_MAX, maximum);
501 return;
502 }
503
504 /* Translate output block count to input count */
505 translate_output_to_input(reshape_info, block_start, block_count,
506 input_block_start, input_block_count);
507 truncate_input_vectors(reshape_info, input_block_start, input_block_count);
508
509 /* Get number of min and max values */
510 get_num_minmax_values(reshape_info, block_start, block_count,
511 &num_min_values, &num_max_values);
512
513 /* Loop over image-min and image-max getting block min and max */
514
515 for (iloop=0; iloop < 2; iloop++) {
516
517 /* Get varid and pointer to min or max value */
518 switch (iloop) {
519 case 0:
520 num_values = num_min_values;
521 varname = MIimagemin;
522 extreme = minimum;
523 sign = -1.0;
524 default_extreme = 0.0;
525 break;
526 case 1:
527 num_values = num_max_values;
528 varname = MIimagemax;
529 extreme = maximum;
530 sign = +1.0;
531 default_extreme = 1.0;
532 break;
533 }
534
535 /* Get values from file */
536 if (num_values > 0) {
537 varid = ncvarid(inmincid, varname);
538 (void) mitranslate_coords(inmincid,
539 inimgid, input_block_start,
540 varid, minmax_start);
541 (void) mitranslate_coords(inmincid,
542 inimgid, input_block_count,
543 varid, minmax_count);
544 (void) mivarget(reshape_info->inmincid, varid,
545 minmax_start, minmax_count,
546 NC_DOUBLE, NULL, minmax_buffer);
547 *extreme = minmax_buffer[0];
548 for (ivalue=1; ivalue < num_values; ivalue++) {
549 if ((minmax_buffer[ivalue] * sign) > (*extreme * sign))
550 *extreme = minmax_buffer[ivalue];
551 }
552 }
553 else {
554 *extreme = default_extreme;
555 }
556 if (reshape_info->need_fillvalue &&
557 (reshape_info->fillvalue == NOFILL) &&
558 (0.0 > (*extreme * sign)))
559 *extreme = 0.0;
560
561 }
562 }
563
564 /* ----------------------------- MNI Header -----------------------------------
565 @NAME : truncate_input_vectors
566 @INPUT : reshape_info - information for reshaping volume
567 input_start - start of input hyperslab (or NULL)
568 input_count - count for input hyperslab
569 @OUTPUT : input_start - start of input hyperslab (or NULL)
570 input_count - count for input hyperslab
571 @RETURNS : (nothing)
572 @DESCRIPTION: Input_start and input_count are truncated to specify a
573 legal hyperslab for the input file (if not specified,
574 input_start is assumed to be [0,0,...]).
575 @METHOD :
576 @GLOBALS :
577 @CALLS :
578 @CREATED : October 25, 1994 (Peter Neelin)
579 @MODIFIED :
580 ---------------------------------------------------------------------------- */
truncate_input_vectors(Reshape_info * reshape_info,long * input_start,long * input_count)581 static void truncate_input_vectors(Reshape_info *reshape_info,
582 long *input_start,
583 long *input_count)
584 {
585 int idim;
586 long first, last; /* last is actually last_index+1 */
587
588 /* Check for NULL vectors */
589 if (input_count == NULL) return;
590
591 /* Loop through input dimensions */
592 for (idim=0; idim < reshape_info->input_ndims; idim++) {
593 first = ( (input_start != NULL) ? input_start[idim] : 0 );
594 last = first + input_count[idim];
595 if (first < 0)
596 first = 0;
597 else if (first >= reshape_info->input_size[idim])
598 first = reshape_info->input_size[idim] - 1;
599 if (last < 0)
600 last = 0;
601 else if (last >= reshape_info->input_size[idim])
602 last = reshape_info->input_size[idim];
603 if (input_start != NULL)
604 input_start[idim] = first;
605 input_count[idim] = last - first;
606 }
607
608 }
609
610 /* ----------------------------- MNI Header -----------------------------------
611 @NAME : translate_output_to_input
612 @INPUT : reshape_info - information for reshaping volume
613 output_start - start of output hyperslab (or NULL)
614 output_count - count for output hyperslab
615 @OUTPUT : input_start - start of input hyperslab (or NULL)
616 input_count - count for input hyperslab
617 @RETURNS : (nothing)
618 @DESCRIPTION: Translates an output start and count to an input start and
619 count. If output_start or input_start are NULL, then only the
620 count is translated.
621 @METHOD :
622 @GLOBALS :
623 @CALLS :
624 @CREATED : October 25, 1994 (Peter Neelin)
625 @MODIFIED :
626 ---------------------------------------------------------------------------- */
translate_output_to_input(Reshape_info * reshape_info,long * output_start,long * output_count,long * input_start,long * input_count)627 static void translate_output_to_input(Reshape_info *reshape_info,
628 long *output_start,
629 long *output_count,
630 long *input_start,
631 long *input_count)
632 {
633 int idim, odim;
634
635 /* Check for NULL vectors */
636 if ((input_count == NULL) || (output_count == NULL)) return;
637
638 /* Loop through input dimensions */
639 for (idim=0; idim < reshape_info->input_ndims; idim++) {
640 odim = reshape_info->map_in_to_out[idim];
641 input_count[idim] = ((odim >= 0) ? output_count[odim] : 1);
642 if ((input_start != NULL) && (output_start != NULL)) {
643 input_start[idim] = reshape_info->input_start[idim];
644 if (odim >= 0) {
645 if (reshape_info->input_count[idim] > 0)
646 input_start[idim] += output_start[odim];
647 else
648 input_start[idim] -=
649 (output_start[odim] + output_count[odim] - 1);
650 }
651 }
652 }
653
654 }
655
656 /* ----------------------------- MNI Header -----------------------------------
657 @NAME : translate_input_to_output
658 @INPUT : reshape_info - information for reshaping volume
659 input_start - start of output hyperslab (or NULL)
660 input_count - count for output hyperslab
661 @OUTPUT : output_start - start of input hyperslab (or NULL)
662 output_count - count for input hyperslab
663 @RETURNS : (nothing)
664 @DESCRIPTION: Translates an input start and count to an output start and
665 count. If output_start or input_start are NULL, then only the
666 count is translated.
667 @METHOD :
668 @GLOBALS :
669 @CALLS :
670 @CREATED : October 25, 1994 (Peter Neelin)
671 @MODIFIED :
672 ---------------------------------------------------------------------------- */
translate_input_to_output(Reshape_info * reshape_info,long * input_start,long * input_count,long * output_start,long * output_count)673 static void translate_input_to_output(Reshape_info *reshape_info,
674 long *input_start,
675 long *input_count,
676 long *output_start,
677 long *output_count)
678 {
679 int idim, odim;
680
681 /* Check for NULL vectors */
682 if ((input_count == NULL) || (output_count == NULL)) return;
683
684 /* Loop through output dimensions */
685 for (odim=0; odim < reshape_info->output_ndims; odim++) {
686 idim = reshape_info->map_out_to_in[odim];
687 output_count[odim] = input_count[idim];
688 if ((input_start != NULL) && (output_start != NULL)) {
689 if (reshape_info->input_count[idim] > 0)
690 output_start[odim] = input_start[idim] -
691 reshape_info->input_start[idim];
692 else
693 output_start[odim] = reshape_info->input_start[idim] -
694 (input_start[idim] + input_count[idim] - 1);
695 }
696 }
697
698 }
699
700 /* ----------------------------- MNI Header -----------------------------------
701 @NAME : copy_the_chunk
702 @INPUT : reshape_info - information for reshaping volume
703 chunk_start - start of current block
704 chunk_count - count for current block
705 chunk_data - pointer to enough space for chunk
706 fillvalue - pixel value to zero volume, if necessary.
707 @OUTPUT : (none)
708 @RETURNS : (nothing)
709 @DESCRIPTION: Copies the chunk from the input file to the output file.
710 @METHOD :
711 @GLOBALS :
712 @CALLS :
713 @CREATED : October 25, 1994 (Peter Neelin)
714 @MODIFIED :
715 ---------------------------------------------------------------------------- */
copy_the_chunk(Reshape_info * reshape_info,long chunk_start[],long chunk_count[],void * chunk_data,double fillvalue)716 static void copy_the_chunk(Reshape_info *reshape_info,
717 long chunk_start[],
718 long chunk_count[],
719 void *chunk_data,
720 double fillvalue)
721 {
722 int idim, odim, in_ndims, out_ndims;
723 long input_start[MAX_VAR_DIMS], input_count[MAX_VAR_DIMS];
724 long output_start[MAX_VAR_DIMS], output_count[MAX_VAR_DIMS];
725 long input_imap[MAX_VAR_DIMS], output_imap[MAX_VAR_DIMS];
726 void *output_origin;
727 int datatype_size;
728 long total_size, ipix, first, last;
729 int zero_data, really_copy_the_data;
730 union {
731 char c; short s; long l; float f; double d;
732 } value_buffer;
733
734 /* Get number of dimensions */
735 out_ndims = reshape_info->output_ndims;
736 in_ndims = reshape_info->input_ndims;
737
738 /* Get size of output datatype */
739 datatype_size = nctypelen(reshape_info->output_datatype);
740
741 /* Create input start and count */
742 translate_output_to_input(reshape_info, chunk_start, chunk_count,
743 input_start, input_count);
744
745 /* Find out if we need to zero the volume and if we need to copy any
746 data */
747 zero_data = FALSE;
748 really_copy_the_data = TRUE;
749 total_size = 1;
750 for (idim=0; idim < in_ndims; idim++) {
751 first = input_start[idim];
752 last = input_start[idim] + input_count[idim] - 1;
753 if ((first < 0) || (last >= reshape_info->input_size[idim]))
754 zero_data = TRUE;
755 if ((last < 0) || (first >= reshape_info->input_size[idim]))
756 really_copy_the_data = FALSE;
757 total_size *= input_count[idim];
758 }
759
760 /* Make sure that input vectors are legal and translate them back
761 to output */
762 truncate_input_vectors(reshape_info, input_start, input_count);
763 translate_input_to_output(reshape_info, input_start, input_count,
764 output_start, output_count);
765
766 /* Write out zero data if needed */
767 if (zero_data) {
768 convert_value_from_double(fillvalue,
769 reshape_info->output_datatype,
770 reshape_info->output_is_signed,
771 &value_buffer);
772 for (ipix=0; ipix < total_size; ipix++) {
773 (void) memcpy((char *)chunk_data + ipix*datatype_size,
774 &value_buffer, datatype_size);
775 }
776 (void) ncvarput(reshape_info->outmincid, reshape_info->outimgid,
777 chunk_start, chunk_count, chunk_data);
778 }
779
780 /* Set up hypothetical imap variable for input */
781 for (idim=in_ndims-1; idim >= 0; idim--) {
782 input_imap[idim] = ((idim == in_ndims-1) ?
783 datatype_size :
784 input_imap[idim+1] * input_count[idim+1]);
785 }
786
787 /* Create output imap variable from input one (re-ordering dimensions and
788 flipping). Also work out the chunk origin (point to byte for output
789 [0,0,0...]). */
790 output_origin = chunk_data;
791 for (odim=0; odim < out_ndims; odim++) {
792 idim = reshape_info->map_out_to_in[odim];
793 if (reshape_info->input_count[idim] > 0) {
794 output_imap[odim] = input_imap[idim];
795 }
796 else {
797 output_imap[odim] = -input_imap[idim];
798 output_origin =
799 (void *) ((char *)output_origin -
800 (output_count[odim] - 1) * output_imap[odim]);
801 }
802 }
803
804 /* Should we really copy the data? */
805 if (really_copy_the_data) {
806
807 /* Read in the data */
808 (void) miicv_get(reshape_info->icvid, input_start, input_count,
809 chunk_data);
810
811 /* Write it out */
812 (void) ncvarputg(reshape_info->outmincid, reshape_info->outimgid,
813 output_start, output_count, NULL, output_imap,
814 output_origin);
815
816 }
817
818 }
819
820 /* ----------------------------- MNI Header -----------------------------------
821 @NAME : convert_value_from_double
822 @INPUT : dvalue - double value to convert
823 datatype - type of desired value
824 is_signed - TRUE if desired value is signed
825 @OUTPUT : ptr - pointer to converted value
826 @RETURNS : (nothing)
827 @DESCRIPTION: Converts a value from double to some other value.
828 @METHOD :
829 @GLOBALS :
830 @CALLS :
831 @CREATED : October 25, 1994 (Peter Neelin)
832 @MODIFIED :
833 ---------------------------------------------------------------------------- */
convert_value_from_double(double dvalue,nc_type datatype,int is_signed,void * ptr)834 static void convert_value_from_double(double dvalue,
835 nc_type datatype, int is_signed,
836 void *ptr)
837 {
838 switch (datatype) {
839 case NC_BYTE :
840 if (!is_signed) {
841 dvalue = MAX(0, dvalue);
842 dvalue = MIN(UCHAR_MAX, dvalue);
843 *((unsigned char *) ptr) = ROUND(dvalue);
844 }
845 else {
846 dvalue = MAX(SCHAR_MIN, dvalue);
847 dvalue = MIN(SCHAR_MAX, dvalue);
848 *((signed char *) ptr) = ROUND(dvalue);
849 }
850 break;
851 case NC_SHORT :
852 if (!is_signed) {
853 dvalue = MAX(0, dvalue);
854 dvalue = MIN(USHRT_MAX, dvalue);
855 *((unsigned short *) ptr) = ROUND(dvalue);
856 }
857 else {
858 dvalue = MAX(SHRT_MIN, dvalue);
859 dvalue = MIN(SHRT_MAX, dvalue);
860 *((signed short *) ptr) = ROUND(dvalue);
861 }
862 break;
863 case NC_INT :
864 if (!is_signed) {
865 dvalue = MAX(0, dvalue);
866 dvalue = MIN(UINT_MAX, dvalue);
867 *((unsigned int *) ptr) = ROUND(dvalue);
868 }
869 else {
870 dvalue = MAX(INT_MIN, dvalue);
871 dvalue = MIN(INT_MAX, dvalue);
872 *((signed int *) ptr) = ROUND(dvalue);
873 }
874 break;
875 case NC_FLOAT :
876 dvalue = MAX(-FLT_MAX,dvalue);
877 *((float *) ptr) = MIN(FLT_MAX,dvalue);
878 break;
879 case NC_DOUBLE :
880 *((double *) ptr) = dvalue;
881 break;
882 }
883 }
884