#include "fileio.h"
#include "utils.h"
#include "complexdata.h" // for unit test

#include <tjutils/tjtypes.h>
#include <tjutils/tjindex.h>
#include <tjutils/tjtest.h>

#include <tjutils/tjlog_code.h>

const char* FileIO::get_compName() {return "FileIO";}
LOGGROUNDWORK(FileIO)


///////////////////////////////////////////////////////////////////////////////


void FileFormat::register_format() {
  svector suffs=suffix();
  for(unsigned int i=0; i<suffs.size(); i++) formats[suffs[i]].push_back(this);
}


int FileFormat::read(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol& prot) {
  Log<FileIO> odinlog("FileFormat","read");
  ODINLOG(odinlog,errorLog) << this->description() << "::read not implemented" << STD_endl;
  return -1;
}

int FileFormat::read(FileIO::ProtocolDataMap& pdmap, const STD_string& filename, const FileReadOpts& opts, const Protocol& protocol_template) {
  Data<float,4> data;
  Protocol prot(protocol_template);
  int result=this->read(data, filename, opts, prot);
  if(result<0) return -1;
  if(result>0) pdmap[prot].reference(data); // Ignore zero-size data, e.g. DICOMDIR
  return result;
}


int FileFormat::write(const Data<float,4>& data, const STD_string& filename, const FileWriteOpts& opts, const Protocol& prot) {
  Log<FileIO> odinlog("FileFormat","write");
  ODINLOG(odinlog,errorLog) << this->description() << "::write not implemented" << STD_endl;
  return -1;
}


