#include "fileio.h"
#include "utils.h"

#include <limits>

#ifdef DICOMSUPPORT

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

// Because we have 2 config.h files, one from ODIN and one from DCMTK,
//  Undefine doubly defined before include of DCMTK headers
#undef PACKAGE_BUGREPORT
#undef PACKAGE_NAME
#undef PACKAGE_STRING
#undef PACKAGE_TARNAME
#undef PACKAGE_VERSION
//#include <dcmtk/config/cfunix.h>
#include <dcmtk/dcmdata/dcfilefo.h>
#include <dcmtk/dcmdata/dcdatset.h>
#include <dcmtk/dcmdata/dcmetinf.h>
#include <dcmtk/dcmdata/dcdeftag.h>
#include <dcmtk/dcmdata/dcuid.h>
#include <dcmtk/dcmdata/dcvrus.h>
#include <dcmtk/dcmdata/dccodec.h>
#include <dcmtk/dcmdata/dcdict.h>
#include <dcmtk/dcmdata/dcddirif.h>
#include <dcmtk/dcmimgle/dcmimage.h>

// fix for use of possesive 's' across different versions of dcmtk
#ifdef DCM_PatientsName
#define DCM_PAT_NAME DCM_PatientsName
#else
#define DCM_PAT_NAME DCM_PatientName
#endif

#ifdef DCM_PatientsBirthDate
#define DCM_PAT_BIRTHDATE DCM_PatientsBirthDate
#else
#define DCM_PAT_BIRTHDATE DCM_PatientBirthDate
#endif

#ifdef DCM_PatientsSex
#define DCM_PAT_SEX DCM_PatientsSex
#else
#define DCM_PAT_SEX DCM_PatientSex
#endif

#ifdef DCM_PatientsWeight
#define DCM_PAT_WEIGHT DCM_PatientsWeight
#else
#define DCM_PAT_WEIGHT DCM_PatientWeight
#endif

#ifdef DCM_PatientsSize
#define DCM_PAT_SIZE DCM_PatientsSize
#else
#define DCM_PAT_SIZE DCM_PatientSize
#endif

#ifdef DCM_PatientsAge
#define DCM_PAT_AGE DCM_PatientsAge
#else
#define DCM_PAT_AGE DCM_PatientAge
#endif

#ifdef DCM_PerformingPhysiciansName
#define DCM_PERF_PHYS DCM_PerformingPhysiciansName
#else
#define DCM_PERF_PHYS DCM_PerformingPhysicianName
#endif

#ifdef DCM_ReferringPhysiciansName
#define DCM_REF_PHYS DCM_ReferringPhysiciansName
#else
#define DCM_REF_PHYS DCM_ReferringPhysicianName
#endif






///////////////////////////////////////
enum dicom_field_importance{optional=0,warning,error};
bool check_status(const char* func, const char* call, const OFCondition& status, const dicom_field_importance severity=error) {
  Log<FileIO> odinlog("DicomFormat","check_status");
  if(status.bad()) {
    logPriority loglevel=noLog;
    switch(severity){
      case optional:  loglevel=normalDebug;break;
      case warning:   loglevel=warningLog;break;
      case error:     loglevel=errorLog;break;
    }
    ODINLOG(odinlog,loglevel) << func << "(" << call << ")" << ": " << status.text() << STD_endl;
    return true;
  }
  return false;
}

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

bool check_dict(const char* func) {
  Log<FileIO> odinlog("DicomFormat",func);
  if (!dcmDataDict.isDictionaryLoaded()) {
    ODINLOG(odinlog,errorLog) << "No data dictionary loaded, check environment variable " << DCM_DICT_ENVIRONMENT_VARIABLE << STD_endl;
    svector dirtoks=tokens(DCM_DICT_DEFAULT_PATH, ':');
    for(unsigned int itok=0; itok<dirtoks.size(); itok++) {
      if(filesize(dirtoks[itok].c_str())<0) {
        ODINLOG(odinlog,errorLog) << "Dictionary file "<< dirtoks[itok] << " of the current dcmtk installation does not exist, please check local dcmtk configuration" << STD_endl;
      }
    }
    return true;
  }
  return false;
}

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

bool insert_uint16_hack(DcmDataset *dataset, const DcmTagKey &key, u16bit value) { // Hack to insert uint16 numbers on pre-3.5.4 dcmtk versions, code taken from dcmtk forum
// the preprocessor directive  of the following line does not work with DCMTK 3.6.0 where OFFIS_DCMTK_VERSION_NUMBER is a string
//#if OFFIS_DCMTK_VERSION_NUMBER > 353
  dataset->putAndInsertUint16(key,value);
//#else
//  DcmElement *us = new DcmUnsignedShort(DcmTag(key, EVR_US));
//  if (us != NULL) {
//     if(check_status("insert_uint16_hack","putUint16",us->putUint16(value))) {return false;}
//     if(check_status("insert_uint16_hack","insert",dataset->insert(us, OFTrue /*replaceOld*/))) {return false;}
//  } else return false;
//#endif
  return true;
 }

template<typename T> void copy(const DiPixel *pix, Data<float,4> &dst,const TinyVector<int,4> &shape,int mosaic_mult){
  Log<FileIO> odinlog("DicomFormat","copy");

  T* src=(T*)pix->getData();

  const Range all=Range::all();

  if(shape(1)>1) { // Sort Siemens Mosaic data
    dst.resize(shape);
    Data<float,4> indata;
    TinyVector<int,4> inshape(mosaic_mult,shape(2),mosaic_mult,shape(3));
    convert_from_ptr(indata,src,inshape);

    for(int n=0; n<mosaic_mult; n++) {
      for(int m=0; m<mosaic_mult; m++) {
        const int slice=n*mosaic_mult+m;
        if(slice<shape(1)) dst(0,slice,all,all)=indata(n,all,m,all);
      }
    }

  } else {
    convert_from_ptr(dst,src,shape);
  }
}

