1 /* ----------------------------- MNI Header -----------------------------------
2 @NAME       :
3 @DESCRIPTION: Primitive" interface to minc files, using minctoraw and rawtominc programs
4               does not require linking with minc2 library.
5               Created during the days when minc2 didn't compile on 64bit architecture
6 @COPYRIGHT  :
7               Copyright 2007 Vladimir Fonov, McConnell Brain Imaging Centre,
8               Montreal Neurological Institute, McGill University.
9               Permission to use, copy, modify, and distribute this
10               software and its documentation for any purpose and without
11               fee is hereby granted, provided that the above copyright
12               notice appear in all copies.  The author and McGill University
13               make no representations about the suitability of this
14               software for any purpose.  It is provided "as is" without
15               express or implied warranty.
16 ---------------------------------------------------------------------------- */
17 #include <stdio.h>
18 #include <iostream>
19 #include <sstream>
20 #include <string.h>
21 //minc stuff
22 #include <math.h>
23 #include <limits.h>
24 #include "minc_1_rw.h"
25 
26 namespace minc
27 {
dim_info(int l,double sta,double spa,dimensions d,bool hd)28   dim_info::dim_info(int l, double sta,
29                      double spa,dimensions d,
30                      bool hd):
31       length(l),start(sta),step(spa),dim(d),have_dir_cos(hd)
32   {
33     switch(dim)
34     {
35       case dim_info::DIM_X:name=MIxspace;break;
36       case dim_info::DIM_Y:name=MIyspace;break;
37       case dim_info::DIM_Z:name=MIzspace;break;
38       case dim_info::DIM_TIME:name=MItime;break;
39       case dim_info::DIM_VEC:name=MIvector_dimension;break;
40       default: REPORT_ERROR("Unknown Dimension!");
41     }
42   }
43 
minc_1_base()44   minc_1_base::minc_1_base():
45     _dims(3,0),
46     _map_to_std(5,-1),
47     _mincid(MI_ERROR),
48     _imgid(MI_ERROR),
49     _icvid(MI_ERROR),
50     _cur(MAX_VAR_DIMS,0),
51     _slab(MAX_VAR_DIMS,1),
52     _last(false),
53     _datatype(MI_ORIGINAL_TYPE),
54     _io_datatype(MI_ORIGINAL_TYPE),
55     _slab_len(0),
56     _positive_directions(true),
57     _minc2(false)
58   {
59     _icvid=miicv_create();
60   }
61 
62   //! destructor, closes minc file
~minc_1_base()63   minc_1_base::~minc_1_base()
64   {
65     close();
66   }
67 
close(void)68   void minc_1_base::close(void)
69   {
70     if(_icvid!=MI_ERROR)
71     {
72       miicv_free(_icvid);
73       _icvid=MI_ERROR;
74     }
75     if(_mincid!=MI_ERROR) miclose(_mincid);
76     _mincid=MI_ERROR;
77   }
78 
79 
history(void) const80   std::string minc_1_base::history(void) const
81   {
82     nc_type datatype;
83     int att_length;
84 #ifndef WIN32
85     int op=ncopts;
86 #endif //WIN32
87     //ncopts=0;
88     if ((ncattinq(_mincid, NC_GLOBAL, MIhistory, &datatype,&att_length) == MI_ERROR) ||
89         (datatype != NC_CHAR))
90     {
91       //ncopts=op;
92       return "";
93     }
94     char* str = new char[att_length+1];
95     str[0] = '\0';
96     miattgetstr(_mincid, NC_GLOBAL, (char*)MIhistory, att_length,str);
97     //ncopts=op;
98     std::string r(str);
99     delete [] str;
100     return r;
101   }
102 
103   //code from mincinfo
var_number(void) const104   int minc_1_base::var_number(void) const
105   {
106     int nvars;
107     if(ncinquire(_mincid, NULL, &nvars, NULL, NULL)!=MI_ERROR)
108       return nvars;
109     return 0;
110   }
111 
var_name(int no) const112   std::string minc_1_base::var_name(int no) const
113   {
114     char name[MAX_NC_NAME];
115     if(ncvarinq(_mincid, no, name, NULL, NULL, NULL, NULL)!=MI_ERROR)
116       return name;
117   }
118 
att_number(const char * var_name) const119   int minc_1_base::att_number(const char *var_name) const
120   {
121     int varid;
122     int natts;
123     if (*var_name=='\0') {
124         varid = NC_GLOBAL;
125     } else {
126       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
127         return 0;
128     }
129     return att_number(varid);
130   }
131 
var_id(const char * var_name) const132   int minc_1_base::var_id(const char *var_name) const
133   {
134     return ncvarid(_mincid, var_name);
135   }
136 
var_length(const char * var_name) const137   long minc_1_base::var_length(const char *var_name) const
138   {
139     int varid=var_id(var_name);
140     if(varid!=MI_ERROR)
141       return var_length(varid);
142     else
143       return 0;
144   }
145 
var_length(int var_id) const146   long minc_1_base::var_length(int var_id) const
147   {
148     int vardims;
149 
150     if(ncvarinq(_mincid, var_id, NULL, NULL, &vardims, NULL, NULL)!=MI_ERROR)
151     {
152       if(vardims==0) return 1;
153       int *dims=new int[vardims];
154       if(ncvarinq(_mincid, var_id, NULL, NULL, NULL, dims, NULL)!=MI_ERROR)
155       {
156         long varlength=1;
157         if(ncdiminq(_mincid,dims[0],NULL,&varlength)!=MI_ERROR)
158         {
159           delete[] dims;
160           return varlength;
161         }
162         delete[] dims;
163         return 1;
164       } else {
165         delete[] dims;
166         return 1;
167       }
168     }
169     return 0;
170   }
171 
172 
173   //! get the number of attributes associated with variable
att_number(int var_no) const174   int minc_1_base::att_number(int var_no) const
175   {
176     int natts;
177     if(ncvarinq(_mincid, var_no, NULL, NULL, NULL, NULL, &natts)!=MI_ERROR)
178       return natts;
179     return 0;
180   }
181 
182 
att_name(const char * var_name,int no) const183   std::string minc_1_base::att_name(const char *var_name,int no) const
184   {
185     int varid;
186     int attid;
187     char name[MAX_NC_NAME];
188     if (*var_name=='\0')
189         varid = NC_GLOBAL;
190     else
191     {
192       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
193         return "";
194     }
195     return att_name(varid,no);
196   }
197 
att_name(int varid,int no) const198   std::string minc_1_base::att_name(int varid,int no) const
199   {
200     int attid;
201     char name[MAX_NC_NAME];
202     if(ncattname(_mincid, varid, no, name)==MI_ERROR)
203       return "";
204 
205     return name;
206   }
207 
att_value_string(const char * var_name,const char * att_name) const208   std::string minc_1_base::att_value_string(const char *var_name,const char *att_name) const
209   {
210     int varid;
211     int attid;
212     char name[MAX_NC_NAME];
213     if (*var_name=='\0')
214         varid = NC_GLOBAL;
215     else
216     {
217       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
218         return "";
219     }
220     return att_value_string(varid,att_name);
221   }
222 
att_value_string(int varid,const char * att_name) const223   std::string minc_1_base::att_value_string(int varid,const char *att_name) const
224   {
225     int att_length;
226     nc_type datatype;
227 
228     //TODO: make this handle other (double?) data types correctly
229     if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) ||
230         (datatype != NC_CHAR))
231     {
232       //ncopts=op;
233       return "";
234     }
235     char* str = new char[att_length+1];
236     str[0] = '\0';
237     miattgetstr(_mincid, varid, (char *)att_name, att_length, str);
238     //ncopts=op;
239     std::string r(str);
240     delete [] str;
241     return r;
242   }
243 
att_value_double(const char * var_name,const char * att_name) const244   std::vector<double> minc_1_base::att_value_double(const char *var_name,const char *att_name) const
245   {
246     int varid;
247     if (*var_name=='\0')
248         varid = NC_GLOBAL;
249     else
250     {
251       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
252         return std::vector<double>(0);
253     }
254     return att_value_double(varid,att_name);
255   }
256 
att_value_short(const char * var_name,const char * att_name) const257   std::vector<short> minc_1_base::att_value_short(const char *var_name,const char *att_name) const
258   {
259     int varid;
260     if (*var_name=='\0')
261       varid = NC_GLOBAL;
262     else
263     {
264       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
265         return std::vector<short>(0);
266     }
267     return att_value_short(varid,att_name);
268   }
269 
att_value_byte(const char * var_name,const char * att_name) const270   std::vector<unsigned char> minc_1_base::att_value_byte(const char *var_name,const char *att_name) const
271   {
272     int varid;
273     if (*var_name=='\0')
274       varid = NC_GLOBAL;
275     else
276     {
277       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
278         return std::vector<unsigned char>(0);
279     }
280     return att_value_byte(varid,att_name);
281   }
282 
283 
att_value_int(const char * var_name,const char * att_name) const284   std::vector<int> minc_1_base::att_value_int(const char *var_name,const char *att_name) const
285   {
286     int varid;
287     if (*var_name=='\0')
288       varid = NC_GLOBAL;
289     else
290     {
291       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
292         return std::vector<int>(0);
293     }
294     return att_value_int(varid,att_name);
295   }
296 
att_value_int(int varid,const char * att_name) const297   std::vector<int> minc_1_base::att_value_int(int varid,const char *att_name) const
298   {
299     int att_length;
300     nc_type datatype;
301 
302     if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) ||
303          (datatype != NC_INT))
304     {
305       //ncopts=op;
306       return std::vector<int>(0);
307     }
308     std::vector<int> r(att_length);
309     miattget(_mincid, varid, (char*)att_name, NC_INT, att_length,&r[0], NULL) ;
310     //ncopts=op;
311     return r;
312   }
313 
314 
att_value_double(int varid,const char * att_name) const315   std::vector<double> minc_1_base::att_value_double(int varid,const char *att_name) const
316   {
317     int att_length;
318     nc_type datatype;
319 
320     //TODO: make this handle other (double?) data types correctly
321     if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) ||
322         (datatype != NC_DOUBLE))
323     {
324       //ncopts=op;
325       return std::vector<double>(0);
326     }
327     std::vector<double> r(att_length);
328     miattget(_mincid, varid, (char*)att_name, NC_DOUBLE, att_length,&r[0], NULL) ;
329     //ncopts=op;
330     return r;
331   }
332 
att_value_short(int varid,const char * att_name) const333   std::vector<short> minc_1_base::att_value_short(int varid,const char *att_name) const
334   {
335     int att_length;
336     nc_type datatype;
337 
338     //TODO: make this handle other (double?) data types correctly
339     if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) ||
340          (datatype != NC_SHORT))
341     {
342       //ncopts=op;
343       return std::vector<short>(0);
344     }
345     std::vector<short> r(att_length);
346     miattget(_mincid, varid, (char*)att_name, NC_SHORT, att_length,&r[0], NULL) ;
347     //ncopts=op;
348     return r;
349   }
350 
att_value_byte(int varid,const char * att_name) const351   std::vector<unsigned char> minc_1_base::att_value_byte(int varid,const char *att_name) const
352   {
353     int att_length;
354     nc_type datatype;
355 
356     //TODO: make this handle other (double?) data types correctly
357     if ((ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR) ||
358          (datatype != NC_BYTE))
359     {
360       //ncopts=op;
361       return std::vector<unsigned char>(0);
362     }
363     std::vector<unsigned char> r(att_length);
364     miattget(_mincid, varid, (char*)att_name, NC_BYTE, att_length,&r[0], NULL) ;
365     //ncopts=op;
366     return r;
367   }
368 
369 
att_type(const char * var_name,const char * att_name) const370   nc_type minc_1_base::att_type(const char *var_name,const char *att_name) const
371   {
372     int varid;
373     if (*var_name=='\0')
374         varid = NC_GLOBAL;
375     else
376     {
377       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
378         return MI_ORIGINAL_TYPE;
379     }
380     return att_type(varid,att_name);
381   }
382 
att_type(int varid,const char * att_name) const383   nc_type minc_1_base::att_type(int varid,const char *att_name) const
384   {
385     int att_length;
386     nc_type datatype;
387 
388     if(ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR)
389       return MI_ORIGINAL_TYPE;
390     return datatype;
391   }
392 
393 
att_length(const char * var_name,const char * att_name) const394   int minc_1_base::att_length(const char *var_name,const char *att_name) const
395   {
396     int varid;
397     if (*var_name=='\0')
398         varid = NC_GLOBAL;
399     else
400     {
401       if((varid = ncvarid(_mincid, var_name))==MI_ERROR)
402         return 0;
403     }
404     return att_length(varid,att_name);
405   }
406 
att_length(int varid,const char * att_name) const407   int minc_1_base::att_length(int varid,const char *att_name) const
408   {
409     int att_length;
410     nc_type datatype;
411 
412     if(ncattinq(_mincid, varid, (char *)att_name, &datatype,&att_length) == MI_ERROR)
413       return 0;
414 
415     return att_length;
416   }
417 
418 
minc_1_reader(const minc_1_reader &)419   minc_1_reader::minc_1_reader(const minc_1_reader&):_metadate_only(false),_have_temp_file(false),_read_prepared(false)
420   {
421   }
422 
minc_1_reader()423   minc_1_reader::minc_1_reader():_metadate_only(false),_have_temp_file(false),_read_prepared(false)
424   {
425   }
426 
427 
428   //based on the code from mincextract
open(const char * path,bool positive_directions,bool metadate_only,bool rw)429   void minc_1_reader::open(const char *path,bool positive_directions/*=true*/,bool metadate_only/*=false*/,bool rw/*=false*/)
430   {
431 #ifndef WIN32
432     ncopts = 0;
433 #endif
434     _metadate_only=metadate_only;
435     _read_prepared=false;
436     int element_size;
437     int idim;
438     int nstart, ncount;
439     _positive_directions=positive_directions;
440     //ncopts = 0;
441 
442     // Open the file
443 
444     // Expand file
445     if(_metadate_only)
446     {
447       int created_tempfile;
448       char * tempfile = miexpand_file((char*)path, NULL, true, &created_tempfile);
449       if (tempfile == NULL) REPORT_ERROR("Error expanding minc file");
450       _tempfile=tempfile;
451       path=_tempfile.c_str();
452       _have_temp_file=created_tempfile;
453     }
454     _mincid = miopen((char*)path, rw?NC_WRITE:NC_NOWRITE);
455 
456     if(_mincid == MI_ERROR) REPORT_ERROR("Can't open minc file for reading!");
457 #ifdef MINC2
458    if (MI2_ISH5OBJ(_mincid)) {
459      _minc2 = true;
460    }
461 #endif
462 
463     /* Inquire about the image variable */
464     _imgid = ncvarid(_mincid, MIimage);
465     if(_imgid == MI_ERROR) REPORT_ERROR("Can't get Image ID");
466     ncvarinq(_mincid, _imgid, NULL, NULL, &_ndims, mdims, NULL);
467     //get image data type... not used for now
468     miget_datatype(_mincid, _imgid, &_datatype, &_is_signed);
469     //dir_cos.SetIdentity();
470     int dim_cnt=0;
471     miget_image_range(_mincid, _image_range);
472     //go through dimensions , calculating parameters for reshaping into ZYX array if needed
473     _info.resize(_ndims);
474     _world_matrix.resize(_ndims*4,0.0);
475     _voxel_matrix.resize(_ndims*4,0);
476     //_dir_cos.resize(_ndims*_ndims,0.0);
477 
478     for(int i=_ndims-1;i>=0;i--)
479     {
480       //_world_matrix[i*(_ndims+1)]=1.0;
481       //_voxel_matrix[i*(_ndims+1)]=1.0;
482       char dimname[MAX_NC_NAME];
483       long dimlength;
484       //get dimensions info
485       ncdiminq(_mincid, mdims[i], dimname, &dimlength);
486       _info[i].name=dimname;
487       _info[i].length=dimlength;
488       _info[i].have_dir_cos=false;
489       int axis=-1;
490       unsigned int sz=0;
491 
492       if(!strcmp(dimname,MIxspace))
493       {
494         _dims[0]=dimlength;axis=0;
495         _info[i].dim=dim_info::DIM_X;
496         _map_to_std[1]=i;
497       } else if(!strcmp(dimname,MIyspace)) {
498         _dims[1]=dimlength;axis=1;
499         _info[i].dim=dim_info::DIM_Y;
500         _map_to_std[2]=i;
501       } else if(!strcmp(dimname,MIzspace)) {
502         _dims[2]=dimlength;axis=2;
503         _info[i].dim=dim_info::DIM_Z;
504         _map_to_std[3]=i;
505       } else if(!strcmp(dimname,MIvector_dimension)) {
506          axis=-1;
507         _info[i].dim=dim_info::DIM_VEC;
508         _map_to_std[0]=i;
509       } else if(!strcmp(dimname,MItime)) {
510          axis=3;
511         _info[i].dim=dim_info::DIM_TIME;
512         _map_to_std[4]=i;
513       } else  {
514         REPORT_ERROR ("Unknown dimension");
515         _info[i].dim=dim_info::DIM_UNKNOWN;
516       }
517 
518       if(_info[i].dim!=dim_info::DIM_VEC)
519       {
520         //ncopts = 0;
521         int dimid = ncvarid(_mincid, dimname);
522         //ncopts = NC_VERBOSE | NC_FATAL;
523         if (dimid == MI_ERROR) continue;
524 
525         // Get dimension attributes
526         //ncopts = 0;
527         miattget1(_mincid, dimid, (char*)MIstep, NC_DOUBLE, &_info[i].step);
528         if(_info[i].step == 0.0)
529            _info[i].step = 1.0;
530         miattget1(_mincid, dimid, (char*)MIstart, NC_DOUBLE, &_info[i].start);
531 
532         if(_positive_directions && _info[i].step<0.0)
533         {
534           _info[i].start+=_info[i].step*(dimlength-1);
535           _info[i].step=-_info[i].step;
536         }
537 
538         if(miattget(_mincid, dimid, (char*)MIdirection_cosines, NC_DOUBLE, 3, &_info[i].dir_cos[0], NULL)!= MI_ERROR)
539         {
540           _info[i].have_dir_cos=true;
541 
542           /* Normalize the direction cosine */
543           double len=sqrt(_info[i].dir_cos[0]*_info[i].dir_cos[0]+
544                           _info[i].dir_cos[1]*_info[i].dir_cos[1]+
545                           _info[i].dir_cos[2]*_info[i].dir_cos[2]);
546 
547           if(len>1e-6 && fabs(len-1.0)>1e-6) //TODO: use some epsiolon here?
548           {
549             for(int a=0;a<3;a++)
550               _info[i].dir_cos[a]/=len;
551           }
552 
553         }
554         // fill voxel matrix
555         _voxel_matrix[i*4+axis]=1;
556       } else { //vectors don't have spatial component!
557         _info[i].start=0;
558         _info[i].step=0.0;
559         _info[i].dir_cos[0]=_info[i].dir_cos[1]=_info[i].dir_cos[2]=0.0;
560         _info[i].have_dir_cos=false;
561       }
562 
563       //fill world matrix
564       for(int a=0;a<3;a++)
565         _world_matrix[i*4+a]=_info[i].dir_cos[a]*_info[i].step;
566       if(axis==3) //time
567         _world_matrix[i*4+3]=_info[i].step;
568       else
569         _world_matrix[i*4+3]=0.0;
570     }
571     //ncopts = NC_VERBOSE | NC_FATAL;
572 
573     // now let's find out the slice dimensions
574     int idmax = ncvarid(_mincid, MIimagemax);
575     _slice_dimensions=0;
576     if(idmax != MI_ERROR)
577     {
578       int nmax_dims;
579       int mmax_dims[MAX_VAR_DIMS];
580       ncvarinq(_mincid, _imgid, NULL, NULL, &nmax_dims, mmax_dims, NULL);
581       if(nmax_dims>0)
582         _slice_dimensions=_ndims-nmax_dims;
583     }
584 
585     if(_slice_dimensions<=0)
586     {
587       if(_info[_ndims-1].dim==dim_info::DIM_VEC || _info[_ndims-1].dim==dim_info::DIM_TIME)
588         _slice_dimensions=std::min(_ndims,3);
589       else
590         _slice_dimensions=std::min(_ndims,2);
591     }
592     std::fill(_slab.begin(),_slab.end(),1);
593     _slab_len=1;
594     for(int i=0;i<_slice_dimensions;i++)
595     {
596       _slab[_ndims-i-1]=_info[_ndims-i-1].length;
597       _slab_len*=_info[_ndims-i-1].length;
598     }
599   }
600 
close(void)601   void minc_1_reader::close(void)
602   {
603     if(_have_temp_file)
604       remove(_tempfile.c_str());
605 
606     _have_temp_file=false;
607 
608     minc_1_base::close();
609   }
610 
~minc_1_reader()611   minc_1_reader::~minc_1_reader()
612   {
613     close();
614   }
615 
minc_1_writer()616   minc_1_writer::minc_1_writer():
617       _set_image_range(false),_set_slice_range(false),
618       _calc_min_max(true),_write_prepared(false)
619   {
620   }
621 
minc_1_writer(const minc_1_writer &)622   minc_1_writer::minc_1_writer(const minc_1_writer&):
623       _set_image_range(false),_set_slice_range(false),
624       _calc_min_max(true),_write_prepared(false)
625   {
626   }
627 
628 
open(const char * path,const minc_info & inf,int slice_dimensions,nc_type datatype,int _s)629   void minc_1_writer::open(const char *path,const minc_info& inf,int slice_dimensions,nc_type datatype,int _s)
630   {
631 #ifndef WIN32
632     ncopts = 0;
633 #endif
634 	  _info=inf;
635     //int  mdims[MAX_VAR_DIMS];
636     double vrange[2];
637     _write_prepared=false;
638 
639     _mincid = micreate((char*)path, NC_CLOBBER/*|MI2_CREATE_V2*/); //TODO: add environment variable checking
640 #ifdef MINC2
641     if (MI2_ISH5OBJ(_mincid)) { //micreate might create MINC2 file if environment variable is set
642       _minc2 = true;
643     }
644 #endif
645     _ndims=_info.size();
646     _datatype=datatype;
647     _slice_dimensions=slice_dimensions;
648     _is_signed=_s;
649     fill(_map_to_std.begin(),_map_to_std.end(),-1);
650     for(int i=_ndims-1;i>=0;i--)
651     {
652       //just a precaution
653       switch(_info[i].dim)
654       {
655         case dim_info::DIM_X:_info[i].name=MIxspace;_map_to_std[1]=i;break;
656         case dim_info::DIM_Y:_info[i].name=MIyspace;_map_to_std[2]=i;break;
657         case dim_info::DIM_Z:_info[i].name=MIzspace;_map_to_std[3]=i;break;
658         case dim_info::DIM_TIME:_info[i].name=MItime;_map_to_std[4]=i;break;
659         default:
660         case dim_info::DIM_VEC:_info[i].name=MIvector_dimension;_map_to_std[0]=i;break;
661         //default: REPORT_ERROR("Unknown Dimension!");
662       }
663       mdims[i]=ncdimdef(_mincid, _info[i].name.c_str(), _info[i].length);
664       if(_info[i].dim!=dim_info::DIM_VEC)
665       {
666         int dimid=micreate_std_variable(_mincid,(char*)_info[i].name.c_str(),NC_INT, 0, NULL);
667         miattputdbl(_mincid, dimid, (char*)MIstep,_info[i].step);
668         miattputdbl(_mincid, dimid, (char*)MIstart,_info[i].start);
669 
670         if(_info[i].have_dir_cos)
671           ncattput(_mincid, dimid, (char*)MIdirection_cosines,NC_DOUBLE, 3, _info[i].dir_cos);
672       }
673     }
674     _slab_len=1;
675     for(int i=0;i<_slice_dimensions;i++)
676     {
677       _slab[_ndims-i-1]=_info[_ndims-i-1].length;
678       _slab_len*=_info[_ndims-i-1].length;
679     }
680 
681     _icmax=_icmin=MI_ERROR;
682     //ncopts = NC_OPTS_VAL;
683     _imgid=micreate_std_variable(_mincid, (char*)MIimage, _datatype, _ndims, mdims);
684     _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX;
685 
686     switch(_datatype)
687     {
688       case NC_DOUBLE:
689         vrange[0]=-DBL_MAX;vrange[1]=DBL_MAX;
690         _is_signed=1;
691         break;
692       case NC_FLOAT:
693         vrange[0]=-FLT_MAX;vrange[1]=FLT_MAX;
694         _is_signed=1;
695         break;
696       case NC_SHORT:
697         if(_is_signed)
698         {
699           vrange[0]=SHRT_MIN;
700           vrange[1]=SHRT_MAX;
701         } else {
702           vrange[0]=0;vrange[1]=USHRT_MAX;
703         }
704         break;
705       case NC_BYTE:
706         if(_is_signed)
707         {
708           vrange[0]=-128;vrange[1]=127;
709         } else {
710           vrange[0]=0;vrange[1]=255;
711         }
712         break;
713       case NC_INT:
714         if(_is_signed)
715         {
716           vrange[0]=INT_MIN;vrange[1]=INT_MAX;
717         }else{
718           vrange[0]=0;vrange[1]=UINT_MAX;
719         }
720         break;
721       default:break;
722     };
723     miattputstr(_mincid, _imgid, (char*)MIcomplete, (char*)MI_FALSE);
724     miattputstr(_mincid, _imgid, (char*)MIsigntype, (char*)(_is_signed?MI_SIGNED:MI_UNSIGNED));
725     ncattput(_mincid, _imgid, (char*)MIvalid_range, NC_DOUBLE, 2, vrange);
726     miset_valid_range(_mincid, _imgid, vrange);
727   }
728 
open(const char * path,const minc_1_base & imitate)729   void minc_1_writer::open(const char *path,const minc_1_base& imitate)
730   {
731 
732     open(path,imitate.info(),imitate.slice_dimensions(),
733          imitate.datatype(),imitate.is_signed());
734 
735     copy_headers(imitate);
736   }
737 
open(const char * path,const char * imitate_file)738   void minc_1_writer::open(const char *path,const char *imitate_file)
739   {
740     minc_1_reader rdr;
741     //open minc file in metadate mode
742     rdr.open(imitate_file,false,true);
743     open(path,rdr);
744     //copy_headers(rdr);
745   }
746 
setup_write_float()747   void minc_1_writer::setup_write_float()
748   {
749     _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX;
750 
751     switch(_datatype)
752     {
753       case NC_DOUBLE:
754         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
755         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
756 
757         _set_image_range=true;
758         _set_slice_range=false;
759         break;
760 
761       case NC_FLOAT:
762         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
763         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
764 
765         _set_image_range=true;
766         _set_slice_range=false;
767         break;
768 
769       case NC_SHORT:
770 
771           _set_image_range=false;
772           _set_slice_range=true;
773 
774           _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
775           _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
776         break;
777 
778       case NC_BYTE:
779         _set_image_range=false;
780         _set_slice_range=true;
781 
782         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
783         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
784         break;
785       case NC_INT:
786         _set_image_range=false;
787         _set_slice_range=true;
788 
789         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
790         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
791         break;
792 
793       default:
794         break;
795     };
796     ncendef(_mincid);
797 
798     if(_datatype==NC_DOUBLE || _datatype==NC_FLOAT)
799     {
800       miicv_setstr(_icvid, MI_ICV_SIGN,   (char*)MI_SIGNED);
801       miicv_setint(_icvid, MI_ICV_TYPE,    NC_FLOAT);
802       miicv_setint(_icvid, MI_ICV_DO_NORM,    true);
803       miicv_setint(_icvid, MI_ICV_USER_NORM, true);
804       miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -FLT_MAX);
805       miicv_setdbl(_icvid, MI_ICV_VALID_MAX, FLT_MAX);
806       _calc_min_max=true;
807 
808     } else { //do something smart here?
809       miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
810       miicv_setint(_icvid, MI_ICV_TYPE, NC_FLOAT);
811       miicv_setint(_icvid, MI_ICV_DO_NORM, false);
812       //miicv_setint(_icvid, MI_ICV_USER_NORM, false);
813       miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -FLT_MAX);
814       miicv_setdbl(_icvid, MI_ICV_VALID_MAX, FLT_MAX);
815       _calc_min_max=true;
816 
817     }
818     miicv_attach(_icvid, _mincid, _imgid);
819     _io_datatype=NC_FLOAT;
820     _write_prepared=true;
821   }
822 
setup_write_double()823   void minc_1_writer::setup_write_double()
824   {
825     _image_range[0]=DBL_MAX;_image_range[1]=-DBL_MAX;
826 
827     switch(_datatype)
828     {
829       case NC_DOUBLE:
830         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
831         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
832 
833         _set_image_range=true;
834         _set_slice_range=false;
835         break;
836 
837       case NC_FLOAT:
838         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
839         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
840 
841         _set_image_range=true;
842         _set_slice_range=false;
843         break;
844 
845       case NC_SHORT:
846 
847           _set_image_range=false;
848           _set_slice_range=true;
849 
850           _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
851           _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
852         break;
853 
854       case NC_BYTE:
855         _set_image_range=false;
856         _set_slice_range=true;
857 
858         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
859         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
860         break;
861       case NC_INT:
862         _set_image_range=false;
863         _set_slice_range=true;
864 
865         _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
866         _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, _ndims-_slice_dimensions, mdims);
867         break;
868 
869       default:
870         break;
871     };
872     ncendef(_mincid);
873 
874     if(_datatype==NC_DOUBLE)
875     {
876       miicv_setstr(_icvid, MI_ICV_SIGN,   (char*)MI_SIGNED);
877       miicv_setint(_icvid, MI_ICV_TYPE,    NC_DOUBLE);
878       miicv_setint(_icvid, MI_ICV_DO_NORM,    true);
879       miicv_setint(_icvid, MI_ICV_USER_NORM, true);
880       miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX);
881       miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX);
882       _calc_min_max=true;
883 
884     } else if(_datatype==NC_FLOAT)  {
885       miicv_setstr(_icvid, MI_ICV_SIGN,   (char*)MI_SIGNED);
886       miicv_setint(_icvid, MI_ICV_TYPE,    NC_DOUBLE);
887       miicv_setint(_icvid, MI_ICV_DO_NORM,    true);
888       miicv_setint(_icvid, MI_ICV_USER_NORM, true);
889       miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX);
890       miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX);
891       _calc_min_max=true;
892 
893     } else { //do something smart here?
894       miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
895       miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE);
896       miicv_setint(_icvid, MI_ICV_DO_NORM, false);
897       //miicv_setint(_icvid, MI_ICV_USER_NORM, false);
898       miicv_setdbl(_icvid, MI_ICV_VALID_MIN, -DBL_MAX);
899       miicv_setdbl(_icvid, MI_ICV_VALID_MAX, DBL_MAX);
900       _calc_min_max=true;
901 
902     }
903     miicv_attach(_icvid, _mincid, _imgid);
904     _io_datatype=NC_DOUBLE;
905     _write_prepared=true;
906   }
907 
setup_write_short(bool n)908   void minc_1_writer::setup_write_short(bool n)
909   {
910     _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
911     _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
912     _set_image_range=true;
913     _set_slice_range=false;
914 
915     ncendef(_mincid);
916 
917     miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT);
918     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
919     //miicv_setstr(_icvid, MI_ICV_SIGN, true);
920     /* Set range of values */ //TODO: set this to something sensible?
921     miicv_setint(_icvid, MI_ICV_VALID_MIN, SHRT_MIN);
922     miicv_setint(_icvid, MI_ICV_VALID_MAX, SHRT_MAX);
923 
924     /* No normalization so that pixels are scaled to the slice */
925     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
926     miicv_setint(_icvid, MI_ICV_DO_RANGE, false);
927 
928     miicv_attach(_icvid, _mincid, _imgid);
929 
930     _io_datatype=NC_SHORT;
931     _write_prepared=true;
932   }
933 
setup_write_ushort(bool n)934   void minc_1_writer::setup_write_ushort(bool n)
935   {
936     _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
937     _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
938     _set_image_range=true;
939     _set_slice_range=false;
940 
941     ncendef(_mincid);
942     miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT);
943     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
944     miicv_setstr(_icvid, MI_ICV_SIGN, false);
945 
946     /* Set range of values */ //TODO: set this to something sensible?
947     miicv_setint(_icvid, MI_ICV_VALID_MIN, 0);
948     miicv_setint(_icvid, MI_ICV_VALID_MAX, USHRT_MAX);
949 
950     /* No normalization so that pixels are scaled to the slice */
951     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
952     miicv_setint(_icvid, MI_ICV_DO_RANGE, false);
953 
954     miicv_attach(_icvid, _mincid, _imgid);
955     _io_datatype=NC_SHORT;
956     _write_prepared=true;
957   }
958 
setup_write_byte(bool n)959   void minc_1_writer::setup_write_byte(bool n)
960   {
961     _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
962     _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
963     _set_image_range=true;
964     _set_slice_range=false;
965 
966     ncendef(_mincid);
967     miicv_setint(_icvid, MI_ICV_TYPE, NC_BYTE);
968     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
969 
970     /* Set range of values */ //TODO: set this to something sensible?
971     miicv_setint(_icvid, MI_ICV_VALID_MIN, 0);
972     miicv_setint(_icvid, MI_ICV_VALID_MAX, UCHAR_MAX);
973 
974     /* No normalization so that pixels are scaled to the slice */
975     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
976     miicv_setint(_icvid, MI_ICV_DO_RANGE, false);
977 
978     miicv_attach(_icvid, _mincid, _imgid);
979 
980     _io_datatype=NC_BYTE;
981     _write_prepared=true;
982   }
983 
setup_write_int(bool n)984   void minc_1_writer::setup_write_int(bool n)
985   {
986     _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
987     _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
988     _set_image_range=true;
989     _set_slice_range=false;
990 
991     ncendef(_mincid);
992     miicv_setint(_icvid, MI_ICV_TYPE, NC_INT);
993     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
994 
995     /* Set range of values */ //TODO: set this to something sensible?
996     miicv_setint(_icvid, MI_ICV_VALID_MIN, INT_MIN);
997     miicv_setint(_icvid, MI_ICV_VALID_MAX, INT_MAX);
998 
999     /* No normalization so that pixels are scaled to the slice */
1000     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
1001     miicv_setint(_icvid, MI_ICV_DO_RANGE, false);
1002 
1003     miicv_attach(_icvid, _mincid, _imgid);
1004 
1005 
1006     _io_datatype=NC_INT;
1007     _write_prepared=true;
1008   }
1009 
setup_write_uint(bool n)1010   void minc_1_writer::setup_write_uint(bool n)
1011   {
1012     _icmax=micreate_std_variable(_mincid, (char*)MIimagemax, NC_DOUBLE, 0, NULL);
1013     _icmin=micreate_std_variable(_mincid, (char*)MIimagemin, NC_DOUBLE, 0, NULL);
1014     _set_image_range=true;
1015     _set_slice_range=false;
1016 
1017     ncendef(_mincid);
1018     miicv_setint(_icvid, MI_ICV_TYPE, NC_INT);
1019     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
1020 
1021     /* Set range of values */ //TODO: set this to something sensible?
1022     miicv_setint(_icvid, MI_ICV_VALID_MIN, 0);
1023     miicv_setint(_icvid, MI_ICV_VALID_MAX, UINT_MAX);
1024 
1025     /* No normalization so that pixels are scaled to the slice */
1026     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
1027     miicv_setint(_icvid, MI_ICV_DO_RANGE, false);
1028 
1029     miicv_attach(_icvid, _mincid, _imgid);
1030     _io_datatype=NC_INT;
1031     _write_prepared=true;
1032   }
1033 
_setup_dimensions(void)1034   void minc_1_reader::_setup_dimensions(void)
1035   {
1036     if(_metadate_only)
1037       REPORT_ERROR("Minc file in metadate only mode!");
1038     if(_positive_directions)
1039     {
1040       /* We want to ensure that images have X, Y and Z dimensions in the
1041       positive direction, giving patient left on left and for drawing from
1042       bottom up. If we wanted patient right on left and drawing from
1043       top down, we would set to MI_ICV_NEGATIVE. */
1044 
1045       miicv_setint(_icvid, MI_ICV_DO_DIM_CONV, true);
1046       //TODO: make sure to change only x,y,z conversions here
1047       //miicv_setint(_icvid, MI_ICV_XDIM_DIR, 3);
1048       //we want to convert only X,Y,Z dimensions if they are present
1049       int num=(_map_to_std[1]>=0?1:0)+(_map_to_std[2]>=0?1:0)+(_map_to_std[3]>=0?1:0);
1050       miicv_setint(_icvid, MI_ICV_NUM_IMGDIMS, num);
1051 
1052       if(_map_to_std[1]>=0)
1053       {
1054         miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[1],-1);
1055         miicv_setint(_icvid, MI_ICV_XDIM_DIR, MI_ICV_POSITIVE);
1056       }
1057 
1058       if(_map_to_std[2]>=0)
1059       {
1060         miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[2],-1);
1061         miicv_setint(_icvid, MI_ICV_YDIM_DIR, MI_ICV_POSITIVE);
1062       }
1063 
1064       if(_map_to_std[3]>=0)
1065       {
1066         miicv_setint(_icvid, MI_ICV_DIM_SIZE+_map_to_std[3],-1);
1067         miicv_setint(_icvid, MI_ICV_ZDIM_DIR, MI_ICV_POSITIVE);
1068       }
1069     }
1070     miicv_setint(_icvid, MI_ICV_DO_SCALAR, false);
1071   }
1072 
setup_read_float(void)1073   void minc_1_reader::setup_read_float(void)
1074   {
1075     if(_metadate_only)
1076       REPORT_ERROR("Minc file in metadate only mode!");
1077     miicv_setint(_icvid, MI_ICV_TYPE, NC_FLOAT);
1078     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1079     miicv_setint(_icvid, MI_ICV_USER_NORM, true);
1080     /* Make sure that any out of range values are mapped to lowest value
1081       of type (for input only) */
1082     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1083 
1084     _setup_dimensions();
1085     miicv_attach(_icvid, _mincid, _imgid);
1086     _io_datatype=NC_FLOAT;
1087     _read_prepared=true;
1088   }
1089 
setup_read_double(void)1090   void minc_1_reader::setup_read_double(void)
1091   {
1092     if(_metadate_only)
1093       REPORT_ERROR("Minc file in metadate only mode!");
1094     miicv_setint(_icvid, MI_ICV_TYPE, NC_DOUBLE);
1095     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1096     miicv_setint(_icvid, MI_ICV_USER_NORM, true);
1097     /* Make sure that any out of range values are mapped to lowest value
1098       of type (for input only) */
1099     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1100 
1101     _setup_dimensions();
1102     miicv_attach(_icvid, _mincid, _imgid);
1103     _io_datatype=NC_DOUBLE;
1104     _read_prepared=true;
1105   }
1106 
setup_read_short(bool n)1107   void minc_1_reader::setup_read_short(bool n)
1108   {
1109     if(_metadate_only)
1110       REPORT_ERROR("Minc file in metadate only mode!");
1111 
1112     miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT);
1113     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
1114     /* Set range of values */
1115     miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]);
1116     miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]);
1117 
1118     /* Do normalization so that all pixels are on same scale */
1119     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1120     //miicv_setint(_icvid, MI_ICV_USER_NORM, true);
1121     /* Make sure that any out of range values are mapped to lowest value
1122       of type (for input only) */
1123     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1124 
1125     _setup_dimensions();
1126 
1127     miicv_attach(_icvid, _mincid, _imgid);
1128     _io_datatype=NC_SHORT;
1129     _read_prepared=true;
1130   }
1131 
setup_read_ushort(bool n)1132   void minc_1_reader::setup_read_ushort(bool n)
1133   {
1134     if(_metadate_only)
1135       REPORT_ERROR("Minc file in metadate only mode!");
1136 
1137     miicv_setint(_icvid, MI_ICV_TYPE, NC_SHORT);
1138     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
1139     /* Set range of values */
1140     miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]);
1141     miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]);
1142 
1143     /* Do normalization so that all pixels are on same scale */
1144     miicv_setint(_icvid, MI_ICV_DO_NORM, false);
1145     //miicv_setint(_icvid, MI_ICV_USER_NORM, true);
1146     /* Make sure that any out of range values are mapped to lowest value
1147       of type (for input only) */
1148     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1149 
1150     _setup_dimensions();
1151 
1152     miicv_attach(_icvid, _mincid, _imgid);
1153     _io_datatype=NC_SHORT;
1154     _read_prepared=true;
1155   }
1156 
setup_read_byte(bool n)1157   void minc_1_reader::setup_read_byte(bool n)
1158   {
1159     if(_metadate_only)
1160       REPORT_ERROR("Minc file in metadate only mode!");
1161     miicv_setint(_icvid, MI_ICV_TYPE, NC_BYTE);
1162     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
1163     /* Set range of values */
1164     miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]);
1165     miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]);
1166 
1167    /* Do normalization so that all pixels are on same scale */
1168     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1169     /* Make sure that any out of range values are mapped to lowest value
1170       of type (for input only) */
1171     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1172 
1173     _setup_dimensions();
1174     miicv_attach(_icvid, _mincid, _imgid);
1175     _io_datatype=NC_BYTE;
1176     _read_prepared=true;
1177   }
1178 
setup_read_int(bool n)1179   void minc_1_reader::setup_read_int(bool n)
1180   {
1181     if(_metadate_only)
1182       REPORT_ERROR("Minc file in metadate only mode!");
1183     miicv_setint(_icvid, MI_ICV_TYPE, NC_INT);
1184     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_SIGNED);
1185     /* Set range of values */
1186     miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]);
1187     miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]);
1188 
1189    /* Do normalization so that all pixels are on same scale */
1190     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1191     /* Make sure that any out of range values are mapped to lowest value
1192       of type (for input only) */
1193     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1194 
1195     _setup_dimensions();
1196     miicv_attach(_icvid, _mincid, _imgid);
1197     _io_datatype=NC_INT;
1198     _read_prepared=true;
1199   }
1200 
setup_read_uint(bool n)1201   void minc_1_reader::setup_read_uint(bool n)
1202   {
1203     if(_metadate_only)
1204       REPORT_ERROR("Minc file in metadate only mode!");
1205     miicv_setint(_icvid, MI_ICV_TYPE, NC_INT);
1206     miicv_setstr(_icvid, MI_ICV_SIGN, (char*)MI_UNSIGNED);
1207     /* Set range of values */
1208     miicv_setdbl(_icvid, MI_ICV_VALID_MIN, _image_range[0]);
1209     miicv_setdbl(_icvid, MI_ICV_VALID_MAX, _image_range[1]);
1210 
1211    /* Do normalization so that all pixels are on same scale */
1212     miicv_setint(_icvid, MI_ICV_DO_NORM, true);
1213     /* Make sure that any out of range values are mapped to lowest value
1214       of type (for input only) */
1215     miicv_setint(_icvid, MI_ICV_DO_FILLVALUE, true);
1216 
1217     _setup_dimensions();
1218     miicv_attach(_icvid, _mincid, _imgid);
1219     _io_datatype=NC_INT;
1220     _read_prepared=true;
1221   }
1222 
read(void * buffer)1223   void minc_1_reader::read(void* buffer)
1224   {
1225     if(!_read_prepared)
1226       REPORT_ERROR("Not ready to read, use setup_read_XXXX");
1227 
1228     miicv_get(_icvid, &_cur[0], &_slab[0], buffer);
1229   }
1230 
write(void * buffer)1231   void minc_1_writer::write(void* buffer)
1232   {
1233     if(!_write_prepared)
1234       REPORT_ERROR("Not ready to write, use setup_write_XXXX");
1235 
1236     double r_min= DBL_MAX; //slab minimal value
1237     double r_max=-DBL_MAX; //slab maximal value
1238     //int irmin=0,irmax=0;
1239     if(_calc_min_max )
1240     {
1241       if(_io_datatype==NC_FLOAT)
1242       {//
1243         float *tmp=(float*)buffer;
1244         for(int i=0;i<_slab_len;i++)
1245         {
1246           if(r_min>tmp[i]) r_min=tmp[i];//irmin=i;
1247           if(r_max<tmp[i]) r_max=tmp[i];//irmax=i;
1248         }
1249       } else if(_io_datatype==NC_DOUBLE) {
1250           double *tmp=(double*)buffer;
1251           for(int i=0;i<_slab_len;i++)
1252           {
1253             if(r_min>tmp[i]) r_min=tmp[i];//irmin=i;
1254             if(r_max<tmp[i]) r_max=tmp[i];//irmax=i;
1255           }
1256       } else if(_io_datatype==NC_SHORT && !_is_signed) {
1257         unsigned short *tmp=(unsigned short *)buffer;
1258         for(int i=0;i<_slab_len;i++)
1259         {
1260           if(r_min>tmp[i]) r_min=tmp[i];
1261           if(r_max<tmp[i]) r_max=tmp[i];
1262         }
1263       } else if(_io_datatype==NC_SHORT && _is_signed) {
1264         short *tmp=(short *)buffer;
1265         for(int i=0;i<_slab_len;i++)
1266         {
1267           if(r_min>tmp[i]) r_min=tmp[i];
1268           if(r_max<tmp[i]) r_max=tmp[i];
1269         }
1270       } else if(_io_datatype==NC_BYTE) {
1271         unsigned char *tmp=(unsigned char *)buffer;
1272         for(int i=0;i<_slab_len;i++)
1273         {
1274           if(r_min>tmp[i]) r_min=tmp[i];
1275           if(r_max<tmp[i]) r_max=tmp[i];
1276         }
1277       } else if(_io_datatype==NC_INT && _is_signed) {
1278         int *tmp=(int *)buffer;
1279         for(int i=0;i<_slab_len;i++)
1280         {
1281           if(r_min>tmp[i]) r_min=tmp[i];
1282           if(r_max<tmp[i]) r_max=tmp[i];
1283         }
1284       } else if(_io_datatype==NC_INT && !_is_signed) {
1285         unsigned int *tmp=(unsigned int *)buffer;
1286         for(int i=0;i<_slab_len;i++)
1287         {
1288           if(r_min>tmp[i]) r_min=tmp[i];
1289           if(r_max<tmp[i]) r_max=tmp[i];
1290         }
1291       }
1292       /*
1293       if(r_min<-10000||r_max>10000)
1294       {
1295         std::cerr<<r_min<<":"<<r_max<<" ";
1296         for(int i=0;i<4;i++)
1297           std::cerr<<_cur[i]<<",";
1298         std::cerr<<"  "<<r_min<<":"<<r_max;
1299         std::cerr<<" "<<irmin<<":"<<irmax<<"   "<<_slab_len;
1300         std::cerr<<std::endl;
1301       }*/
1302 
1303       if(_set_slice_range)
1304       {
1305         miicv_detach(_icvid);
1306         miicv_setdbl(_icvid, MI_ICV_VALID_MIN, r_min);
1307         miicv_setdbl(_icvid, MI_ICV_VALID_MAX, r_max);
1308         miicv_attach(_icvid, _mincid, _imgid);
1309       }
1310 
1311       if(_set_slice_range)
1312       {
1313         mivarput1(_mincid, _icmin, &_cur[0], NC_DOUBLE, NULL, &r_min);
1314         mivarput1(_mincid, _icmax, &_cur[0], NC_DOUBLE, NULL, &r_max);
1315       }
1316 
1317       if(_image_range[0]>r_min) _image_range[0]=r_min;
1318       if(_image_range[1]<r_max) _image_range[1]=r_max;
1319     }
1320     miicv_put(_icvid, &_cur[0], &_slab[0], buffer);
1321   }
1322 
close(void)1323   void minc_1_writer::close(void)
1324   {
1325     if(_set_image_range)
1326     {
1327       mivarput1(_mincid, _icmin, 0, NC_DOUBLE, NULL, &_image_range[0]);
1328       mivarput1(_mincid, _icmax, 0, NC_DOUBLE, NULL, &_image_range[1]);
1329       miset_valid_range(_mincid, _imgid, _image_range);
1330       _set_image_range=false;
1331     }
1332     minc_1_base::close();
1333   }
1334 
1335 
~minc_1_writer()1336   minc_1_writer::~minc_1_writer()
1337   {
1338     close();
1339   }
1340 
copy_headers(const minc_1_base & src)1341   void minc_1_writer::copy_headers(const minc_1_base& src)
1342   {
1343 
1344     //code copied from mincresample
1345     int nexcluded, excluded_vars[10];
1346     int varid;
1347 
1348     /* Create the list of excluded variables */
1349     nexcluded = 0;
1350     //ncopts = 0;
1351 
1352     if ((varid=ncvarid(src.mincid(), MIxspace)) != MI_ERROR)
1353       excluded_vars[nexcluded++] = varid;
1354     if ((varid=ncvarid(src.mincid(), MIyspace)) != MI_ERROR)
1355       excluded_vars[nexcluded++] = varid;
1356     if ((varid=ncvarid(src.mincid(), MIzspace)) != MI_ERROR)
1357       excluded_vars[nexcluded++] = varid;
1358     if ((varid=ncvarid(src.mincid(), MItime)) != MI_ERROR)
1359       excluded_vars[nexcluded++] = varid;
1360     if ((varid=ncvarid(src.mincid(), MIimage)) != MI_ERROR)
1361       excluded_vars[nexcluded++] = varid;
1362     if ((varid=ncvarid(src.mincid(), MIimagemax)) != MI_ERROR)
1363       excluded_vars[nexcluded++] = varid;
1364     if ((varid=ncvarid(src.mincid(), MIimagemin)) != MI_ERROR)
1365       excluded_vars[nexcluded++] = varid;
1366     if ((varid=ncvarid(src.mincid(), "rootvariable")) != MI_ERROR)
1367       excluded_vars[nexcluded++] = varid;
1368 //#if MINC2
1369     if ((varid=ncvarid(src.mincid(), MIvector_dimension)) != MI_ERROR)
1370       excluded_vars[nexcluded++] = varid;
1371 //#endif /* MINC2 */
1372     //ncopts = NC_VERBOSE | NC_FATAL;
1373     /* Copy all other variable definitions */
1374     micopy_all_var_defs(src.mincid(), _mincid, nexcluded, excluded_vars);
1375   }
1376 
1377   //! append a line into minc history
append_history(const char * append_history)1378   void minc_1_writer::append_history(const char *append_history)
1379   {
1380     nc_type datatype;
1381     int att_length;
1382     //ncopts=0;
1383     if ((ncattinq(_mincid, NC_GLOBAL, MIhistory, &datatype,&att_length) == MI_ERROR) ||
1384         (datatype != NC_CHAR))
1385       att_length = 0;
1386     att_length += strlen(append_history) + 1;
1387     char* str = new char[att_length];
1388     str[0] = '\0';
1389     miattgetstr(_mincid, NC_GLOBAL, (char*)MIhistory, att_length,str);
1390     //ncopts=NC_VERBOSE | NC_FATAL;
1391     strcat(str, append_history);
1392     miattputstr(_mincid, NC_GLOBAL, (char*)MIhistory, str);
1393     delete [] str;
1394   }
1395 
create_var_id(const char * varname)1396   int minc_1_base::create_var_id(const char *varname)
1397   {
1398     int old_ncopts = ncopts; ncopts = 0;
1399     int res=var_id(varname);
1400     if(res==MI_ERROR) //need to create a variable
1401       res=micreate_group_variable(_mincid,(char*)varname);//ncvardef(_mincid,varname,NC_INT,0,0);
1402     if(res==MI_ERROR) //need to create a variable
1403       res=ncvardef(_mincid,varname,NC_INT,0,0);
1404     ncopts = old_ncopts;
1405     return res;
1406   }
1407 
insert(const char * varname,const char * attname,double val)1408   void minc_1_base::insert(const char *varname,const char *attname,double val)
1409   {
1410     ncattput(_mincid, create_var_id(varname),attname, NC_DOUBLE, 1, (void *) &val);
1411   }
1412 
insert(const char * varname,const char * attname,const char * val)1413   void minc_1_base::insert(const char *varname,const char *attname,const char* val)
1414   {
1415     ncattput(_mincid, create_var_id(varname),attname, NC_CHAR, strlen(val) + 1, (void *) val);
1416   }
1417 
insert(const char * varname,const char * attname,const std::vector<double> & val)1418   void minc_1_base::insert(const char *varname,const char *attname,const std::vector<double> &val)
1419   {
1420     ncattput(_mincid, create_var_id(varname),attname, NC_DOUBLE, val.size(), (void *) &val[0]);
1421   }
1422 
insert(const char * varname,const char * attname,const std::vector<int> & val)1423   void minc_1_base::insert(const char *varname,const char *attname,const std::vector<int> &val)
1424   {
1425     ncattput(_mincid, create_var_id(varname),attname, NC_INT, val.size(), (void *) &val[0]);
1426   }
1427 
insert(const char * varname,const char * attname,const std::vector<short> & val)1428   void minc_1_base::insert(const char *varname,const char *attname,const std::vector<short> &val)
1429   {
1430     ncattput(_mincid, create_var_id(varname),attname, NC_SHORT, val.size(), (void *) &val[0]);
1431   }
1432 
insert(const char * varname,const char * attname,const std::vector<unsigned char> & val)1433   void minc_1_base::insert(const char *varname,const char *attname,const std::vector<unsigned char> &val)
1434   {
1435     ncattput(_mincid, create_var_id(varname),attname, NC_BYTE, val.size(), (void *) &val[0]);
1436   }
1437 
1438 };
1439