int FileFormat::write(const FileIO::ProtocolDataMap& pdmap, const STD_string& filename, const FileWriteOpts& opts) {
  Log<FileIO> odinlog("FileFormat","write");
  int result=0;

  svector fnames=create_unique_filenames(filename, pdmap, opts.fnamepar);
  ODINLOG(odinlog,normalDebug) << "fnames=" << fnames.printbody() << STD_endl;

  int ifile=0;
  for(FileIO::ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {
    int writeresult=this->write(pdit->second, fnames[ifile], opts, pdit->first);
    if(writeresult<0) return writeresult;
    result+=writeresult;
    ifile++;
  }
  return result;
}




STD_string FileFormat::analyze_suffix(const STD_string& filename) {
  LDRfileName fname(filename);
  return fname.get_suffix();
}

FileFormat* FileFormat::get_format(const STD_string& filename, const STD_string& override_suffix) {
  Log<FileIO> odinlog("FileFormat","get_format");

  STD_string suffix;
  if(override_suffix==AUTODETECTSTR) {
    suffix=analyze_suffix(filename);
  } else {
    suffix=override_suffix;
  }
  ODINLOG(odinlog,normalDebug) << "filename:override_suffix:suffix=" << filename << ":" << override_suffix << ":" << suffix << STD_endl;

  if(formats.find(suffix)!=formats.end()) {
    const FormatList& formatlist=formats[suffix];
    ODINLOG(odinlog,normalDebug) << "#formats(" << suffix << ")=" << formatlist.size() << STD_endl;
    if(formatlist.size()>1) {
       ODINLOG(odinlog,errorLog) << "Ambiguous file extension >" << analyze_suffix(filename) << "<" << STD_endl;
       ODINLOG(odinlog,errorLog) << "Use -wf/-rf option with unique identifier (e.g. -wf analyze)" << STD_endl;
       return 0;
    }
    return *(formatlist.begin());
  }
  return 0;
}


svector FileFormat::possible_formats() {
  svector result; result.resize(formats.size());
  int index=0;
  for(FormatMap::const_iterator it=formats.begin(); it!=formats.end(); ++it) {
    result[index]=it->first;
    index++;
  }
  return result;
}

STD_string FileFormat::formats_str(const STD_string& indent) {
  STD_string result;
  for(FormatMap::const_iterator it=formats.begin(); it!=formats.end(); ++it) {
    const FormatList& formatlist=it->second;
    for(FormatList::const_iterator it2=formatlist.begin(); it2!=formatlist.end(); ++it2) {
      result+=indent+(it->first)+" \t ("+(*it2)->description();
      svector dialects=(*it2)->dialects();
      if(dialects.size()) result+=", dialects: "+dialects.printbody();
      result+=")\n";
    }
  }
  return result;
}


void FileFormat::format_error(const STD_string& filename) {
  Log<FileIO> odinlog("FileFormat","format_error");
  ODINLOG(odinlog,errorLog) << "File extension >" << analyze_suffix(filename) << "< of file >" << filename << "< not recognized" << STD_endl;
  ODINLOG(odinlog,errorLog) << "Recognized file extensions (and formats) are" << STD_endl << formats_str("") << STD_endl;
}


svector FileFormat::create_unique_filenames(const STD_string& filename, const FileIO::ProtocolDataMap& pdmap, const STD_string& par) {
  Log<FileIO> odinlog("FileFormat","create_unique_filenames");

  unsigned int nnames=pdmap.size();
  svector result; result.resize(nnames);

  ODINLOG(odinlog,normalDebug) << "nnames=" << nnames << STD_endl;

  if(nnames==1) {
    result[0]=filename;
    return result;
  }

  STD_string serDesc;
  int serNum;

  // Get max series number
  int maxseries=0;
  for(FileIO::ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {
    pdit->first.study.get_Series(serDesc,serNum);
    maxseries=STD_max(maxseries,serNum);
  }

  svector pars;
  pars=tokens(par);
  int npars=pars.size();

  // Create vector of strings "S<serNumber>[_<serDescription>]"
  STD_map<STD_string,unsigned int> namescount;
  svector unique_names; unique_names.resize(nnames);
  unsigned int iname=0;
  for(FileIO::ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {
    pdit->first.study.get_Series(serDesc,serNum);
    STD_string str("S"+itos(serNum,maxseries));
    if(serDesc!="") str+=STD_string("_")+serDesc;
    
    if(npars) {
      for(int ipar=0; ipar<npars; ipar++) {
        str+=STD_string("_")+rmblock(pdit->first.printval(pars[ipar], true), "\n", ""); // truncate after first newline
      }
    }

    //replace characters not suitable for file names
    for(unsigned int i=0;i<str.length();i++) {
      const char cr=str[i];
      if(cr==' ' || cr==':' || cr==';' || cr=='/' || cr=='&' || cr=='*') str[i]='_';
    }

    if(namescount.find(str)==namescount.end()) namescount[str]=1;
    else namescount[str]++;

    unique_names[iname]=str;
    iname++;
  }
  ODINLOG(odinlog,normalDebug) << "unique_names=" << unique_names.printbody() << STD_endl;


  for(STD_map<STD_string,unsigned int>::const_iterator it=namescount.begin(); it!=namescount.end(); ++it) {
    ODINLOG(odinlog,normalDebug) << "namescount(" << it->first << ")=" << it->second << STD_endl;
  }
 

  // Make sure filenames are unique by appending enumeration, if neccessary
  STD_map<STD_string,unsigned int> uniqueindex; // separate count for each non-unique filename
  for(unsigned int i=0; i<nnames; i++) {
    STD_string str=unique_names[i];
    int ncount=namescount[str];
    if(ncount>1) {
      if(uniqueindex.find(str)==uniqueindex.end()) uniqueindex[str]=0;
      else uniqueindex[str]++;
      unique_names[i]+="_"+itos(uniqueindex[str],ncount);
    }
  }


  // Decompose file name to insert unique name before first known prefix
  LDRfileName fname(filename);
  svector fmts=FileIO::autoformats();
  LDRfileName bname=fname.get_basename_nosuffix();
  STD_string fnamesuff(fname.get_suffix());
  STD_string suffixstr;
  if(fnamesuff!="") suffixstr="."+fname.get_suffix();
  ODINLOG(odinlog,normalDebug) << "suffixstr=" << suffixstr << STD_endl;
  bool known_fmt=true;
  while(known_fmt) { // split filename from right to left as long as known file extensions are detected
    STD_string suff=bname.get_suffix();
    known_fmt=false;
    for(unsigned int i=0; i<fmts.size(); i++) {
      if(suff==fmts[i]) {
        known_fmt=true;
        break;
      }
    }
    if(known_fmt) {
      suffixstr="."+suff+suffixstr;
      bname=bname.get_basename_nosuffix();
    }
  }
  ODINLOG(odinlog,normalDebug) << "bname/suffixstr=" << bname << "/" << suffixstr << STD_endl;

  // Compose result
  for(unsigned int i=0; i<nnames; i++) {
    result[i]=fname.get_dirname()+SEPARATOR_STR;
    if(bname!="") result[i]+=bname+"_";
    result[i]+=unique_names[i]+suffixstr;
  }

  return result;
}



float FileFormat::voxel_extent(const Geometry& geometry, direction direction, int size) {
  Log<FileIO> odinlog("FileFormat","voxel_extent");
  float result;
  if(direction==sliceDirection) {
    if(geometry.get_Mode()==voxel_3d) {
      result=secureDivision(geometry.get_FOV(direction),  size);
    } else {
      if(geometry.get_nSlices()>1) result=geometry.get_sliceDistance();
      else                         result=geometry.get_sliceThickness();
    }
  } else {
    result=secureDivision(geometry.get_FOV(direction),  size);
  }

  ODINLOG(odinlog,normalDebug) << "result(" << directionLabel[direction] << ")=" << result << STD_endl;

  return result;
}


FileFormat::FormatMap FileFormat::formats;

//////////////////////////////////////////////////////////////////////////

void register_asc_format();
void register_dicom_format();
void register_gzip_format();
void register_hfss_format();
void register_interfile_format();
void register_ismrmrd_format();
void register_ser_format();
void register_mhd_format();
void register_mat_format();
void register_nifti_format();
void register_png_format();
void register_raw_format();
void register_Iris3D_format();
void register_vtk_format();

//////////////////////////////////////////////////////////////

struct FileFormatCreator : public StaticHandler<FileFormatCreator> {

  // functions to initialize/delete static members by the StaticHandler template class
  static void init_static() {
    register_asc_format();
    register_dicom_format();
    register_gzip_format();
    register_interfile_format();
    register_ismrmrd_format();
    register_ser_format();
    register_mhd_format();
    register_mat_format();
    register_nifti_format();
    register_png_format();
    register_Iris3D_format();
    register_raw_format();
    register_hfss_format();
    register_vtk_format();
  }
  static void destroy_static() {}

};

EMPTY_TEMPL_LIST bool StaticHandler<FileFormatCreator>::staticdone=false;

//////////////////////////////////////////////////////////////////////////


int FileIO::autoread(ProtocolDataMap& pdmap, const STD_string& filename, const FileReadOpts& opts, const Protocol& protocol_template,ProgressMeter* progmeter) {
  Log<FileIO> odinlog("FileIO","autoread");

  if(!checkdir(filename.c_str()) && filesize(filename.c_str())<=0) {
    ODINLOG(odinlog,errorLog) << "File " << filename << " not found or is empty" << STD_endl;
    return -1;
  }

  FileFormatCreator ffc; // for some reason, we needed a named object here instead of a plain constructor


  // override file suffix if jdx option is given
  STD_string suffix=opts.format;
  if(opts.ldr!="") suffix="smp";

  int result=-1;

  if(checkdir(filename.c_str())) {
    result=read_dir(pdmap, filename, opts, protocol_template,progmeter);
    if(result<0) return -1;

  } else {

    FileFormat* ff=FileFormat::get_format(filename,suffix);
    if(!ff) {
      FileFormat::format_error(filename);
      return -1;
    }

    ODINLOG(odinlog,loglevel()) << "Reading format " << ff->description() << STD_endl;
    result=ff->read(pdmap, filename, opts, protocol_template);
    if(result<0) {
      ODINLOG(odinlog,errorLog) << "Cannot read file " << filename << STD_endl;
      return -1;
    }
  }

  ODINLOG(odinlog,normalDebug) << "slicedist=" << pdmap.begin()->first.geometry.get_sliceDistance() << STD_endl;

  ProtocolDataMap pdmapin(pdmap);
  pdmap.clear();unsigned short cnt=0;
  for(ProtocolDataMap::const_iterator pdit=pdmapin.begin(); pdit!=pdmapin.end(); ++pdit) {

    if(!pdit->second.size()) {
      ODINLOG(odinlog,errorLog) << "Zero-size data while reading " << filename << STD_endl;
      return -1;
    }

    TinyVector<int,4> shape=pdit->second.shape();
    Protocol prot(pdit->first);

    // overwrite extent in protocol
    prot.seqpars.set_NumOfRepetitions(shape(timeDim));
    if(prot.geometry.get_Mode()==slicepack) prot.geometry.set_nSlices(shape(sliceDim));
    else prot.seqpars.set_MatrixSize(sliceDirection,shape(sliceDim));
    prot.seqpars.set_MatrixSize(phaseDirection,shape(phaseDim));
    prot.seqpars.set_MatrixSize(readDirection,shape(readDim));

    // Adjust total duration of sequence, only if uninitialized
    ODINLOG(odinlog,normalDebug) << "ExpDuration(TR)=" << prot.seqpars.get_ExpDuration() << STD_endl;
    prot.seqpars.set_ExpDuration(STD_max(prot.seqpars.get_ExpDuration(), shape(0)*prot.seqpars.get_RepetitionTime()/(60.0*1000.0)));

    // Set label of protocol to filename
    STD_string postfix;
    if(pdmapin.size()>1) postfix=":"+itos(cnt,pdmapin.size()-1);
    prot.set_label("Protocol of "+LDRfileName(filename).get_basename()+postfix);


    const Data<float,4>& data=pdit->second;
    if(!opts.fmap && data.is_filemapped()) {
      ODINLOG(odinlog,normalDebug) << "filemapped" << STD_endl;
      pdmap[prot].reference(data.copy());
    } else {
      pdmap[prot].reference(data);
    }

    STD_string cntstr;
    if(pdmapin.size()>1) cntstr="#"+itos(cnt)+" ";

    ODINLOG(odinlog,loglevel()) << "Read dataset " << cntstr << "from file " << filename << "  with dimensions " << shape << STD_endl;

    cnt++;
  }
  ODINLOG(odinlog,normalDebug) << "mem(read): " << Profiler::get_memory_usage() << STD_endl;


  // Extract only a certain dataset if requested
  int dset2extract=-1;
  if(opts.dset!="") {
    dset2extract=atoi(opts.dset.c_str());
    if(dset2extract>=int(pdmap.size()) || dset2extract<0) {
      ODINLOG(odinlog,errorLog) << "Dataset to extract (" << dset2extract << ") out of range (" << pdmap.size() << ")" << STD_endl;
      return -1;
    }
  }

  if(dset2extract>=0) {
    ODINLOG(odinlog,loglevel()) << "Extracting dataset with index " << dset2extract << STD_endl;

    ProtocolDataMap pdmap_copy(pdmap);
    pdmap.clear();

    int iset=0;
    for(ProtocolDataMap::const_iterator pdit=pdmap_copy.begin(); pdit!=pdmap_copy.end(); ++pdit) {
      if(dset2extract==iset) pdmap[pdit->first].reference(pdit->second);
      iset++;
    }
  }

  // filter on certain parameters in protocol
  if(opts.filter!="") {
    svector toks=tokens(opts.filter,'=');
    if(toks.size()==2) {
      ODINLOG(odinlog,loglevel()) << "Filtering data where parameter " << toks[0] << " contains " << toks[1] << STD_endl;

      ProtocolDataMap pdmap_copy(pdmap);
      pdmap.clear();

      for(ProtocolDataMap::const_iterator pdit=pdmap_copy.begin(); pdit!=pdmap_copy.end(); ++pdit) {
        if(pdit->first.printval(toks[0]).find(toks[1])!=STD_string::npos) pdmap[pdit->first].reference(pdit->second);
      }
    }
  }

  ODINLOG(odinlog,normalDebug) << "mem(result): " << Profiler::get_memory_usage() << STD_endl;


  return result;

}


//////////////////////////////////////////////////////////////////////////

struct ImageKey : public UniqueIndex<ImageKey> {
  double acquisition_time;
  double slice_loc;
  STD_string filename;
  STD_string datatype;

  ImageKey(double tt, double sl, const STD_string& fn, const STD_string& dt) : acquisition_time(tt), slice_loc(sl), filename(fn), datatype(dt)  {}

  bool operator < (const ImageKey& ik) const {
    if(slice_loc!=ik.slice_loc) return (slice_loc<ik.slice_loc); // slice location has highest priority
    if(acquisition_time!=ik.acquisition_time) return (acquisition_time<ik.acquisition_time); // sort by acquisition time stamp
    if(filename!=ik.filename) return (filename<ik.filename); // sort alphabetically
    return (get_index()<ik.get_index()); // last resort
  }

  // functions for UniqueIndex
  static const char* get_typename() {return "ImageKey";}
  static unsigned int get_max_instances() {return 0;}
};

//////////////////////////////////////////////////

typedef STD_map<ImageKey,Data<float,2> > KeyImageMap;

int FileIO::read_dir(ProtocolDataMap& pdmap, const STD_string& dirname, const FileReadOpts& opts, const Protocol& protocol_template,ProgressMeter* progmeter) {
  Log<FileIO> odinlog("FileIO","read_dir");

  ODINLOG(odinlog,normalDebug) << "mem(start): " << Profiler::get_memory_usage() << STD_endl;

  Range all=Range::all();

  int result=0;

  svector dirfiles=browse_dir(dirname,false,true);
  unsigned int nfiles=dirfiles.size();
  if(!nfiles) {
    ODINLOG(odinlog,errorLog) << "No files found in directory " << dirname << STD_endl;
    return -1;
  }

  STD_map<Protocol, KeyImageMap> oneimagemap;

  // Use local options to skip some of the global autoread options
  FileReadOpts opts_copy=opts;
  opts_copy.dset="";

  if(progmeter) progmeter->new_task( nfiles,"FileIO::read_dir");

  for(unsigned int i=0; i<nfiles; i++) {

    ProtocolDataMap onepdmap;

    STD_string dirfilename=dirname+SEPARATOR_STR+dirfiles[i];

    if(checkdir(dirfilename.c_str())){
      ODINLOG(odinlog,infoLog) << "Ignoring subdirectory " << dirfiles[i] << " while reading dir " << dirname << STD_endl;
      continue;
    }

    bool do_trace_cache=get_trace_status(); set_trace_status(false); // Do not trace reading of single files
    int readresult=FileIO::autoread(onepdmap,dirfilename,opts_copy,protocol_template);
    set_trace_status(do_trace_cache);
    if(progmeter && progmeter->increase_counter(dirfilename.c_str())) break;//cancel if the ProgressMeter says so

    if(readresult<0) {
      ODINLOG(odinlog,warningLog) << "Ignoring unrecognized file  " << dirfiles[i] << " while reading dir"  << STD_endl;
    } else if(readresult>0) { // ignore files unrecognized by the plugin (e.g. DICOMDIR)

      // Iterate over protocol-data pairs
      for(ProtocolDataMap::const_iterator pdit=onepdmap.begin(); pdit!=onepdmap.end(); ++pdit) {
        Protocol onefileprotcopy(pdit->first);
        const Data<float,4>& onefiledata=pdit->second;

        onefileprotcopy.use_acqstart_during_comparison(opts.framesplit);

        KeyImageMap& imgmap=oneimagemap[onefileprotcopy]; // Get/insert key-image map for this protocol


        // Check x- and y-size
        KeyImageMap::const_iterator it=imgmap.begin();
        if(it!=imgmap.end()) { // only check size if there already is an entry in the key-image map
          if(it->second.extent(0)!=onefiledata.extent(phaseDim)) {
            ODINLOG(odinlog,errorLog) << "Files in dir have different ysize"  << STD_endl;
            return -1;
          }
          if(it->second.extent(1)!=onefiledata.extent(readDim)) {
            ODINLOG(odinlog,errorLog) << "Files in dir have different xsize"  << STD_endl;
            return -1;
          }
        }


        STD_string dtype_prot=onefileprotcopy.system.get_data_type();

        // Insert repetitions and slices separately into map
        for(int irep=0; irep<onefiledata.extent(timeDim); irep++) {
          for(int islice=0; islice<onefiledata.extent(sliceDim); islice++) {

            double acqstart=onefileprotcopy.seqpars.get_AcquisitionStart()+double(irep)*onefileprotcopy.seqpars.get_RepetitionTime();
            double sliceloc=onefileprotcopy.geometry.get_sliceOffsetVector()[islice];
            ODINLOG(odinlog,normalDebug) << "acqstart/sliceloc=" << acqstart << "/" << sliceloc << STD_endl;
            Data<float,2>& dataref=imgmap[ImageKey(acqstart,sliceloc,dirfilename,dtype_prot)];
            dataref.resize(onefiledata.extent(phaseDim),onefiledata.extent(readDim));
            dataref=onefiledata(irep,islice,all,all);

            if(!progmeter) {
              ODINLOG(odinlog,loglevel()) << "Read slice from file " << dirfiles[i] << " "  <<  onefiledata.shape() << " with acquisition_start/slice_location=" << acqstart << "/" << sliceloc << STD_endl;
            }
          }
        }
      } // end iterating over protocol-data pairs
    }
  } // end iterating over files

  ODINLOG(odinlog,loglevel()) << "Found " << oneimagemap.size() << " protocol(s)" << STD_endl;

  ODINLOG(odinlog,normalDebug) << "mem(read): " << Profiler::get_memory_usage() << STD_endl;

  
  
  // Iterate over separate protocols to combine data
  int protindex=0;
  for(STD_map<Protocol, KeyImageMap>::iterator it=oneimagemap.begin(); it!=oneimagemap.end(); ++it) {
    KeyImageMap::const_iterator imgit;

    Protocol prot(it->first);
    KeyImageMap& imgmap=it->second; // non-const to be cleared after retrieving data from it

    TinyVector<int,4> shape;

    imgit=imgmap.begin();
    if(imgit==imgmap.end()) {
      ODINLOG(odinlog,errorLog) << "Empty imgmap"  << STD_endl;
      return -1;
    }

    shape(readDim)=imgit->second.extent(1);
    shape(phaseDim)=imgit->second.extent(0);

    // Get slice offsets and min acq time by iterating over images
    STD_list<float> slice_offsets;
    STD_list<STD_string> datatypes;
    double min_acqstart=imgit->first.acquisition_time;
    double max_acqstart=imgit->first.acquisition_time;
    for(imgit=imgmap.begin(); imgit!=imgmap.end(); ++imgit) {
      slice_offsets.push_back(imgit->first.slice_loc);
      datatypes.push_back(imgit->first.datatype);
      min_acqstart=STD_min(min_acqstart,imgit->first.acquisition_time);
      max_acqstart=STD_max(max_acqstart,imgit->first.acquisition_time);
    }


    // Find datatype with highest resolution in final dataset and use it for all datasets
    datatypes.sort();
    datatypes.unique();
    STD_string dtype=TypeTraits::type2label((float)0); // default  
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((u8bit) 0))!=datatypes.end()) dtype=TypeTraits::type2label((u8bit)0);
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((s8bit) 0))!=datatypes.end()) dtype=TypeTraits::type2label((s8bit)0);
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((u16bit)0))!=datatypes.end()) dtype=TypeTraits::type2label((u16bit)0);
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((s16bit)0))!=datatypes.end()) dtype=TypeTraits::type2label((s16bit)0);
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((u32bit)0))!=datatypes.end()) dtype=TypeTraits::type2label((u32bit)0);
    if(find(datatypes.begin(),datatypes.end(),TypeTraits::type2label((s32bit)0))!=datatypes.end()) dtype=TypeTraits::type2label((s32bit)0);
    if(datatypes.size()>1) {
      ODINLOG(odinlog,loglevel()) << "Using single datatype " << dtype << " for images with datatypes=" << svector(list2vector(datatypes)).printbody() << STD_endl;
    }
    prot.system.set_data_type(dtype);


    slice_offsets.sort();
    slice_offsets.unique();
    shape(sliceDim)=slice_offsets.size();
    int ndatasets=imgmap.size();
    shape(timeDim)=ndatasets/shape(sliceDim); // interpret remainder as time dimension
    int remainder=ndatasets%shape(sliceDim);
    if(remainder) {
      int missing=shape(sliceDim)-remainder;
      shape(timeDim)++;
      ODINLOG(odinlog,warningLog) << "Sparse data with " << missing << " slices missing, ndatasets(" << ndatasets << ")!=nrep(" << shape(timeDim) << ")*nslices(" << shape(sliceDim) << ")" << STD_endl;
    }

    // Insert slices into data according to their ordering
    Data<float,4> data(shape);
    ODINLOG(odinlog,normalDebug) << "mem(data): " << Profiler::get_memory_usage() << STD_endl;
    int irep=0;
    int islice=0;
    STD_map<double,int> sliceindexmap; // gets sorted according to temporal index
    for(imgit=imgmap.begin(); imgit!=imgmap.end(); ++imgit) {
      if(irep==0) sliceindexmap[imgit->first.acquisition_time]=islice;
      data(irep,islice,all,all)=imgit->second;
      result++;
      irep++;
      if(irep>=shape(timeDim)) {
        irep=0;
        islice++;
        if(islice>=shape(sliceDim)) islice=0;
      }
    }


    // set new acquisition start and duration
    prot.seqpars.set_AcquisitionStart(min_acqstart);
    ODINLOG(odinlog,normalDebug) << "min_acqstart: " << min_acqstart << " max_acqstart: " << max_acqstart << " dur: "<< max_acqstart - min_acqstart+ prot.seqpars.get_RepetitionTime()/1000 << STD_endl;
    double expdur=(max_acqstart - min_acqstart+ prot.seqpars.get_RepetitionTime()/1000 ) /60.0;
    ODINLOG(odinlog,normalDebug) << "ExpDuration(acqstart)=" << expdur << STD_endl;
    prot.seqpars.set_ExpDuration(expdur);

    unsigned int nslices=shape(sliceDim);

    // overwrite slice settings in protocol
    prot.geometry.set_nSlices(nslices);

    // store slice order in protocol for following steps
    ODINLOG(odinlog,normalDebug) << "nslices/sliceindexmap.size()=" << nslices << "/" << sliceindexmap.size() << STD_endl;
    if( nslices>1 && nslices<=sliceindexmap.size()) {
      LDRintArr sliceorder;
      sliceorder.redim(nslices); sliceorder.set_label("SliceOrder"); // separately instead of using constructor for GCC3.2
      STD_map<double,int>::const_iterator it=sliceindexmap.begin();
      for(unsigned int i=0; i<nslices; i++) {
        sliceorder[i]=it->second;
        ++it;
      }
      ODINLOG(odinlog,normalDebug) << "sliceorder=" << sliceorder.printbody() << STD_endl;
      prot.methpars.append_copy(sliceorder);
    }

    fvector slicevec=list2vector(slice_offsets);
    ODINLOG(odinlog,normalDebug) << "slicevec=" << slicevec.printbody() << STD_endl;

    if(slicevec.size()>1) {
      float slicedist=fabs(slicevec[1]-slicevec[0]);
      ODINLOG(odinlog,normalDebug) << "slicedist=" << slicedist << STD_endl;
      prot.geometry.set_sliceDistance(slicedist);
    }
    if(slicevec.size()>0) prot.geometry.set_offset(sliceDirection,0.5*(slicevec[0]+slicevec[slicevec.size()-1]));

    pdmap[prot].reference(data); // Finally, add data to protocol-data map

    imgmap.clear(); // save memory
    protindex++;
  } // end iterating over protocols

  ODINLOG(odinlog,normalDebug) << "mem(comb): " << Profiler::get_memory_usage() << STD_endl;

  return result;
}