template<typename BASE,typename DST> DST endian(const BASE *b){
  Log<FileIO> odinlog("DicomFormat","endian");
  DST ret=0;
#if __BYTE_ORDER == __LITTLE_ENDIAN
  for(short i=0;i<(short)sizeof(DST);i++) {
#elif __BYTE_ORDER == __BIG_ENDIAN
  for(short i=(short)sizeof(DST)-1;i>=0;i--) {
#else
#error __BYTE_ORDER undefined
#endif
    ODINLOG(odinlog,normalDebug) << "b[" << i << "]" << int(b[i]) << STD_endl;
    ret+=b[i]<<i*8;
  }
  return ret;
}


svector fetch_from_MR_CSA_Header(DcmElement* csaHeader, STD_string label){
  Log<FileIO> odinlog("DicomFormat","fetch_from_MR_CSA_Header");

  svector ret;

  unsigned int csalength=csaHeader->getLength();
  ODINLOG(odinlog,normalDebug) << "csaHeader->getLength()=" << csalength << STD_endl;
  if(!csalength) return ret;

  Uint8 *array;
  csaHeader->getUint8Array(array);
  for(STD_string::size_type pos=0;pos<=csalength;){
    STD_string dummy((char*)array+pos);
    STD_string::size_type found=dummy.find(label);
    if(found==STD_string::npos) {
      pos+=dummy.length()+1;
    } else {
      ODINLOG(odinlog,normalDebug) << "found >" << label << "< at pos=" << pos << STD_endl;
      pos+=found+0x40; //len of the name is 0x20 plus 0x20 "I dont know what the hell this is supposed to be"
      int markval=array[pos]; // supposed to be 0x01 at this position for usable values
      if(markval!=0x01) {
        ODINLOG(odinlog,normalDebug) << "markval=" << markval << STD_endl;
        break; // integrity check
      }
      /*Sint32 &vm=*((Sint32*)array+pos);*/pos+=sizeof(Sint32);
      /*char *vr=(char*)array+pos;*/pos+=0x4;
      /*Sint32 syngodt=endian<Uint8,Uint32>(array+pos);*/pos+=sizeof(Sint32);
      Sint32 nitems=endian<Uint8,Uint32>(array+pos);pos+=sizeof(Sint32);
      ODINLOG(odinlog,normalDebug) << "nitems=" << nitems << STD_endl;
      if(nitems) {
        pos+=sizeof(Sint32); //77
        Sint32 len=0;
        for(unsigned short n=0;n<nitems;n++){
          len=endian<Uint8,Uint32>(array+pos);pos+=sizeof(Sint32);//the length of this element
          pos+=3*sizeof(Sint32); //whatever
          if(!len)continue;
          const unsigned s=ret.size();ret.resize(s+1);
          ret[s]=STD_string((char*)array+pos);//store the text if the is some
          pos+=((len+sizeof(Sint32)-1)/sizeof(Sint32))*sizeof(Sint32);//increment pos by len aligned to sizeof(Sint32)
          if(pos>csalength) { // integrity check
             ODINLOG(odinlog,normalDebug) << "pos out of bounds" << STD_endl;
             break;
          }
        }
        ODINLOG(odinlog,normalDebug) << "Found \""+label+"\" in CSA header: \"" << ret.printbody() << "\"" << STD_endl;
      } else {
        ODINLOG(odinlog,normalDebug) << "Found empty \""+label+"\" in CSA header" << STD_endl;
      }
      return ret;
    }
  }
  ODINLOG(odinlog,normalDebug) << "Could not find usable field \""+label+"\" in CSA header" << STD_endl;

  return ret;
}

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

void timestr2seconds(const OFString& timestr, LONGEST_INT& sec, double& fract) {
  sec=0;
  fract=0;
  if(timestr.length()<13) return;
  //@todo backward compatibility missing (http://idlastro.gsfc.nasa.gov/idl_html_help/Value_Representations.html#TM)
  LONGEST_INT hh=atoi(timestr.substr(0,2).c_str());
  LONGEST_INT mm=atoi(timestr.substr(2,2).c_str());
  LONGEST_INT ss=atoi(timestr.substr(4,2).c_str());
  sec=ss + 60*mm + 3600*hh;
  fract=atof(timestr.substr(6,7).c_str());
}

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

struct DicomFormat : public FileFormat {
  STD_string description() const {return "DICOM";}
  svector suffix() const  {
    svector result; result.resize(4);
    result[0]="dcm";
    result[1]="mag";
    result[2]="ph";
    result[3]="ima";
    return result;
  }
  svector dialects() const {
    svector result; result.resize(1);
    result[0]="siemens";
    return result;
  }


  int read(Data<float,4>& data, const STD_string& filename, const FileReadOpts& opts, Protocol& prot) {
    Log<FileIO> odinlog("DicomFormat","read");

    OFString valstr;

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

    const char* func="DicomFormat::read";
    if(check_dict(func)) return -1;


    if(LDRfileName(filename).get_basename()=="DICOMDIR") return 0; // ignore DICOMDIR

#ifdef OFLOG_H
    OFLog::configure(OFLogger::WARN_LOG_LEVEL); // reset logging (in dcmtk 3.6.0 info logging is somehow enabled sometimes)
#endif

    DcmFileFormat fileformat;
    DcmDataset* dset=0;
    if(check_status(func,"loadFile",fileformat.loadFile(filename.c_str()))) {
      ODINLOG(odinlog,errorLog) << "Unable to get dataset" << STD_endl;
      return -1;
    }
    dset=fileformat.getDataset();
    if(!dset) {
      ODINLOG(odinlog,errorLog) << "fileformat.getDataset() failed" << STD_endl;
      return -1;
    }


    // The meta-info stuff not provided by dcmtk
    DcmMetaInfo* metainfo=fileformat.getMetaInfo();


    STD_string syntax;
    if(!check_status(func,"TransferSyntaxUID",metainfo->findAndGetOFStringArray(DCM_TransferSyntaxUID, valstr),warning)) syntax=valstr.c_str();
    ODINLOG(odinlog,normalDebug) << "syntax=" << syntax << STD_endl;
    if(syntax==UID_JPEG2000LosslessOnlyTransferSyntax || syntax==UID_JPEG2000TransferSyntax) {
      ODINLOG(odinlog,errorLog) << "Transfer syntax not supported in free version of dcmtk" << STD_endl;
      return -1;
    }


/*  change representation (does not work)
    if(syntax==UID_JPEGProcess14SV1TransferSyntax) {
      if(!dset->chooseRepresentation(EXS_LittleEndianExplicit, NULL).good()) {
        ODINLOG(odinlog,errorLog) << "Cannot change representation" << STD_endl;
        return -1;
      }
    }
*/


    DcmElement* csaHeader=0;
    if(dset->findAndGetElement(DcmTagKey(0x0029,0x1010),csaHeader).good()) {
      ODINLOG(odinlog,normalDebug) << "found Siemens-specific CSA Header" << STD_endl;

/*
      // write headers to file for debugging
      Uint8* array=0;
      csaHeader->getUint8Array(array);
      if(array) {
        FILE* fp=fopen("/tmp/csa", "w");
        fwrite(array, sizeof(Uint8), csaHeader->getLength(),fp);
        fclose(fp);
      }
*/
     }

     DcmElement* csaSeriesHeader=0;
     if(dset->findAndGetElement(DcmTagKey(0x0029,0x1020),csaSeriesHeader).good()) {
       ODINLOG(odinlog,normalDebug) << "found Siemens-specific CSA Series Header" << STD_endl;

       Uint8* serarray=0;
       csaSeriesHeader->getUint8Array(serarray);


      if(serarray) {

/*      // XProtocol stuff
        const char* findpattern="<XProtocol>";
        unsigned int findlength=STD_string(findpattern).length();

        STD_list<unsigned int> poslist;
        for(unsigned int ipos=0; ipos<=(csaSeriesHeader->getLength()-findlength); ipos++) {
          bool found=true;
          for(unsigned int ilength=0; ilength<findlength; ilength++) {
            if(serarray[ipos+ilength]!=findpattern[ilength]) {
              found=false;
              break;
            }
          }
          if(found) {
            poslist.push_back(ipos);
            ODINLOG(odinlog,normalDebug) << "found " << findpattern << " at pos=" << ipos << STD_endl;
          }
        }
*/

/*
          // write headers to file for debugging
          FILE* fp=fopen("/tmp/csaser", "w");
          fwrite(serarray, sizeof(Uint8), csaSeriesHeader->getLength(),fp);
          fclose(fp);
*/
      }
    }




    DicomImage img(dset, EXS_Unknown, CIF_IgnoreModalityTransformation); // apply modality transform manually to get higher floating point accuracy

    EI_Status status=img.getStatus();
    if(status!=EIS_Normal) {
      const char* errmsg="Unknown Error";
      if(status==EIS_NoDataDictionary)  errmsg="NoDataDictionary";
      if(status==EIS_InvalidDocument)   errmsg="InvalidDocument";
      if(status==EIS_MissingAttribute)  errmsg="MissingAttribute";
      if(status==EIS_InvalidValue)      errmsg="InvalidValue";
      if(status==EIS_NotSupportedValue) errmsg="NotSupportedValue";
      if(status==EIS_MemoryFailure)     errmsg="MemoryFailure";
      if(status==EIS_InvalidImage)      errmsg="InvalidImage";
      if(status==EIS_OtherError)        errmsg="OtherError";
      ODINLOG(odinlog,errorLog) << "Unable to load DICOM image " << filename << " - " << errmsg << STD_endl;
      return -1;
    }


    TinyVector<int,4> shape;
    shape=1; // default

    shape(2)=img.getHeight();
    shape(3)=img.getWidth();
    ODINLOG(odinlog,normalDebug) << "shape=" << shape << STD_endl;


    STD_string image_type;
    if(!check_status(func,"ImageType",dset->findAndGetOFStringArray(DCM_ImageType, valstr),warning)) image_type=valstr.c_str();
    ODINLOG(odinlog,normalDebug) << "image_type=" << image_type << STD_endl;

    // Hack for Siemens Mosaic
    int mosaic_mult=1;
    Sint16 mosaic_slices=1;

    //try to find 0x0019,100a (images per mosaic)
    if(image_type.find("MOSAIC")!=STD_string::npos) {
      if(dset->findAndGetSint16(DcmTagKey(0x0019,0x100a),mosaic_slices).good()){
        ODINLOG(odinlog,normalDebug) << "Found \"NumberOfImagesInMosaic\" in 0x0019,0x100a: " << mosaic_slices << STD_endl;
      } else {
        if(csaHeader) {
          //othwise find NumberOfImagesInMosaic in Siemens shadow Header 0x0029,0x1010
          svector csa_mosaic=fetch_from_MR_CSA_Header(csaHeader,"NumberOfImagesInMosaic");
          if(csa_mosaic.size()) {
            mosaic_slices=atoi(csa_mosaic[0].c_str());
            ODINLOG(odinlog,normalDebug) << "Found \"NumberOfImagesInMosaic\" in CSADataInfo: " << mosaic_slices << STD_endl;
          }
        }
      }
      if(mosaic_slices>1) {
        mosaic_mult=(int)ceil(sqrt((double)mosaic_slices));
        ODINLOG(odinlog,normalDebug) << "mosaic_mult=" << mosaic_mult << STD_endl;
        ODINLOG(odinlog,normalDebug) << "mosaic_slices=" << mosaic_slices << STD_endl;
        shape(3)/=mosaic_mult;
        shape(2)/=mosaic_mult;
        shape(1)=mosaic_slices;
      } else {
        ODINLOG(odinlog,warningLog) << "\"NumberOfImagesInMosaic\" not available. Mosaic won't be split up" << STD_endl;
        shape(1)=1;
      }
    }


    //Get Diffusion Data
    Float64 Bvalue;
    OFCondition foundB=dset->findAndGetFloat64(DCM_DiffusionBValue,Bvalue); //in case someone actually used the right Tag
    if(foundB.bad() && opts.dialect=="siemens"){//fallback for siemens
      OFString buff;
      foundB=dset->findAndGetOFString(DcmTagKey(0x0019,0x100c), buff);
      Bvalue=atof(buff.c_str());
    }
    //@todo fallback for GE/Philips
    if(foundB.good()){//we got the b_value
      LDRtriple bvector(0,0,0,"Diffusion_bVector");
      if(Bvalue){
        const Float64 *dir;unsigned long count;
        OFCondition foundD=dset->findAndGetFloat64Array(DCM_DiffusionGradientOrientation,dir,&count); //in case someone actually used the right Tag
        if(foundD.bad() && opts.dialect=="siemens") foundD=dset->findAndGetFloat64Array( DcmTagKey(0x0019,0x100e),dir,&count);//siemens fallback
        if(foundD.good()) {
          for(unsigned long i=0;i<count && i<3;i++)
            bvector[i]=dir[i]*Bvalue;
          ODINLOG(odinlog,normalDebug) <<  "Found diffusion b-Value: [" << bvector.printbody() << "] in >" << filename << "<"<< STD_endl;
        } else
          ODINLOG(odinlog,warningLog) <<  "no diffusion direction found in >" << filename << "<"<< STD_endl;
      }
      prot.methpars.append_copy(bvector);
    }

    //detect pixel datatype and copy pixel data
    //@todo: Transpose data and geometry in plane according to DCM_InPlanePhaseEncodingDirection
    if(!img.isMonochrome()){
      ODINLOG(odinlog,errorLog) <<  "only gray value images are supportet"<< STD_endl;
      return -1;
    }

    const DiPixel *pix=img.getInterData();
    switch(pix->getRepresentation()){
      case EPR_Uint8:
        copy<u8bit>  (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((u8bit)0));
        break;
      case EPR_Uint16:
        copy<u16bit> (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((u16bit)0));
        break;
      case EPR_Uint32:
        copy<u32bit> (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((u32bit)0));
        break;
      case EPR_Sint8:
        copy<s8bit>  (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((s8bit)0));
        break;
      case EPR_Sint16:
        copy<s16bit> (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((s16bit)0));
        break;
      case EPR_Sint32:
        copy<s32bit> (pix,data,shape,mosaic_mult);
        prot.system.set_data_type(TypeTraits::type2label((s32bit)0));
        break;
      default: ODINLOG(odinlog,errorLog) <<  "pixel representation not supported"<< STD_endl;
    }

#ifdef ODIN_DEBUG
    double min,max;
    img.getMinMaxValues(min,max);
    ODINLOG(odinlog,normalDebug) << "image range is " << min << " - " << max << STD_endl;
    ODINLOG(odinlog,normalDebug) << "datatype set to " << prot.system.get_data_type() << STD_endl;
#endif //ODIN_DEBUG

    // rescale data, if necessary
    float intercept=0;
    float slope=1.0;
    if(!check_status(func,"RescaleIntercept",dset->findAndGetOFString(DCM_RescaleIntercept, valstr),optional)) {
      intercept=atof(valstr.c_str());
    }
    if(!check_status(func,"RescaleSlope",dset->findAndGetOFString(DCM_RescaleSlope, valstr),optional)) {
      slope=atof(valstr.c_str());
    }
    if(intercept!=0.0 || slope!=1.0) {
      ODINLOG(odinlog,normalDebug) << "rescaling data with intercept/slope=" << intercept << "/" << slope << STD_endl;
      data=data*slope+intercept;
      prot.system.set_data_type(TypeTraits::type2label((float)0)); // use float to preserve correct values
    }



    STD_string unknown=ODIN_UNKNOWN_STR;

    // Date/time stuff
    STD_string scandate;
    STD_string scantime;
    if(!check_status(func,"SeriesDate",dset->findAndGetOFString(DCM_SeriesDate, valstr),warning)) {
      scandate=valstr.c_str();
    }
    OFString seriestimestr;
    if(!check_status(func,"SeriesTime",dset->findAndGetOFString(DCM_SeriesTime, seriestimestr),warning)) {
      svector toks=tokens(seriestimestr.c_str(), '.');
      if(toks.size()) scantime=toks[0];
    }
    prot.study.set_DateTime(scandate,scantime);


    // Get patient stuff
    STD_string pname(unknown);
    if(!check_status(func,"PatientsName",dset->findAndGetOFString(DCM_PAT_NAME, valstr),warning)) {
      pname=valstr.c_str();
      pname=replaceStr(pname,"^"," ");
      pname=valstr.c_str();
    }
    STD_string pid(unknown);
    if(!check_status(func,"PatientID",dset->findAndGetOFString(DCM_PatientID, valstr),warning)) {
      pid=valstr.c_str();
    }
    STD_string pdbirth(ODIN_DATE_LENGTH,'0');
    if(!check_status(func,"PatientsBirthDate",dset->findAndGetOFString(DCM_PAT_BIRTHDATE, valstr),warning)) {
      pdbirth=valstr.c_str();
    }
    char psex='M';
    if(!check_status(func,"PatientsSex",dset->findAndGetOFString(DCM_PAT_SEX, valstr),warning)) {
      if(valstr.length()) psex=valstr[0];
    }
    float pweight=50.0;
    if(!check_status(func,"PatientsWeight",dset->findAndGetOFString(DCM_PAT_WEIGHT, valstr),warning)) {
      if(valstr.length()) pweight=atof(valstr.c_str());
    }
    float psize=2000.0;
    if(!check_status(func,"PatientsSize",dset->findAndGetOFString(DCM_PAT_SIZE, valstr),warning)) {
      if(valstr.length()) psize=1000.0*atof(valstr.c_str());
    }
    prot.study.set_Patient(pid, pname, pdbirth, psex, pweight, psize);


    // Get study stuff
    STD_string studydescr(unknown);
    if(!check_status(func,"StudyDescription",dset->findAndGetOFString(DCM_StudyDescription, valstr),warning)) {
      studydescr=valstr.c_str();
    }
    STD_string physician(unknown);
    if(!check_status(func,"PerformingPhysiciansName",dset->findAndGetOFString(DCM_PERF_PHYS, valstr),optional)) {
      physician=valstr.c_str();
    }
    prot.study.set_Context(studydescr,physician);


    // Get series stuff
    STD_string seriesdescr(unknown);
    if(!check_status(func,"SeriesDescription",dset->findAndGetOFString(DCM_SeriesDescription, valstr),warning)) {
      seriesdescr=valstr.c_str();
    }
    int seriesno=1;
    if(!check_status(func,"SeriesNumber",dset->findAndGetOFString(DCM_SeriesNumber, valstr),warning)) {
      seriesno=atoi(valstr.c_str());
    }
    if(opts.dialect=="siemens") {
      // Quick hack to produce different datasets for different channels on Siemens
      if(dset->findAndGetOFString(DcmTagKey (0x0051, 0x100f), valstr).good()) {
        seriesdescr+="_"+STD_string(valstr.c_str());
      }
    }
    prot.study.set_Series(seriesdescr,seriesno);


    // get Pixel Spacing and FOV
    float FOVread=200.0,psize_read;
    if(!check_status(func,"PixelSpacing(read)",dset->findAndGetOFString(DCM_PixelSpacing, valstr, 0), warning)) {
      psize_read=atof(valstr.c_str());
      FOVread=psize_read*shape(3);
      prot.geometry.set_FOV(readDirection, FOVread);
    }else
      psize_read=FOVread/shape(3);

    float FOVphase=200.0,psize_phase;
    if(!check_status(func,"PixelSpacing(phase)",dset->findAndGetOFString(DCM_PixelSpacing, valstr, 1), warning)) {
      psize_phase=atof(valstr.c_str());
      FOVphase=psize_phase*shape(2);
      prot.geometry.set_FOV(phaseDirection, FOVphase);
    }else
      psize_phase=FOVphase/shape(2);


    // Get the slice orientation, i.e. the vectors of read and phase direction
    dvector readvec(3);  readvec[0]=1.0;
    dvector phasevec(3); phasevec[1]=1.0;
    dvector slicevec(3); slicevec[2]=1.0;

    // The following does not work with dcmtk since we can only access the first 3 elements of ImageOrientationPatient
    dvector offset(3); offset = -0.5*(FOVread-psize_read)*readvec - 0.5*(FOVphase-psize_phase)*phasevec;
    for(int i=0; i<3; i++) {
      if(!check_status(func,("ImageOrientationPatient("+itos(i)+")").c_str(),dset->findAndGetOFString(DCM_ImageOrientationPatient, valstr,i))) {
        readvec[i]=atof(valstr.c_str());
      }
      // Does not work in dcmtk
      if(!check_status(func,("ImageOrientationPatient("+itos(3+i)+")").c_str(),dset->findAndGetOFString(DCM_ImageOrientationPatient, valstr,3+i))) {
        phasevec[i]=atof(valstr.c_str());
      }
      if(!check_status(func,("ImagePositionPatient("+itos(i)+")").c_str(),dset->findAndGetOFString(DCM_ImagePositionPatient, valstr,i))) {
        offset[i]=atof(valstr.c_str());
      }
    }

    //remove the additional mosaic offset
    //eg. if there is a 10x10 Mosaic, substract the half size of 9 Images from the offset
    for(int i=0;i<3;i++) {
      offset[i]+=(mosaic_mult-1)*FOVread*.5*readvec[i] + (mosaic_mult-1)*FOVphase*.5*phasevec[i];
    }

    //Following DICOM standard ImagePositionPatient is the _middle_ of the first vector so we would
    //use (FOV-voxelsize)/2 to get the center.
    dvector center_offset= offset + 0.5*(FOVread-psize_read)*readvec + 0.5*(FOVphase-psize_phase)*phasevec;

    //But siemens seems to store the outer edge of the first voxel. So...
    if(opts.dialect=="siemens") center_offset= offset + 0.5*FOVread*readvec + 0.5*FOVphase*phasevec;

    dvector transform(3); transform[0]=-1.0; transform[1]=-1.0; transform[2]=1.0; // default: x- and y-axis are inverted when converting from DICOM coord system to ODIN coord system (as defined in geometry.h)

    if(!check_status(func,"PatientPosition",dset->findAndGetOFString(DCM_PatientPosition, valstr),optional)) {
      bool known_position=false;
      if(valstr=="HFP") {
        transform[0]=-transform[0];
        transform[1]=-transform[1];
        known_position=true;
      }
      if(valstr=="HFS") { // use default
        known_position=true;
      }
      if(!known_position) {
        ODINLOG(odinlog,normalDebug) << "Unknown PatientPosition=" << valstr << STD_endl;
      }
    }

    readvec*=transform;
    phasevec*=transform;
    center_offset*=transform;


    slicevec=Data<double,1>(vector_product(Data<double,1>(readvec), Data<double,1>(phasevec))); // has correct handness for ODIN

    ODINLOG(odinlog,normalDebug) << "psize_read/phase=" << psize_read << "/" << psize_phase << STD_endl;
    ODINLOG(odinlog,normalDebug) << "offset=" << offset.printbody() << STD_endl;
    ODINLOG(odinlog,normalDebug) << "center_offset=" << center_offset.printbody() << STD_endl;
    ODINLOG(odinlog,normalDebug) << "readvec=" << readvec.printbody() << STD_endl;
    ODINLOG(odinlog,normalDebug) << "phasevec=" << phasevec.printbody() << STD_endl;
    ODINLOG(odinlog,normalDebug) << "slicevec=" << slicevec.printbody() << STD_endl;

    prot.geometry.set_orientation_and_offset(readvec, phasevec, slicevec, center_offset);

    if(!check_status(func,"SliceThickness",dset->findAndGetOFString(DCM_SliceThickness, valstr),warning)) {
      prot.geometry.set_sliceThickness(atof(valstr.c_str()));
    }

    float space=0;
    if(!check_status(func,"SpacingBetweenSlices",dset->findAndGetOFString(DCM_SpacingBetweenSlices, valstr),optional)) {
      space+=atof(valstr.c_str());
      prot.geometry.set_sliceDistance(space);
    }

    // Calculate slice offset, do not use DCM_SliceLocation since it is based on variable handness
    // also take in account, that for mosaic the slice position is given for the first image of the mosaic
    float sliceproj=(slicevec*center_offset).sum();
    prot.geometry.set_offset(sliceDirection,sliceproj+(mosaic_slices-1)*space/2);


    STD_string modality("MR");
    if(!check_status(func,"Modality",dset->findAndGetOFString(DCM_Modality, valstr),warning)) {
      modality=valstr.c_str();
    }
    
    if(modality=="MR") { // avoid useles error messages with non-MR modalities
    
      // Get system stuff
      if(!check_status(func,"ImagingFrequency",dset->findAndGetOFString(DCM_ImagingFrequency, valstr),warning)) {
        prot.system.set_B0_from_freq(1000.0*atof(valstr.c_str()));
      }
      if(!check_status(func,"TransmitCoilName",dset->findAndGetOFString(DCM_TransmitCoilName, valstr),warning)) {
        prot.system.set_transmit_coil_name(valstr.c_str());
      }

      if(csaHeader) {
        svector imacoils=fetch_from_MR_CSA_Header(csaHeader,"ImaCoilString");
        prot.system.set_receive_coil_name(imacoils.printbody());
      }


      // Sequence stuff
      if(!check_status(func,"PixelBandwidth",dset->findAndGetOFString(DCM_PixelBandwidth, valstr),warning)) {
        prot.seqpars.set_AcqSweepWidth(0.001*atof(valstr.c_str())*shape(3));
      }
      if(!check_status(func,"RepetitionTime",dset->findAndGetOFString(DCM_RepetitionTime, valstr),warning)) {
        prot.seqpars.set_RepetitionTime(atof(valstr.c_str()));
      }
      if(!check_status(func,"EchoTime",dset->findAndGetOFString(DCM_EchoTime, valstr),warning)) {
        prot.seqpars.set_EchoTime(atof(valstr.c_str()));
      }
      if(!check_status(func,"FlipAngle",dset->findAndGetOFString(DCM_FlipAngle, valstr),warning)) {
        prot.seqpars.set_FlipAngle(atof(valstr.c_str()));
      }

      prot.seqpars.set_PhysioTrigger(dset->findAndGetOFString(DCM_TriggerTime, valstr).good()); // presence of TriggerTime indicates gating, at least on Siemens

      if(!check_status(func,"SequenceName",dset->findAndGetOFString(DCM_SequenceName, valstr),warning)) {
        STD_string seqname(valstr.c_str());

        // Add sequence variant/options
        STD_string seqopts;
        if(!check_status(func,"ScanningSequence",dset->findAndGetOFString(DCM_ScanningSequence, valstr),warning)) {
          seqopts+=valstr.c_str();
        }
        if(!check_status(func,"SequenceVariant",dset->findAndGetOFStringArray(DCM_SequenceVariant, valstr),warning)) {
          if(seqopts!="" && valstr!="") seqopts+="_";
          seqopts+=replaceStr(valstr.c_str(),"\\","_");
        }
        if(!check_status(func,"ScanOptions",dset->findAndGetOFStringArray(DCM_ScanOptions, valstr),warning)) {
          if(seqopts!="" && valstr!="") seqopts+="_";
          seqopts+=replaceStr(valstr.c_str(),"\\","_");
        }
        if(seqopts!="") seqname+="_"+seqopts;
        prot.seqpars.set_Sequence(seqname);
      }
    }


    // Timestamp temporal sorting of multi-file DICOMs
    prot.seqpars.set_AcquisitionStart(-1.0); // inidicates no useful value
    if(dset->findAndGetOFString(DCM_FrameReferenceTime, valstr).good()){ // for dynamic PET data
      ODINLOG(odinlog,normalDebug) << "Using AcquisitionStart=FrameReferenceTime=" << valstr << STD_endl;
      prot.seqpars.set_AcquisitionStart(atof(valstr.c_str()));
    } else {
      if(dset->findAndGetOFString(DCM_AcquisitionTime, valstr).good()){// Try timestamp field first
        ODINLOG(odinlog,normalDebug) << "Using AcquisitionStart=AcquisitionTime=" << valstr << STD_endl;
        LONGEST_INT sec, sec_series;
        double fract, fract_series;
        timestr2seconds(valstr, sec, fract);
        timestr2seconds(seriestimestr, sec_series, fract_series);
        ODINLOG(odinlog,normalDebug) << "sec / fract = " << sec << " / " << fract << STD_endl;
        ODINLOG(odinlog,normalDebug) << "sec_series / fract_series = " << sec_series << " / " << fract_series << STD_endl;
        prot.seqpars.set_AcquisitionStart( 1000.0*(sec-sec_series) + 1000.0*(fract-fract_series) );  // subtracting series time to avoid overflow of double precision
        ODINLOG(odinlog,normalDebug) << "AcquisitionStart=" << ftos(prot.seqpars.get_AcquisitionStart(),10) << STD_endl;

      } else {
        if(dset->findAndGetOFString(DCM_InstanceNumber, valstr).good()) { // fallback
          ODINLOG(odinlog,normalDebug) << "Using AcquisitionStart=InstanceNumber=" << valstr << STD_endl;
          prot.seqpars.set_AcquisitionStart(atof(valstr.c_str()));
        }
      }
    }

    if(modality=="PT") {
      if(!check_status(func,"ActualFrameDuration",dset->findAndGetOFString(DCM_ActualFrameDuration, valstr),optional)) {
        prot.seqpars.set_RepetitionTime(atof(valstr.c_str()));
        ODINLOG(odinlog,normalDebug) << "Start / TR(PET) = " << prot.seqpars.get_AcquisitionStart() << " / " << prot.seqpars.get_RepetitionTime() << STD_endl;
      }
    }

    return shape(0)*shape(1);
  }


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


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


    Range all=Range::all();

    STD_string delimiter("\\");

    const char* func="DicomFormat::write";
    if(check_dict(func)) return -1;

    LDRfileName fname(filename);

#ifdef OFLOG_H
    OFLog::configure(OFLogger::WARN_LOG_LEVEL); // reset logging (in dcmtk 3.6.0 info logging is somehow enabled sometimes)
#endif

    // Create directory to store DICOMS
    DicomDirInterface dicomdir;
    LDRfileName dstdir; dstdir.set_dir(true);
    dstdir=fname.get_dirname()+SEPARATOR_STR+fname.get_basename_nosuffix()+STD_string("_dcm");
    if( createdir(dstdir.c_str()) && (!checkdir(dstdir.c_str())) ) {
      ODINLOG(odinlog,errorLog) << "Cannot create directory " << dstdir << STD_endl;
      return -1;
    }
    STD_string dcmdirfname=dstdir+SEPARATOR_STR+"DICOMDIR";
    if(check_status(func,"createNewDicomDir",dicomdir.createNewDicomDir(DicomDirInterface::AP_GeneralPurpose,dcmdirfname.c_str()))) return -1;


    // Iterate over protocol-data pairs
    unsigned int index=0;
    for(FileIO::ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit) {

      Protocol prot(pdit->first); // writable copy to adjust FOV and geometry mode
      const Data<float,4>& data=pdit->second;

      if(prot.geometry.get_Mode()==voxel_3d) {
        ODINLOG(odinlog,normalDebug) << "Forcing slicepack mode" << STD_endl;
        float slicespacing=secureDivision(prot.geometry.get_FOV(sliceDirection), data.extent(sliceDim));
        prot.geometry.set_Mode(slicepack);
        prot.geometry.set_nSlices(data.extent(sliceDim));
        prot.geometry.set_sliceDistance(slicespacing);
        prot.geometry.set_sliceThickness(slicespacing);
      }

      DcmFileFormat fileformat;

      // The meta-info stuff not provided by dcmtk
      DcmMetaInfo* metainfo=fileformat.getMetaInfo();
      if(check_status(func,"MediaStorageSOPClassUID", metainfo->putAndInsertString(DCM_MediaStorageSOPClassUID,UID_MRImageStorage))) return -1;

      // The dataset stuff
      DcmDataset *dataset=fileformat.getDataset();

      // The static UIDs
      char uid[100];
      if(check_status(func,"SOPClassUID",         dataset->putAndInsertString(DCM_SOPClassUID,         dcmGenerateUniqueIdentifier(uid, SITE_INSTANCE_UID_ROOT)))) return -1;
      if(check_status(func,"FrameOfReferenceUID", dataset->putAndInsertString(DCM_FrameOfReferenceUID, dcmGenerateUniqueIdentifier(uid, SITE_INSTANCE_UID_ROOT)))) return -1;
      if(check_status(func,"SeriesInstanceUID",   dataset->putAndInsertString(DCM_SeriesInstanceUID,   dcmGenerateUniqueIdentifier(uid, SITE_INSTANCE_UID_ROOT)))) return -1;

      // Common stuff
      if(check_status(func,"ImageType",            dataset->putAndInsertString(DCM_ImageType,("ORIGINAL"+delimiter+"PRIMARY"+delimiter+"OTHER").c_str())))  return -1;
      if(check_status(func,"SpecificCharacterSet", dataset->putAndInsertString(DCM_SpecificCharacterSet,"ISO_IR 100")))  return -1;
      if(check_status(func,"SOPClassUID",          dataset->putAndInsertString(DCM_SOPClassUID,UID_MRImageStorage))) return -1;

      // Pixel data stuff
      if(check_status(func,"SamplesPerPixel",          dataset->putAndInsertUint16(DCM_SamplesPerPixel, 1))) return -1;
      if(check_status(func,"PhotometricInterpretation",dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2"))) return -1;
      if(check_status(func,"PixelRepresentation",      dataset->putAndInsertString(DCM_PixelRepresentation, "0")))  return -1;

      // Sequence type stuff
      if(check_status(func,"ScanningSequence",         dataset->putAndInsertString(DCM_ScanningSequence,    "RM"))) return -1; // research mode
      if(check_status(func,"SequenceVariant",          dataset->putAndInsertString(DCM_SequenceVariant,     "NONE"))) return -1;
      if(check_status(func,"ScanOptions",              dataset->putAndInsertString(DCM_ScanOptions,         ""))) return -1; // bogus value
      if(check_status(func,"MRAcquisitionType",        dataset->putAndInsertString(DCM_MRAcquisitionType,   "2D"))) return -1;
      if(check_status(func,"EchoTrainLength",          dataset->putAndInsertString(DCM_EchoTrainLength,     "1"))) return -1; // bogus value


      // date/time of scan, will be used for all DICOM timestamps
      STD_string scandate;
      STD_string scantime;
      prot.study.get_DateTime(scandate,scantime);


      // Patient stuff
      STD_string pname;
      STD_string pid;
      STD_string pdbirth;
      char psex;
      float pweight, psize;
      prot.study.get_Patient(pid, pname, pdbirth, psex, pweight, psize);
      pdbirth[ODIN_DATE_LENGTH]='\0';
      pname=replaceStr(pname," ","^");
      if(pid=="") pid="Zero"; // dicomdir.addDicomFile does not like zero-size string
      if(check_status(func,"PatientsName",      dataset->putAndInsertString(DCM_PAT_NAME,      pname.c_str()))) return -1;
      if(check_status(func,"PatientID",         dataset->putAndInsertString(DCM_PatientID,     pid.c_str()))) return -1;
      if(check_status(func,"PatientsBirthDate", dataset->putAndInsertString(DCM_PAT_BIRTHDATE, pdbirth.c_str()))) return -1;
      if(check_status(func,"PatientsSex",       dataset->putAndInsertString(DCM_PAT_SEX,       STD_string(1,psex).c_str()))) return -1;
      if(check_status(func,"PatientsWeight",    dataset->putAndInsertString(DCM_PAT_WEIGHT,    ftos(pweight).c_str()))) return -1;
      if(check_status(func,"PatientsSize",      dataset->putAndInsertString(DCM_PAT_SIZE,      ftos(0.001*psize).c_str()))) return -1;

      // Calculate patient age in years
      int ibirth=atoi(pdbirth.c_str());
      int iscan=atoi(scandate.c_str());
      int ageyears=int(1.0e-4*(iscan-ibirth)); // use months/days as digits after comma to round down age
      ODINLOG(odinlog,normalDebug) << "ageyears=" << ageyears << STD_endl;
      if(check_status(func,"PatientAge",    dataset->putAndInsertString(DCM_PAT_AGE,    (itos(ageyears)+"Y").c_str()))) return -1;


      // Study stuff
      STD_string studydescr;
      STD_string physician;
      prot.study.get_Context(studydescr, physician);
      studydescr=replaceStr(studydescr," ","^");
      physician=replaceStr(physician," ","^");

      // Mangle date, PatientID and study description into unique Study UID
      STD_string studyinstanceuid;
      studyinstanceuid+=scandate;
      for(unsigned int i=0; i<pid.length(); i++) studyinstanceuid+="."+itos(pid[i]);
      for(unsigned int i=0; i<studydescr.length(); i++) studyinstanceuid+="."+itos(studydescr[i]);
      if(studyinstanceuid.length()>64) studyinstanceuid=studyinstanceuid.substr(0,64);

      if(check_status(func,"StudyInstanceUID",         dataset->putAndInsertString(DCM_StudyInstanceUID, studyinstanceuid.c_str()))) return -1;
      if(check_status(func,"StudyDate",                dataset->putAndInsertString(DCM_StudyDate,        scandate.c_str()))) return -1;
      if(check_status(func,"StudyTime",                dataset->putAndInsertString(DCM_StudyTime,        scantime.c_str()))) return -1;
      if(check_status(func,"PerformingPhysiciansName", dataset->putAndInsertString(DCM_PERF_PHYS,        physician.c_str()))) return -1;
      if(check_status(func,"ReferringPhysiciansName",  dataset->putAndInsertString(DCM_REF_PHYS,         physician.c_str()))) return -1;
      if(check_status(func,"StudyID",                  dataset->putAndInsertString(DCM_StudyID,          "1"))) return -1; // notyet defined in ODIN
      if(check_status(func,"StudyDescription",         dataset->putAndInsertString(DCM_StudyDescription, studydescr.c_str()))) return -1;
      if(check_status(func,"AccessionNumber",          dataset->putAndInsertString(DCM_AccessionNumber,  ""))) return -1;

      // Series stuff
      STD_string seriesdescr;
      int seriesno;
      prot.study.get_Series(seriesdescr, seriesno);
      seriesdescr=replaceStr(seriesdescr," ","^");
      if(seriesno<=0) seriesno=1;
      if(check_status(func,"Modality",    dataset->putAndInsertString(DCM_Modality,     "MR"))) return -1;
      if(check_status(func,"SeriesDescription",dataset->putAndInsertString(DCM_SeriesDescription, seriesdescr.c_str()))) return -1; // must have value for DICOMDIR
      if(check_status(func,"ProtocolName",dataset->putAndInsertString(DCM_ProtocolName, seriesdescr.c_str()))) return -1;
      if(check_status(func,"SeriesNumber",dataset->putAndInsertString(DCM_SeriesNumber, itos(seriesno).c_str()))) return -1; // must have value for DICOMDIR
      if(check_status(func,"SeriesDate",  dataset->putAndInsertString(DCM_SeriesDate, scandate.c_str()))) return -1;
      if(check_status(func,"SeriesTime",  dataset->putAndInsertString(DCM_SeriesTime, scantime.c_str()))) return -1;

      // Other date/time stuff
      if(check_status(func,"AcquisitionDate",  dataset->putAndInsertString(DCM_AcquisitionDate, scandate.c_str()))) return -1;
      if(check_status(func,"ContentDate",      dataset->putAndInsertString(DCM_ContentDate, scandate.c_str()))) return -1;
      if(check_status(func,"ContentTime",      dataset->putAndInsertString(DCM_ContentTime, scantime.c_str()))) return -1;

      // Frame of Reference stuff
      if(check_status(func,"PositionReferenceIndicator",dataset->putAndInsertString(DCM_PositionReferenceIndicator, ""))) return -1;



      int xsize=data.extent(readDim);
      int ysize=data.extent(phaseDim);
      if(opts.dialect=="siemens") {
        // Siemens screws up with some non-quadratic image sizes, so we will use next-power-of-2-sized quadratic images
        int quadsize = 1 << (int)(ceil( log(double(STD_max(data.extent(phaseDim),data.extent(readDim))))/ log(2.0)));
        ODINLOG(odinlog,normalDebug) << "Fixing matrix size for siemens: " << xsize << "x" << ysize << " -> " << quadsize << "x" << quadsize << STD_endl;
        xsize=ysize=quadsize;
      }
      int xoffset=(xsize-data.extent(readDim))/2;
      int yoffset=(ysize-data.extent(phaseDim))/2;
      ODINLOG(odinlog,normalDebug) << "xsize/ysize/xoffset/yoffset=" << xsize << "/" << ysize << "/" << xoffset << "/" << yoffset << STD_endl;

      // adjust FOV in protocol to new matrix size
      float FOVread= secureDivision(xsize,data.extent(readDim))*prot.geometry.get_FOV(readDirection);
      float FOVphase=secureDivision(ysize,data.extent(phaseDim))*prot.geometry.get_FOV(phaseDirection);
      prot.geometry.set_FOV(readDirection,  FOVread);
      prot.geometry.set_FOV(phaseDirection, FOVphase);

      // Matrix size
      if(check_status(func,"Rows",   dataset->putAndInsertUint16(DCM_Rows,    ysize))) return -1;
      if(check_status(func,"Columns",dataset->putAndInsertUint16(DCM_Columns, xsize))) return -1;

      // seqpars stuff
      if(check_status(func,"RepetitionTime",  dataset->putAndInsertString(DCM_RepetitionTime,   ftos(prot.seqpars.get_RepetitionTime()).c_str()))) return -1;
      if(check_status(func,"TemporalResolution",dataset->putAndInsertString(DCM_TemporalResolution,   ftos(prot.seqpars.get_RepetitionTime()).c_str()))) return -1;
      if(check_status(func,"EchoTime",        dataset->putAndInsertString(DCM_EchoTime,         ftos(prot.seqpars.get_EchoTime()).c_str()))) return -1;
//    if(check_status(func,"NumberOfAverages",dataset->putAndInsertString(DCM_NumberOfAverages, itos(prot.seqpars.get_NumOfRepetitions()).c_str()))) return -1;
      float pixelbw=secureDivision(1000.0*prot.seqpars.get_AcqSweepWidth(),prot.seqpars.get_MatrixSize(readDirection));
      if(check_status(func,"PixelBandwidth",  dataset->putAndInsertString(DCM_PixelBandwidth,   ftos(pixelbw).c_str()))) return -1;
      if(check_status(func,"FlipAngle",       dataset->putAndInsertString(DCM_FlipAngle,        ftos(prot.seqpars.get_FlipAngle()).c_str()))) return -1;
      if(check_status(func,"SequenceName",    dataset->putAndInsertString(DCM_SequenceName,     prot.seqpars.get_Sequence().c_str()))) return -1;

      // system stuff
      if(check_status(func,"TransmitCoilName",     dataset->putAndInsertString(DCM_TransmitCoilName,      prot.system.get_transmit_coil_name().c_str()))) return -1;
      if(check_status(func,"ImagedNucleus",        dataset->putAndInsertString(DCM_ImagedNucleus,         prot.system.get_main_nucleus().c_str()))) return -1;
      if(check_status(func,"MagneticFieldStrength",dataset->putAndInsertString(DCM_MagneticFieldStrength, ftos(0.001*prot.system.get_B0()).c_str()))) return -1;
      if(check_status(func,"ImagingFrequency",     dataset->putAndInsertString(DCM_ImagingFrequency,      ftos(0.001*prot.system.get_nuc_freq()).c_str()))) return -1;
      STD_string swstring=STD_string(PACKAGE)+"-"+VERSION;
      if(check_status(func,"SoftwareVersions",     dataset->putAndInsertString(DCM_SoftwareVersions,      swstring.c_str()))) return -1;

      // Manufacturer stuff
      if(check_status(func,"Manufacturer",dataset->putAndInsertString(DCM_Manufacturer, prot.system.get_platform_str().c_str()))) return -1;
      STD_string fieldstr=ftos(0.001*prot.system.get_B0(),1)+" Tesla"; // Use field strength to label station
      if(check_status(func,"StationName",dataset->putAndInsertString(DCM_StationName, fieldstr.c_str()))) return -1;

      // geometry stuff
      float xspacing=secureDivision(FOVread, xsize);
      float yspacing=secureDivision(FOVphase, ysize);
      STD_string pixstr=ftos(xspacing)+delimiter+ftos(yspacing);
      if(check_status(func,"PixelSpacing",dataset->putAndInsertString(DCM_PixelSpacing, pixstr.c_str()))) return -1;
      if(check_status(func,"NumberOfSlices",dataset->putAndInsertString(DCM_NumberOfSlices, itos(prot.geometry.get_nSlices()).c_str()))) return -1;
      if(check_status(func,"SliceThickness",dataset->putAndInsertString(DCM_SliceThickness, ftos(prot.geometry.get_sliceThickness()).c_str()))) return -1;
      if(check_status(func,"SpacingBetweenSlices",dataset->putAndInsertString(DCM_SpacingBetweenSlices, ftos(prot.geometry.get_sliceDistance()).c_str()))) return -1;


      // DTI
      LDRbase* bptr=prot.get_parameter("Diffusion_bVector");
      if(bptr) {
        LDRtriple* bvec=0;
        bvec=bptr->cast(bvec);
        if(bvec) {
          fvector fbvec(*bvec);
          float bval=norm3(fbvec[0],fbvec[1],fbvec[2]);
          ODINLOG(odinlog,normalDebug) << "bval=" << bval << STD_endl;
          if(check_status(func,"DiffusionBValue",dataset->putAndInsertString(DCM_DiffusionBValue, ftos(bval).c_str()))) return -1;

          if(bval) fbvec/=bval;
          STD_string borient=ftos(fbvec[0])+delimiter+ftos(fbvec[1])+delimiter+ftos(fbvec[2]);
          if(check_status(func,"DiffusionGradientOrientation", dataset->putAndInsertString(DCM_DiffusionGradientOrientation, borient.c_str())))  return -1;

        }
      }

      // for SliceLocation
      dvector slice_offset=prot.geometry.get_sliceOffsetVector();

      // Image orientation
      // Note: x- and y-axis are inverted when converting from ODIN coord system (as defined in geometry.h) to DICOM coord system
      int digits=12; // Need high accuracy, otherwise DICOMs won't import on Siemens...
      dvector linvec=prot.geometry.get_readVector();
      dvector colvec=prot.geometry.get_phaseVector();
      STD_string patort=ftos(-linvec[0],digits)+delimiter+
                        ftos(-linvec[1],digits)+delimiter+
                        ftos( linvec[2],digits)+delimiter+
                        ftos(-colvec[0],digits)+delimiter+
                        ftos(-colvec[1],digits)+delimiter+
                        ftos( colvec[2],digits);
      if(check_status(func,"ImageOrientationPatient", dataset->putAndInsertString(DCM_ImageOrientationPatient, patort.c_str())))  return -1;

      // fixed position in odin: head first supine
      if(check_status(func,"PatientPosition", dataset->putAndInsertString(DCM_PatientPosition, "HFS")))  return -1;

      // use get_cornerPoints() for ImagePositionPatient
      // account for DICOM's convention that the position is center of the 1st voxel:
      prot.geometry.set_FOV(readDirection,  FOVread - xspacing);
      prot.geometry.set_FOV(phaseDirection, FOVphase - yspacing);
      darray cornerpoints=prot.geometry.get_cornerPoints(Geometry(),0);
      ODINLOG(odinlog,normalDebug) << "cornerpoints" << STD_string(cornerpoints.get_extent()) << "=" << cornerpoints.printbody() << STD_endl;

      Data<u16bit,4> u16bitdata;
      data.convert_to(u16bitdata,!opts.noscale);

      if(opts.dialect=="siemens") {
        // clip data for viewer on host
        u16bit thresh=std::numeric_limits<u16bit>::max()-10;
        clip_max(u16bitdata,thresh);
      }

      u16bit min_stored=min(u16bitdata);
      u16bit max_stored=max(u16bitdata);

      unsigned int bits_used=16;
      if(opts.noscale) {
        bits_used=STD_max((unsigned int)1, (unsigned int)ceil(log(max_stored)/log(2)));
      }
      ODINLOG(odinlog,normalDebug) << "bits_used=" << bits_used << STD_endl;
      if(check_status(func,"BitsAllocated", dataset->putAndInsertString(DCM_BitsAllocated, "16"))) return -1;
      if(check_status(func,"BitsStored",    dataset->putAndInsertString(DCM_BitsStored, itos(bits_used).c_str()))) return -1;
      if(check_status(func,"HighBit",       dataset->putAndInsertString(DCM_HighBit, itos(bits_used-1).c_str()))) return -1;

      u16bit windowwidth=max_stored-min_stored;
      u16bit windowcenter=min_stored+windowwidth/2;
      ODINLOG(odinlog,normalDebug) << "windowwidth/windowcenter=" << windowwidth << "/" << windowcenter << STD_endl;
      if(check_status(func,"WindowWidth",  dataset->putAndInsertString(DCM_WindowWidth, itos(windowwidth) .c_str()))) return -1;
      if(check_status(func,"WindowCenter", dataset->putAndInsertString(DCM_WindowCenter, itos(windowcenter).c_str()))) return -1;


      // preserving (floating point) values
      if(!opts.noscale) {
        double minval_out=min(data);
        double maxval_out=max(data);
        double slope=secureDivision(maxval_out-minval_out, max_stored-min_stored);
        double intercept=maxval_out-slope*max_stored;
        if(check_status(func,"RescaleType",dataset->putAndInsertString(DCM_RescaleType, "US"))) return -1;
        if(check_status(func,"RescaleSlope",dataset->putAndInsertString(DCM_RescaleSlope, ftos(slope,6,alwaysExp).c_str()))) return -1;
        if(check_status(func,"RescaleIntercept",dataset->putAndInsertString(DCM_RescaleIntercept, ftos(intercept,6,alwaysExp).c_str()))) return -1;
      }


      Data<u16bit,2> oneslice_padded(ysize,xsize);

      for(int iframe=0; iframe<data.extent(timeDim); iframe++) {
        for(int iplane=0; iplane<data.extent(sliceDim); iplane++) {

          // The instance UID and number, separate for each image
          if(check_status(func,"SOPInstanceUID",dataset->putAndInsertString(DCM_SOPInstanceUID, dcmGenerateUniqueIdentifier(uid, SITE_INSTANCE_UID_ROOT)))) return -1;
          if(check_status(func,"InstanceNumber",dataset->putAndInsertString(DCM_InstanceNumber, itos(index+1).c_str()))) return -1;

          // timestamp
          ODINLOG(odinlog,normalDebug) << "acqstart/TR=" << prot.seqpars.get_AcquisitionStart() << "/" << prot.seqpars.get_RepetitionTime() << STD_endl;
          double timestamp_s=0.001*(prot.seqpars.get_AcquisitionStart()+double(iframe)*prot.seqpars.get_RepetitionTime());
          int hours=int(timestamp_s)/3600;
          timestamp_s-=hours*3600;
          int minutes=int(timestamp_s)/60;
          timestamp_s-=minutes*60;
          int seconds=int(timestamp_s);
          timestamp_s-=seconds;
          int fract=int(timestamp_s*1.0e6);
          char timestr[1000];
          snprintf(timestr, 1000, "%02i%02i%02i.%06i", hours, minutes, seconds, fract);
          if(check_status(func,"AcquisitionTime",dataset->putAndInsertString(DCM_AcquisitionTime, timestr))) return -1;


          // Set slice location of this image
          // Note: x- and y-axis are inverted when converting from ODIN coord system (as defined in geometry.h) to DICOM coord system
          STD_string patpos=ftos(-cornerpoints(iplane,0,0,0,0))+delimiter+ftos(-cornerpoints(iplane,0,0,0,1))+delimiter+ftos(cornerpoints(iplane,0,0,0,2));
          ODINLOG(odinlog,normalDebug) << "patpos(" << iframe << "," << iplane << ")=" << patpos << STD_endl;
          if(check_status(func,"ImagePositionPatient", dataset->putAndInsertString(DCM_ImagePositionPatient, patpos.c_str())))  return -1;

          if(check_status(func,"SliceLocation",dataset->putAndInsertString(DCM_SliceLocation, ftos(slice_offset[iplane]).c_str()))) return -1;

          oneslice_padded=0;
          oneslice_padded(Range(yoffset,yoffset+data.extent(phaseDim)-1), Range(xoffset,xoffset+data.extent(readDim)-1))=u16bitdata(iframe,iplane,all,all);

          if(!insert_uint16_hack(dataset, DCM_SmallestImagePixelValue, min_stored)) return -1;
          if(!insert_uint16_hack(dataset, DCM_LargestImagePixelValue,  max_stored)) return -1;

          if(check_status(func,"PixelData",dataset->putAndInsertUint16Array(DCM_PixelData, oneslice_padded.c_array(), oneslice_padded.numElements()))) return -1;

          STD_string basename(itos(index,10000000)); // Filenames with eight digits
          STD_string singlefname=dstdir+SEPARATOR_STR+basename; // only upper-case letters, no dot, max 8 chars allowed here
          ODINLOG(odinlog,normalDebug) << "singlefname(" << index << ")=" << singlefname << STD_endl;
#if OFFIS_DCMTK_VERSION_NUMBER > 354
          if(check_status(func,"validateMetaInfo",fileformat.validateMetaInfo(EXS_LittleEndianExplicit, EWM_createNewMeta))) return -1; // avoid warnings about differences in SOPInstanceUID
#endif
          if(check_status(func,"saveFile",fileformat.saveFile(singlefname.c_str(), EXS_LittleEndianExplicit))) return -1; // file must be stored before adding to DCOMDIR

          if(check_status(
             func,
             (STD_string("addDicomFile \"")+dstdir+"/"+basename+"\"").c_str(),
             dicomdir.addDicomFile(basename.c_str(),dstdir.c_str())
                         )) return -1;

          index++;
        }
      }

    } // End iterating over protocol-data pairs

    if(check_status(func,"writeDicomDir",dicomdir.writeDicomDir())) return -1;


    return index;
  }

};


#endif




void register_dicom_format() {
#ifdef DICOMSUPPORT
 static DicomFormat df;
 df.register_format();
#endif
}