//////////////////////////////////////////////////////////////////////////

int FileIO::autowrite(const ProtocolDataMap& pdmap, const STD_string& filename, const FileWriteOpts& opts) {
  Log<FileIO> odinlog("FileIO","autowrite");

  if(filename=="") {
    ODINLOG(odinlog,errorLog) << "Empty file name" << STD_endl;
    return -1;
  }

  FileFormatCreator ffc;

  ODINLOG(odinlog,normalDebug) << "filename:opts.format:opts.datatype=" << filename << ":" << opts.format << ":" << opts.datatype << STD_endl;

  FileFormat* ff=FileFormat::get_format(filename,opts.format);
  if(!ff) {
    FileFormat::format_error(filename);
    return -1;
  }


  if(opts.wprot!="") {
    svector protfnames=FileFormat::create_unique_filenames(opts.wprot, pdmap, opts.fnamepar);
    unsigned int ifile=0;
    for(ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {
      ODINLOG(odinlog,loglevel()) << "Storing protocol in file " << protfnames[ifile] << STD_endl;
      pdit->first.write(protfnames[ifile]);
      ifile++;
    }
  }

  // local copy to avoid recursion
  FileWriteOpts opts_copy(opts);
  opts_copy.split=false;


  // TODO: check and correct protocol and data for consistency (nrep, nslices, ny, nx)


  ODINLOG(odinlog,loglevel()) << "Writing format " << ff->description() << STD_endl;
  int result=0;
  if(opts.split) {

    svector splitfnames=FileFormat::create_unique_filenames(filename, pdmap, opts.fnamepar);

    unsigned int ifile=0;
    for(ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {

      STD_string onefilename=splitfnames[ifile];
      ODINLOG(odinlog,normalDebug) << "onefilename=" << onefilename << STD_endl;

      ProtocolDataMap mDummy;
      mDummy[pdit->first].reference(pdit->second);
      int writeresult=ff->write(mDummy, onefilename, opts_copy);
      if(writeresult<0) return -1;
      result+=writeresult;
      ODINLOG(odinlog,loglevel()) << "Wrote dataset to file " << onefilename << STD_endl;
      ifile++;
    }

  } else {
    result=ff->write(pdmap, filename, opts_copy);
    if(result<0) return -1;
    ODINLOG(odinlog,loglevel()) << "Wrote " << pdmap.size() << " dataset(s) to file " << filename << STD_endl;
  }

  return result;
}

//////////////////////////////////////////////////////////////


svector FileIO::autoformats() {
  Log<FileIO> odinlog("FileIO","autoread");
  FileFormatCreator ffc; // for some reason, we needed a named object here instead of a plain constructor
  return FileFormat::possible_formats();
}

//////////////////////////////////////////////////////////////

STD_string FileIO::autoformats_str(const STD_string& indent) {
  FileFormatCreator ffc; // for some reason, we needed a named object here instead of a plain constructor
  return FileFormat::formats_str(indent);
}

//////////////////////////////////////////////////////////////

bool FileIOTrace::do_trace=true;

STD_string FileFormat::select_write_datatype(const Protocol &prot,const FileWriteOpts& opts){
  if(opts.datatype != AUTOTDATAYPESTR)
    return opts.datatype;
  else
    return prot.system.get_data_type();
}

///////////////////////////////////////////////////////////////////////////////
// Unit Test


#ifndef NO_UNIT_TEST


////////// General tests ///////////////////////////////

class FileIOTest : public UnitTest {

 public:
  FileIOTest() : UnitTest("FileIO") {}

 private:
  bool check() const {
    Log<UnitTest> odinlog(this,"check");

    Range all=Range::all();

    FileIO::set_trace_status(false); // disable logging to console

    // Testing whether files are sorted alphabetically according to their file name
    STD_string tmpdir(tempfile());
    if(createdir(tmpdir.c_str())) return false;
    unsigned int testsize=16;
    Data<float,4> testslice(1,1,testsize,testsize);

    int nfiles=22;
    for(int i=0; i<nfiles; i++) {
      testslice=float(i);
      if(testslice.autowrite(tmpdir+SEPARATOR_STR+itos(i,nfiles)+".xml")<0) return false;
    }

    Data<float,4> testdirarr;
    if(testdirarr.autoread(tmpdir)<0) return false;


    TinyVector<int,4> expected_shape(nfiles,1,testsize,testsize);

    if(testdirarr.shape()!=expected_shape) {
      ODINLOG(odinlog,errorLog) << "testdirarr.shape()=" << testdirarr.shape() << ", but expected " << expected_shape << STD_endl;
      return false;
    }

    for(int i=0; i<nfiles; i++) {
      float meanval=mean(testdirarr(i,0,all,all));
      if(fabs(meanval-float(i))>1.0e-3) {
        ODINLOG(odinlog,errorLog) << "meanval(" << i << ")=" << meanval << ", but expected " << float(i) << STD_endl;
        return false;
      }
    }


    // Testing reading of complex data
    STD_string tmpfile(tempfile()+".float");
    ComplexData<1> cdata(testsize);
    cdata=STD_complex(0.0,1.0);
    if(cdata.write(tmpfile)<0) return false;

    FileReadOpts opts;
    Data<float,1> fdata;

    STD_map<STD_string,float> expected;
    expected["abs"]=1.0;
    expected["pha"]=0.5*PII;
    expected["real"]=0.0;
    expected["imag"]=1.0;

    for(STD_map<STD_string,float>::const_iterator it=expected.begin(); it!=expected.end(); ++it) {
      opts.cplx=it->first;
      if(fdata.autoread(tmpfile, opts)<0) return false;
      if(fdata.size()!=testsize) {
        ODINLOG(odinlog,errorLog) << "reading complex raw: size mismatch" << STD_endl;
        return false;
      }
      float meanval=mean(fdata);
      float expected=it->second;
      if(fabs(meanval-expected)>1.0e-3) {
        ODINLOG(odinlog,errorLog) << "reading complex raw: mean(" << opts.cplx.operator STD_string() << ")=" << meanval << ", but expected " << expected << STD_endl;
        return false;
      }
    }


    // Testing write and read of multi-dataset DICOM data
    FileIO::ProtocolDataMap pdmapout;
    Protocol prot1;
    Protocol prot2;

    prot1.seqpars.set_AcquisitionStart(12.0);
    prot1.seqpars.set_RepetitionTime(10.0);
    prot1.seqpars.set_NumOfRepetitions(7);
    prot1.geometry.set_nSlices(3);
    prot1.study.set_Series("series11",11);

    prot2.seqpars.set_AcquisitionStart(5423.0);
    prot2.seqpars.set_RepetitionTime(50.0);
    prot2.seqpars.set_NumOfRepetitions(4);
    prot2.geometry.set_nSlices(9);
    prot2.study.set_Series("series22",22);


    Data<float,4> outdata1(prot1.seqpars.get_NumOfRepetitions(),prot1.geometry.get_nSlices(),56,78); outdata1=111.1;
    Data<float,4> outdata2(prot2.seqpars.get_NumOfRepetitions(),prot2.geometry.get_nSlices(),64,96); outdata2=222.2;

    pdmapout[prot1].reference(outdata1);
    pdmapout[prot2].reference(outdata2);

    FileWriteOpts writeopts;
    writeopts.format="dcm";

    LDRfileName writefname=tempfile();
    ODINLOG(odinlog,normalDebug) << "writefname(dcm)=" << writefname << STD_endl;

    if(FileIO::autowrite(pdmapout, writefname, writeopts)<0) return false;


    FileIO::ProtocolDataMap pdmapin;
    FileReadOpts readopts;
    readopts.format="dcm";
    STD_string readfname=writefname.get_dirname()+SEPARATOR_STR+writefname.get_basename_nosuffix()+STD_string("_dcm");
    ODINLOG(odinlog,normalDebug) << "readfname(dcm)=" << readfname << STD_endl;

    Protocol dummy_protocol_template;
    if(FileIO::autoread(pdmapin, readfname, readopts, dummy_protocol_template)<0) return false;

    if(pdmapout.size()!=pdmapin.size()) {
      ODINLOG(odinlog,errorLog) << "pdmap size mismatch, readfname=" << readfname << ", pdmapout(" << pdmapout.size() << ")!=pdmapin(" << pdmapin.size() << ")" << STD_endl;
      return false;
    }

    FileIO::ProtocolDataMap::const_iterator pdit=pdmapin.begin();
    Data<float,4> indata1(pdit->second);
    if(indata1.shape() != outdata1.shape()) {
      ODINLOG(odinlog,errorLog) << "data1 shape mismatch, readfname=" << readfname << STD_endl;
      return false;
    }

    ++pdit;
    Data<float,4> indata2(pdit->second);
    if(indata2.shape() != outdata2.shape()) {
      ODINLOG(odinlog,errorLog) << "data2 shape mismatch, readfname=" << readfname << STD_endl;
      return false;
    }

    return true;
  }
};


////////// Format-specific tests ///////////////////////////////


void create_fileio_testarr(Data<float,4>& testarr, const TinyVector<int,4>& shape) {
  testarr.resize(shape); testarr=0.0;
  for(unsigned int i=0; i<testarr.numElements(); i++) {
    TinyVector<int,4> indexvec=testarr.create_index(i);
    for(int j=0; j<testarr.rank(); j++) {
      testarr(indexvec)+=pow(-1.0,i)*indexvec(j)*pow(10.0,j-2);
    }
  }
}

STD_string label4unittest(const STD_string& suff, const STD_string& uniquesuff, const STD_string& datatype) {
  STD_string result("FileIO "+suff);
  if(uniquesuff!="" || datatype!="") {
    result+="( ";
    if(uniquesuff!="") result+=uniquesuff+" ";
    if(datatype!="") result+=datatype+" ";
    result+=")";
  }
  return result;
}

template <int XSize, int YSize, typename CompareType, bool OnlyBitmap, bool HasRepetitions, bool ReadDir, bool HasOrientation, bool HasSliceThick>
class FileIOFormatTest : public UnitTest {

 public:
  FileIOFormatTest(const STD_string& suff, const STD_string& uniquesuff="", const STD_string& datatype="")
    : UnitTest(label4unittest(suff,uniquesuff,datatype).c_str()), suffix(suff), uniquesuffix(uniquesuff), dtype(datatype) {}

 private:

  STD_string suffix;
  STD_string uniquesuffix;
  STD_string dtype;

  bool compare_arrays(const STD_string& test, const Data<float,4>& a1, const Data<CompareType,4>& a2) const {
    Log<UnitTest> odinlog(this,"compare_arrays");
    if(a1.shape()!=a2.shape()) {
      ODINLOG(odinlog,errorLog) << test << " failed, shape mismatch:" << STD_endl;
      ODINLOG(odinlog,errorLog) << a1.shape() << " != " << a2.shape() << STD_endl;
      return false;
    }

    Data<CompareType,4> a1copy;
    a1.convert_to(a1copy);

    for(unsigned int i=0; i<a1.numElements(); i++) {
      TinyVector<int,4> indexvec=a1.create_index(i);
      if(a1copy(indexvec)!=a2(indexvec)) {
        ODINLOG(odinlog,errorLog) << test << " failed, value mismatch at index " << indexvec << STD_endl;
        ODINLOG(odinlog,errorLog) << a1copy(indexvec) << " != " << a2(indexvec) << STD_endl;
        return false;
      }
    }
    return true;
  }



  bool check() const {
    Log<UnitTest> odinlog(this,"check");

    FileIO::set_trace_status(false); // disable logging to console

    STD_list<TinyVector<int,4> > shapelst;
    shapelst.push_back(TinyVector<int,4>(1,1,YSize,XSize));
    if(!OnlyBitmap) {
      shapelst.push_back(TinyVector<int,4>(1,4,YSize,XSize));
      if(HasRepetitions) {
        shapelst.push_back(TinyVector<int,4>(3,4,YSize,XSize));
        shapelst.push_back(TinyVector<int,4>(3,1,YSize,XSize));
      }
    }

    for(STD_list<TinyVector<int,4> >::const_iterator it=shapelst.begin(); it!=shapelst.end(); ++it) {

      FileReadOpts readopts;
      FileWriteOpts writeopts;

      STD_string prefix(tempfile()); // separate name for each shape
      STD_string writefname=prefix+"."+suffix;
      STD_string readfname=writefname;
      if(ReadDir) {
        readfname=prefix+"_"+suffix;
        readopts.format=suffix;
      }

      if(uniquesuffix!="") {
        readopts.format=uniquesuffix;
        writeopts.format=uniquesuffix;
      }

      if(dtype!="") {
        writeopts.datatype=dtype;
      }

      Data<float,4> testarr;
      create_fileio_testarr(testarr,*it);


      Data<CompareType,4> readdata;

      ODINLOG(odinlog,normalDebug) << "writefname:writeopts.format:writeopts.datatype=" << writefname << ":" << writeopts.format << ":" << writeopts.datatype << STD_endl;

      // simple write/read
      if(testarr.autowrite(writefname, writeopts)<0) {
        ODINLOG(odinlog,errorLog) << "simple autowrite failed" << STD_endl;
        return false;
      }
      if(readdata.autoread(readfname, readopts)<0) {
        ODINLOG(odinlog,errorLog) << "simple autoread failed" << STD_endl;
        return false;
      }
      if(!compare_arrays("autowrite/autoread("+readfname+")", testarr, readdata)) return false;

      // Testing whether geometry info is preserved
      if(!OnlyBitmap) {
        Protocol prot;
        Geometry& geo=prot.geometry;
        if(HasOrientation) {
          geo.set_orientation(-66.7, 78.2, -124.7);
          geo.set_offset(readDirection,22.7);
          geo.set_offset(phaseDirection,-5.9);
          geo.set_offset(sliceDirection,99.9);
        }
        geo.set_FOV(readDirection,192.6);
        geo.set_FOV(phaseDirection,200.2);
        geo.set_nSlices(testarr.extent(1));
        geo.set_sliceDistance(6.1);
        if(HasSliceThick) {
          geo.set_sliceThickness(3.2);
        } else {
          geo.set_sliceThickness(6.1); // the same as slice distance
        }
        if(testarr.autowrite(writefname, writeopts, &prot)<0) {
          ODINLOG(odinlog,errorLog) << "autowrite with protocol failed" << STD_endl;
          return false;
        }
        Protocol readprot;
        if(readdata.autoread(readfname, readopts, &readprot)<0) {
          ODINLOG(odinlog,errorLog) << "autoread with protocol failed" << STD_endl;
          return false;
        }
        if(!compare_arrays("autowrite/autoread+geo("+readfname+")", testarr, readdata)) return false;
        Protocol protcopy(prot); // Create copy of protocol to use limited accuracy comparison of Protocol
        protcopy.geometry=readprot.geometry; // only geometry may differ
        if(!(prot==protcopy)) { // use limited accuracy comparison of Protocol
          ODINLOG(odinlog,errorLog) << "autowrite/autoread(geo)" << (*it) << " failed: prot.geometry=" << prot.geometry << "readprot.geometry=" << protcopy.geometry << STD_endl;
          return false;
        }
      }
    }

    return true;
  }

};


///////////////////////////////////////////////////////////////////////


void alloc_FileIOTest() { // create test instances

  new FileIOTest;

  new FileIOFormatTest<7, 13, float, true,false,false,false,false>("dat");


#ifdef VTKSUPPORT
  new FileIOFormatTest<7, 13, u8bit, false,false,false,false,false>("vtk");
#endif

#ifdef PNGSUPPORT
  new FileIOFormatTest<7, 13, u8bit, true,false,false,false,false>("png");
#endif

#ifdef HAVE_LIBZ
#ifndef STREAM_REPLACEMENT
  new FileIOFormatTest<7, 13, double, false,true,false,true,true>("jdx.gz");
  new FileIOFormatTest<7, 13, float, false,true,false,true,true>("xml.gz");
#endif
#endif

#ifdef DICOMSUPPORT
  new FileIOFormatTest<16, 16, u16bit, false,true,true,true,true>("dcm"); // needs square matrix size
#endif
#ifdef NIFTISUPPORT
// Disable these tests for now, since they fail when building with cowbuilder,
// but pass whenever building in a clean chroot. Need to have a closer look.
//  new FileIOFormatTest<7, 13, double, false,true,false,true,false>("nii","","float");
//  new FileIOFormatTest<7, 13, double, false,true,false,false,false>("hdr", "analyze");
#ifdef HAVE_LIBZ
#ifndef STREAM_REPLACEMENT
// Disable this test for now, since it fails when building with cowbuilder,
// but pass whenever building in a clean chroot. Need to have a closer look.
//  new FileIOFormatTest<7, 13, double, false,true,false,true,false>("nii.gz");
#endif
#endif
#endif

#ifdef ISMRMRDSUPPORT
  new FileIOFormatTest<7, 13, float, false,true,false,true,false>("h5");
#endif

  new FileIOFormatTest<7, 13, s16bit, false,false,false,false,false>("hdr", "interfile","s16bit");
  new FileIOFormatTest<7, 13, s16bit, false,false,false,false,false>("hdr", "interfile","float");

}
#endif