#!/usr/bin/python

# This file is part of mri_reface
# https://www.nitrc.org/projects/mri_reface/
#
# Copyright 2021-2023 Mayo Foundation for Medical Education and Research.
# This software is accepted by users "as is" and without warranties or
# guarantees of any kind. This software was designed to be used only for
# research purposes, and it is made freely available only for
# non-commercial research use. Contact the authors to obtain information on
# purchasing a separate license for commercial use. Clinical applications
# are not recommended, and this software has NOT been evaluated by the
# United States FDA for any clinical use.
#
# See LICENSE.txt and README.txt for more information.

from __future__ import division, print_function
import sys
import os
import numpy as np
import re
import glob
import copy
import nibabel as nib
import pydicom as dicom
import pydicom._storage_sopclass_uids
from pydicom.sequence import Sequence
from pydicom.dataset import Dataset
from pydicom import DataElement  
try:
    from ADIR.NiiExtensions.NiiExtensionUtility import getDICOMExtension
    have_extension_util = True
except ImportError:
    have_extension_util = False

#experimenting with GNU-zip handling
import gzip

class DCMTag:
    def __init__(self,nVols):
        self.sqTag = ''
        self.sqNum = 0 #checking one PET image with lots of populated sequences showed all sqNum as '0'
        self.sqDescription = ''
        self.tag = ''
        self.description = ''
        self.uniqueValues = [] # across volumes
        self.values = [[] for i in range(nVols)] # per-volume
        self.DataElements = [[] for i in range(nVols)] #per-volume
        self.VM = ''
        self.VR = ''

    def deepKey(self):
        if(self.sqTag):
            return "%s:%d:%s" % (self.sqTag, self.sqNum, self.tag)
        else:
            return "%s" % (self.tag)

    def __repr__(self):
        return "<DCMTag %s:%s     values:%s>" % (self.deepKey(), self.description, self.values)

    def __str__(self):
        return "<DCMTag %s:%s     values:%s>" % (self.deepKey(), self.description, self.values)

def nii2dicom(niifile,outDir,dcmfiles=None,sliceDim=2,preferDicomGeom=True,verbose=False,modality_default='MR',manufacturer_default='Unknown',photometric_default='MONOCHROME2',dry_run=False,uid_root='1.2.826.0.1.3680043.8.498.', dcm_tags=None,seriesNUM_select=None,seriesUID_select=None):
    if(dry_run):
        print('Performing a dry-run, no data will be written.')
    if(not os.path.isdir(outDir) and not dry_run):
        try:
            os.makedirs(outDir)
        #except FileExistsError:  this seems to work for python 3.3 and later
        except OSError as e:  #seems this might work for python 2 https://stackoverflow.com/questions/20790580/python-specifically-handle-file-exists-exception
            if e.errno == errno.EEXIST:
                pass
            else:
                raise

    nii = nib.load(niifile)
    M = nii.header.get_best_affine()
    imgMat3D = nii.dataobj.get_unscaled()
    pixdim = nii.header['pixdim'][1:4]
    nam = os.path.splitext(os.path.basename(niifile))[0]

    # Volumes
    nVols = nii.shape[3] if len(nii.shape)>3 else 1
    if(verbose and nVols>1):
        print("Detected multiple volumes: " + str(nVols))

    # Check if DICOM directory is provide and is a directory
    if(dcmfiles is None):
        input_dicoms = 0
        if(verbose):
            print('No DICOM directory was input. We will produce new DICOM from scratch.')
    elif(not os.path.isdir(dcmfiles)):
        input_dicoms = 0
        print('Error: Input DICOM path is not a directory.')
        return 1
    else:
        if(os.path.isdir(dcmfiles)):
            # store this as _temp because we may need to print location for verbose output.
            dcmfiles_temp = glob.glob(dcmfiles + os.path.sep + '*')
            # Exclude some filetypes that might be cluttering the dcm directory
            dcmfiles_temp = [i for i in dcmfiles_temp if not re.search(r'\.xml$', i) ]
            dcmfiles_temp = [i for i in dcmfiles_temp if not re.search(r'\.json$', i) ]
            input_dicoms = len(dcmfiles_temp)
            if(input_dicoms == 0):    #check if directory is empty
                print('Error: Input DICOM file directory is empty.')
                return 1
            else:
                if(verbose):
                    print('Using ' + str(input_dicoms) + ' dcm file(s) in ' + dcmfiles);
                dcmfiles = dcmfiles_temp   #done with the directory name, saving as dcmfiles for code compatability below
                # Load all dicom, sort by instancenumber
                dcm_by_instance = {}

                #check Instance Number of two slices
                dicom_0 = dicom.dcmread(gzip.open(dcmfiles[0],mode='rb')) if(dcmfiles[0][-3:] == '.gz') else dicom.dcmread(dcmfiles[0])
                if(input_dicoms > 1):
                    dicom_1 = dicom.dcmread(gzip.open(dcmfiles[1],mode='rb')) if(dcmfiles[1][-3:] == '.gz') else dicom.dcmread(dcmfiles[1])
                    concatenated = True if(dicom_0.InstanceNumber == dicom_1.InstanceNumber) else False
                    if(verbose and concatenated):
                        print('Concatenated DICOM image.')
                    if(verbose and not concatenated):
                        print('Non-concatenated DICOM image.')
                else:
                    concatenated = False

                #populate dcm_by_instance
                series_numbers = []
                series_uids    = []
                for dcmfn in dcmfiles:
                    if(dcmfn[-3:] == '.gz'):
                        dicom_i = dicom.dcmread(gzip.open(dcmfn,mode='rb'))
                        if(verbose): print('DICOM file is GNU-zipped, handling of this file is experimental.')
                    else:
                        dicom_i = dicom.dcmread(dcmfn) # We could use stop_before_pixels=True but it is measurably SLOWER. Others have reported the same.
                    if((seriesNUM_select is not None and int(dicom_i.SeriesNumber) != int(seriesNUM_select)) or (seriesUID_select is not None and str(dicom_i.SeriesInstanceUID) != str(seriesUID_select))):
                        continue

                    if(not concatenated):
                        if(int(dicom_i.InstanceNumber) in dcm_by_instance.keys()):
                            print("ERROR: Duplicate Instance Number detected (non-concatenated).")
                            if(dcm_by_instance[dicom_i.InstanceNumber].SeriesNumber != dicom_i.SeriesNumber):
                                print('ERROR: Duplicate Instance Numbers have different Series Numbers, Series Numbers {} and {} detected, other Series Numbers may also be present.'.format(dcm_by_instance[dicom_i.InstanceNumber].SeriesNumber,dicom_i.SeriesNumber))
                            if(dcm_by_instance[dicom_i.InstanceNumber].SeriesInstanceUID != dicom_i.SeriesInstanceUID):
                                print('ERROR: Duplicate Instance Numbers have different Series Instance UID, Series Instance UIDs {} and {} detected, other Series Instance UIDs may also be present.'.format(dcm_by_instance[dicom_i.InstanceNumber].SeriesInstanceUID,dicom_i.SeriesInstanceUID))

                            return 1
                        else:
                            dcm_by_instance[int(dicom_i.InstanceNumber)] = dicom_i
                    else:
                        try:
                            if(int(dicom_i.InConcatenationNumber) in dcm_by_instance.keys()):
                                print("ERROR: Duplicate In-Concatenation Number detected (concatenated).")
                                if(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesNumber != dicom_i.SeriesNumber):
                                    print('ERROR: Duplicate In-Concatenation Numbers have different Series Numbers, Series Numbers {} and {} detected, other Series Number may also be present.'.format(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesNumber,dicom_i.SeriesNumber))
                                if(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesInstanceUID != dicom_i.SeriesInstanceUID):
                                    print('ERROR: Duplicate In-Concatenation Numbers have different Series Instance UIDs, Series Instance UIDs {} and {} detected, other Series Instance UIDs may also be present.'.format(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesInstanceUID,dicom_i.SeriesInstanceUID))
                                raise Exception("ERROR: Duplicate In-Concatenation Number detected (concatenated).")
                            else:
                                dcm_by_instance[int(dicom_i.InConcatenationNumber)] = dicom_i
                        except:
                            if((int(dicom_i.ConcatenationFrameOffsetNumber) + 1) in dcm_by_instance.keys()):
                                print("ERROR: Duplicate Concatenation Frame Offset Number detected (concatenated).")
                                if(dcm_by_instance[int(dicom_i.ConcatenationFrameOffsetNumber) + 1].SeriesNumber != dicom_i.SeriesNumber):
                                    print('ERROR: Duplicate Concatenation Frame Offset Numbers have different Series Numbers, Series Numbers {} and {} detected, other Series Numbers may also be present.'.format(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesNumber,dicom_i.SeriesNumber))
                                if(dcm_by_instance[int(dicom_i.ConcatenationFrameOffsetNumber) + 1].SeriesInstanceUID != dicom_i.SeriesInstanceUID):
                                    print('ERROR: Duplicate Concatenation Frame Offset Numbers have different Series Instance UIDs, Series Instance UIDs {} and {} detected, other Series Instance UIDs may also be present.'.format(dcm_by_instance[dicom_i.InConcatenationNumber].SeriesInstanceUID,dicom_i.SeriesInstanceUID))
                                raise Exception("ERROR: Duplicate Concatenation Frame Offset Number detected (concatenated).")
                            else:
                                dcm_by_instance[int(dicom_i.ConcatenationFrameOffsetNumber) + 1] = dicom_i
                    if(dicom_i.SeriesNumber not in series_numbers):series_numbers.insert(0,dicom_i.SeriesNumber)
                    if(dicom_i.SeriesInstanceUID not in series_uids):series_uids.insert(0,dicom_i.SeriesInstanceUID)
                if(len(series_numbers) > 1):
                    print('ERROR: Multiple Series Numbers detected, Series Numbers {} are present in the original DICOM directory.'.format(series_numbers))
                    return 1
                if(len(series_uids) > 1):
                    print('ERROR: Multiple Series Instance UIDs detected, Series Instance UIDs {} are present in the original DICOM directory.'.format(series_uids))
                    return 1
                if(len(series_numbers) == 0 or len(series_uids) == 0):
                    print('ERROR: All DICOM files have been filtered by seriesNUM_select/seriesUID_select.')
                    return 1
                if(len(dcm_by_instance.keys()) == sorted(dcm_by_instance.keys())[-1] and (len(dcm_by_instance.keys()) == len(dcmfiles) or seriesNUM_select is not None or seriesUID_select is not None)):
                    if(verbose): print("All DICOM slices appear to be present")
                else:
                     print("ERROR: Missing DICOM slice(s) detected, {} files present in DICOM directory, {} unique Instance/Concatenation Numbers present, and {} as the highest Instance/Concatenation Number.".format(len(dcmfiles),len(dcm_by_instance.keys()),dcm_by_instance.keys()[-1]))
                     raise Exception("ERROR: Missing DICOM slice(s) detected, {} files present in DICOM directory, {} unique Instance/Concatenation Numbers present, and {} as the highest Instance/Concatenation Number.".format(len(dcmfiles),len(dcm_by_instance.keys()),dcm_by_instance.keys()[-1]))
                if("Media Storage SOP Class UID" in str(dcm_by_instance[1].file_meta) and 'enhanced' in dcm_by_instance[1].file_meta.MediaStorageSOPClassUID.name.lower()):
                    print('ERROR: Enhanced DICOM provided; This is not supported.')
                    print('Use dcuncat with -unenhance flag to produce individual DICOM slices. If this is not possible, rerun without providing DICOM directory.')
                    return 1

                is_mosaic = False
                mosaic_images = 0
                mosaic_ref_times = []
                if(nVols>1):
                    # Determine if image is MOSAIC
                    is_mosaic = 'MOSAIC' in dcm_by_instance[1].ImageType
                    if('(0019, 100a)' in str(dcm_by_instance[1]) or '(0019,100a)' in str(dcm_by_instance[1])):
                        is_mosaic = True
                        if(type(dcm_by_instance[1][0x19100a].value) == int):
                            mosaic_images = dcm_by_instance[1][0x19100a].value
                    if('(0019, 1029)' in str(dcm_by_instance[1]) or '(0019,1029)' in str(dcm_by_instance[1])):
                        is_mosaic = True
                        if(type(dcm_by_instance[1][0x191029].value) == list):
                            mosaic_ref_times = dcm_by_instance[1][0x191029].value
                    if(('(0051, 1016)' in str(dcm_by_instance[1]) or '(0051,1016)' in str(dcm_by_instance[1])) and ('mosaic' in str.lower(str(dcm_by_instance[1][0x511016].value)))):
                        is_mosaic = True
                    if(verbose and is_mosaic):
                       print('Image is MOSAIC.')

                    # Determine 'coordinate order'
                    slice_diff = np.array(dcm_by_instance[1].ImagePositionPatient) - np.array(dcm_by_instance[2].ImagePositionPatient)
                    if(np.count_nonzero(slice_diff)>0 or dcm_by_instance[1].Manufacturer == 'SIEMENS'): #dcm2niix seems to assume Siemens uses xyzt, the manufacturer test is added from compatability
                        slice_order_4d = 'xyzt'
                    else:
                        slice_order_4d = 'xytz'

                # Create a list of unique values for each tag
                if dcm_tags is None:
                    dcm_tags = {}
                v = 0

                tag_blocklist = [ '(0008, 0018)',  # SOPInstanceUID (file-unique and we will overwrite)
                                  '(0020, 0013)',  # InstanceNumber (file-unique and we will overwrite)
                                  '(0020, 1041)',  # SliceLocation (vendors have inconsistent meanings for this)
                                  '(0028, 0107)',  # LargestImagePixelValue (we will overwrite)
                                  '(0028, 1050)',  # WindowCenter
                                  '(0028, 1051)',  # WindowWidth
                                  '(0088, 0200)',  # IconImageSequence (we are not generating icons)
                                  '(2005, 100e)',  # PhilipsScaleSlope (we're using our own scale/slope from nii)
                                  '(0400, 0561)',  # Original Attribute Sequence, throws rare error in pydicom related to (0000, 0400), which appears to be retired tag
                                  '(7fe0, 0010)',  # PixelData

                                  # Mosaic tags, these seems to be mosaic specific
                                  '(0019, 100a)',  # NumberOfImagesInMosaic,     int
                                  '(0019, 1029)' ] # MosaicRefAcqTime,           list of floats

                if(is_mosaic):
                    tag_blocklist.extend(['(0029, 1010)',  # CSA "image" header includes MOSAIC status; we can only keep or remove the tag, so we must remove it.
                                          '(0029, 0010)',  # CSA "image" header includes MOSAIC status; we can only keep or remove the tag, so we must remove it.
                                          #'(0029, 1020)',  # CSA header "series" header is fine to leave intact.
                                          #'(0029, 0020)',  # CSA header "series" header is fine to leave intact.
                                          '(0051, 1016)']) # "private data tag" also marks MOSAIC for a MOSAIC

                for i, dicom_i in sorted(dcm_by_instance.items()):
                    if(nVols>1):
                        if(slice_order_4d == 'xyzt'):
                            if(i>1 and dicom_i.ImagePositionPatient == dcm_by_instance[1].ImagePositionPatient):
                                v = v + 1
                        else: # xytz
                            v = (i-1) % nVols
                    for elem in dicom_i:
                        # if not in blocklist, and not a group length
                        if(str(elem.tag) not in tag_blocklist and elem.tag.element != 0):
                            if(str(elem.tag) == '(0008, 0008)' and is_mosaic):
                                if('MOSAIC' in elem.value):
                                    elem.value.remove('MOSAIC')
                            newTag = DCMTag(nVols)
                            newTag.tag = elem.tag
                            newTag.VM = elem.VM
                            newTag.VR = elem.VR
                            try:
                                newTag.description = elem.keyword if elem.keyword else elem.description()
                            except:
                                pass
                            if(newTag.deepKey() not in dcm_tags):
                                if(i==1): # If we first encounter the key later than instance 1, it was not consistent
                                    dcm_tags[newTag.deepKey()] = newTag
                                elif(verbose):
                                    print("Found new tag %s:%s in instance %d; It will not be copied." % (newTag.deepKey(), newTag.description, i))
                            if(newTag.deepKey() in dcm_tags):
                                if(elem.value not in dcm_tags[newTag.deepKey()].values[v]):
                                    dcm_tags[newTag.deepKey()].values[v].append(elem.value)
                                    dcm_tags[newTag.deepKey()].DataElements[v].append(elem)
                                if(elem.value not in dcm_tags[newTag.deepKey()].uniqueValues):
                                    dcm_tags[newTag.deepKey()].uniqueValues.append(elem.value)

                            # Sequences get added both as a sequence (above) and as individual items (below)
                            # If the whole sequence never changes across dicom, we can copy the whole sequence at once
                            if(elem.VR == "SQ"):
                                for seq_i in range(0,len(elem.value)):
                                    for elem2 in elem[seq_i]:
                                        if(elem2.tag.element == 0): continue # skip group length tags
                                        newTag = DCMTag(nVols)
                                        newTag.sqTag = elem.tag
                                        newTag.sqNum = seq_i
                                        try:
                                            newTag.sqDescription = elem.keyword if elem.keyword else elem.description()
                                        except:
                                            pass
                                        newTag.tag = elem2.tag
                                        try:
                                            newTag.description = elem2.keyword if elem2.keyword else elem2.description()
                                        except:
                                            pass
                                        if(newTag.deepKey() not in dcm_tags):
                                            dcm_tags[newTag.deepKey()] = newTag
                                        if(elem2.value not in dcm_tags[newTag.deepKey()].values[v]):
                                            dcm_tags[newTag.deepKey()].values[v].append(elem2.value)
                                            dcm_tags[newTag.deepKey()].DataElements[v].append(elem2)
                                        if(elem2.value not in dcm_tags[newTag.deepKey()].uniqueValues):
                                            dcm_tags[newTag.deepKey()].uniqueValues.append(elem2.value)


    # Figure out modality
    if(input_dicoms == 0):
        if(have_extension_util):
            dicomExt = getDICOMExtension(niifile)
        else:
            dicomExt = None

        if(dicomExt is not None):
            modality = dicomExt.get('Modality')
            manufacturer = dicomExt.get('Manufacturer')
        else:
            modality = modality_default
            manufacturer = manufacturer_default
    else:
        manufacturer = dcm_by_instance[1].Manufacturer
        modality = dcm_by_instance[1].Modality

    if(verbose):
        print("Guessed modality: " + modality)
        print("Guessed manufacturer: " + manufacturer)

    if(modality == 'MR'):
        SOP_UID = pydicom._storage_sopclass_uids.MRImageStorage
    elif(modality == 'PT'):
        SOP_UID = pydicom._storage_sopclass_uids.PositronEmissionTomographyImageStorage
    elif(modality == 'CT'):
        SOP_UID = pydicom._storage_sopclass_uids.CTImageStorage
    elif(modality == 'DX'):
        preferDicomGeom = 0
        SOP_UID = pydicom._storage_sopclass_uids.DigitalXRayImageStorageForPresentation

    # Setup output dicom object
    meta = pydicom.Dataset()
    meta.MediaStorageSOPClassUID = SOP_UID
    meta.TransferSyntaxUID = dicom.uid.ExplicitVRLittleEndian

    staticMetaOut = dicom.dataset.FileDataset('out.dcm',{}, file_meta = meta)
    staticMetaOut.is_little_endian = True
    staticMetaOut.is_implicit_VR = False # If true, then private tags all get written as UN
    staticMetaOut.SOPClassUID = SOP_UID
    staticMetaOut.Modality = modality
    staticMetaOut.Manufacturer = manufacturer
    # Copy what we can from old tags
    if(input_dicoms or dcm_tags is not None):
        if(input_dicoms > 1):
            # Determine the dicom slice dimension. IOP vectors a and b are orthogonal and define the rows and columns of the slices.
            # Slices are "stacked" in a direction normal to the slice, we call this the slice dimension, and is the direction without a 1 in the IOP
            iop_a = abs(np.array(dcm_tags['(0020, 0037)'].uniqueValues[0][0:3]))
            iop_b = abs(np.array(dcm_tags['(0020, 0037)'].uniqueValues[0][3:6]))
            dcm_slice_dim = np.argmin(iop_a + iop_b)

        # First, copy all tags that never changed, AND any sequences (to set up the sequences and their members)
        for key,tag in dcm_tags.items():
            if(not tag.sqTag and (len(tag.uniqueValues)==1 or tag.VR == 'SQ')):
                staticMetaOut[tag.tag] = copy.copy(tag.DataElements[0][0])

        # Now, get rid of any tags within sequences that changed > once
        for key,tag in dcm_tags.items():
            try: #if tag.sqTag is not in metaOut, we get an exception.  This is needed for ADIR_PET_Convert
                if(len(tag.uniqueValues)>1 and tag.sqTag and tag.tag in staticMetaOut[tag.sqTag][tag.sqNum]):
                    del staticMetaOut[tag.sqTag][tag.sqNum][tag.tag]
            except:
                pass
    else:
        staticMetaOut.PhotometricInterpretation = photometric_default #'MONOCHROME2'
        staticMetaOut.SamplesPerPixel = 1
        staticMetaOut.ProtocolName = nam
        staticMetaOut.SeriesDescription = nam
        staticMetaOut.SeriesInstanceUID = pydicom.uid.generate_uid(prefix=uid_root)

    staticMetaOut.BitsAllocated = 16
    staticMetaOut.BitsStored = 16
    staticMetaOut.HighBit = 15
    staticMetaOut.PixelRepresentation = 1

    # handle slope/intercept
    slope = nii.dataobj.slope
    inter = nii.dataobj.inter
    if(not np.isfinite(slope)):
        slope = 1
    if(not np.isfinite(inter)):
        inter = 0

    # Cast float to int, if necessary, while preserving as much as we can and warning when we can't
    do_inf = 0
    if(np.issubdtype(nii.get_data_dtype(),np.floating)):
        if(not np.all(np.isfinite(imgMat3D))): # nan to 0, inf to large
            imgMat3D = imgMat3D.copy()
            print("Warning: found non-finite (nan/inf) values: they will be converted to integers")
            imgMat3D[np.isnan(imgMat3D)] = 0
            ids_neginf = np.isneginf(imgMat3D)
            ids_posinf = np.isposinf(imgMat3D)
            do_inf = 1

        if((slope==0 or slope==1) and inter==0): # if there's already a slope and inter, leave it alone
            imgMat3D_scaled = imgMat3D[:] * slope + inter

            # We used to borrow nibabel's functions to calculate a new optimized scale and intercept, as below.
            # That uses the full int16 range, but this can cause some issues downstream, notable with reg_resample
            #`aw = nib.arraywriters.SlopeInterArrayWriter(imgMat3D_scaled, np.int16)
            #`slope = float(aw.slope)
            #`inter = float(aw.inter)

            # Instead, do the same, but don't use the very top/bottom of the range
            # Inspired by _range_scale in niibabel arraywriters.py
            reserve_range = 500 # reserve x values on each end of the range
            in_min = np.amin(imgMat3D_scaled)
            in_max = np.amax(imgMat3D_scaled)
            in_range = in_max - in_min
            out_min, out_max = np.iinfo(np.int16).min, np.iinfo(np.int16).max # change this for different HighBit
            out_min = out_min + reserve_range
            out_max = out_max - reserve_range
            out_range = out_max - out_min
            slope = in_range / out_range
            inter = in_min - out_min * slope

            if(verbose):
                print("Attempting to improve float image to integer conversion with new slope: %f, inter: %f" % (slope, inter))
            imgMat3D = (imgMat3D_scaled[:] - inter)/slope

        if(not np.array_equal(imgMat3D.astype('int16'),imgMat3D)):
            print('Warning: input is float and must be downcast for dicom. For this image, this was not lossless.')
    imgMat3D = imgMat3D.astype('int16')

    if(do_inf):
        imgMat3D[ids_neginf] = np.iinfo('int16').min
        imgMat3D[ids_posinf] = np.iinfo('int16').max

    def trim_DS(ds,l):
        if('e' in str(ds)):
            (m,e) = str(ds).split('e') #split the mantissa and exponent
            (i,f) = m.split('.') # split the mantissa into integer and fractional component
            len_i = len(i) # save length of integet component
            len_e = len(e) # save length of exponent
            m = round(float(m),l-len_i-len_e-2) # shorten the mantissa to specified length minus length of integer componment and exponent minus two (for the period and the 'e')
            return 'e'.join([str(m),e]) #str() will cut this down to 13 characters, decimal string compatability we will have the length of the mantissa be 12 in almost all cases with an exponent. str() doesn't allow overriding
        else:
            m = "{:.100f}".format(ds) # str(ds) loses a bunch of digits, this keeps extras that will be rounded off
            (i,f) = m.split('.') # split the mantissa into integer and fractional component
            len_i = len(i) # save length of integet component
            m = round(float(m),l-len_i-1) # shorten the mantissa to specified length minus length of integer componment minus one (for the period)
            return m

    #Rescale Slope has DS(Decimal String) VR which specifies 16 byte maximum, pydicom doesn't enforce this
    if(slope != 1): 
        staticMetaOut.RescaleSlope = trim_DS(slope,16)
    if(inter != 0):
        staticMetaOut.RescaleIntercept = trim_DS(inter,16)

    # inspired by spm_imatrix
    R = M[0:3,0:3]
    scale_mat = np.linalg.cholesky(np.matmul(np.transpose(R),R)).copy()
    if(np.linalg.det(R) < 0): scale_mat[0,0] = -scale_mat[0,0]
    flip_mat = np.sign(scale_mat)
    flips = np.diag(flip_mat)
    scale_mat = np.diag(np.diag(scale_mat))    
    rot_mat = np.matmul(R,np.linalg.inv(scale_mat))

    nii_dims = [0,0,0]
    nii_dims[0] = np.argmax(abs(rot_mat[0,:]))
    nii_dims[1] = np.argmax(abs(rot_mat[1,:]))
    nii_dims[2] = np.argmax(abs(rot_mat[2,:]))

    if(input_dicoms > 1):
        slice_dim = nii_dims[dcm_slice_dim]

        # Determine forward or reversed slice order
        # Get nii coordinates corresponding to ipp at instance 1
        ipp_1 = np.array(dcm_by_instance[1].ImagePositionPatient)
        # These are the same sign flips down in save_dcm
        ipp_1[0] = -ipp_1[0]
        ipp_1[1] = -ipp_1[1]
        nii_coords_1 = np.round(np.matmul(np.linalg.inv(M),np.append(ipp_1,1)))
        # nii_coords_1[slice_dim] should be roughly either 0 or nii.shape[slice_dim].
        # If it's closer to 0, reversed=1. If closer to nii.shape, reversed=0
        reversed_slices = 0 if np.argmin(abs(np.array([0,nii.shape[slice_dim]]) - nii_coords_1[slice_dim]))==1 else 1
        # Flip it one more time, if the nifti is flipped in this dimension
        if(flips[slice_dim] == -1): reversed_slices = 1 - reversed_slices
    else:
        slice_dim = 2
        reversed_slices = 1

    # if user wanted to override slice_dim, apply the override
    if(sliceDim != -1): slice_dim = sliceDim

    if(nVols>1 and not input_dicoms):
        # dcm2niix assumes different order per-vendor.
        # "Philips Medical Systems", "GE MEDICAL SYSTEMS", "UNKNOWN" all assume xyzt
        slice_order_4d = "xytz"
        # SIEMENS assumes xyzt, and TemporalPosition does not override that assumption
        if(manufacturer == "SIEMENS"): slice_order_4d = "xyzt"

    if(verbose):
        print('Dimension Mapping: world [0, 1, 2] = nifti %s' % (nii_dims))
        print('Dimension flips: %s' % (flips))
        print('Using slice dimension: ' + str(slice_dim))
        if(reversed_slices):
            print("Using decreasing slice order")
        else:
            print("Using increasing slice order")
        if(nVols>1): print("Using multi-volume slice order: %s" % (slice_order_4d))

    if(slice_dim == 0): #generate IOP from the rotation matrix
        IOP = [rot_mat[0][1],rot_mat[1][1],-rot_mat[2][1],rot_mat[0][2],rot_mat[1][2],-rot_mat[2][2]]
    elif(slice_dim == 1):
        IOP = [rot_mat[0][0],rot_mat[1][0],-rot_mat[2][0],rot_mat[0][2],rot_mat[1][2],-rot_mat[2][2]]
    else: #slice_dim == 2
        IOP = [rot_mat[0][0],rot_mat[1][0],-rot_mat[2][0],rot_mat[0][1],rot_mat[1][1],-rot_mat[2][1]]

    for i in range(6): # Must write values as strings if we want to write them exactly
        IOP[i] = str(IOP[i])

    # Implement preferDicomGeom
    tol = .001
    def replaceIfClose(numIn, dicomValues, tol, desc=''):
        diff = abs(float(numIn) - np.array(dicomValues).astype(np.float))
        minDiff = np.amin(diff)
        minPos = np.argmin(diff)
        if(minDiff < tol):
            if(verbose): print("preferDicomGeom: replacing %s %.15f with \"%s\" from DICOM" % (desc,float(numIn),dicomValues[minPos]))
            return(dicomValues[minPos])
        else:
            return(numIn)

    if(slice_dim == 0):
        staticMetaOut.PixelSpacing = [ '%f' % pixdim[2],'%f' % pixdim[1] ]
        staticMetaOut.SliceThickness = '%f' % pixdim[0]
    elif(slice_dim == 1):
        staticMetaOut.PixelSpacing = [ '%f' % pixdim[2],'%f' % pixdim[0] ]
        staticMetaOut.SliceThickness = '%f' % pixdim[1]
    else: #slice_dim == 2 and catch any invalid value passed by user
        staticMetaOut.PixelSpacing = [ '%f' % pixdim[1], '%f' % pixdim[0] ]
        staticMetaOut.SliceThickness = '%f' % pixdim[2]

    if(input_dicoms and preferDicomGeom):
        for i in range(6):
            IOP[i] = replaceIfClose(IOP[i], [dcm_by_instance[1].ImageOrientationPatient[i].original_string],tol,'IOP[%d]' % (i))

        staticMetaOut.PixelSpacing[0] = replaceIfClose(staticMetaOut.PixelSpacing[0],[dcm_by_instance[1].PixelSpacing[0].original_string],tol,'PixelSpacing[0]')
        staticMetaOut.PixelSpacing[1] = replaceIfClose(staticMetaOut.PixelSpacing[1],[dcm_by_instance[1].PixelSpacing[1].original_string],tol,'PixelSpacing[1]')
        staticMetaOut.SliceThickness = replaceIfClose(staticMetaOut.SliceThickness,[dcm_by_instance[1].SliceThickness.original_string],tol,'SliceThickness')

    staticMetaOut.ImageOrientationPatient = [IOP[0],IOP[1],IOP[2],IOP[3],IOP[4],IOP[5]]
    
    # Prepare for main loop
    nSlices = nii.shape[slice_dim];
    instanceNumber = 0
    #Generate required dicom tags, these should exist for nifti but ADIR_PET_Convert needs this
    if(staticMetaOut.get_item('SeriesInstanceUID') is None):
        staticMetaOut.SeriesInstanceUID = pydicom.uid.generate_uid()
        print('WARNING: Series Instance UID not found, generating new UID.')
    if(staticMetaOut.get_item('StudyInstanceUID') is None):
        staticMetaOut.StudyInstanceUID = pydicom.uid.generate_uid()
        print('WARNING: Study Instance UID not found, generating new UID.')
    if(staticMetaOut.get_item('FrameOfReferenceUID') is None):
        staticMetaOut.FrameOfReferenceUID = pydicom.uid.generate_uid()
        print('WARNING: Frame Instance UID not found, generating new UID.')
    if(staticMetaOut.get_item('NumberOfSlices') is None):
        staticMetaOut.NumberOfSlices = nVols * nSlices
    if(staticMetaOut.get_item('NumberOfTimeSlices') is None):
        staticMetaOut.NumberOfTimeSlices = nVols

    # We had to put this in a function so we can run it in xytz or xyzt loop orders
    def save_dcm():
        metaOut = copy.deepcopy(staticMetaOut)
        metaOut.InstanceNumber = instanceNumber
        metaOut.file_meta.MediaStorageSOPInstanceUID = pydicom.uid.generate_uid(prefix=uid_root)
        metaOut.SOPInstanceUID = metaOut.file_meta.MediaStorageSOPInstanceUID
        if(metaOut.get_item('ImageIndex') is None):
            metaOut.ImageIndex = (2*v + 1) * nSlices + 1 - instanceNumber
        # Assign ContentTime as the corresponding MosaicRefAcqTime if applicable
        if(input_dicoms and is_mosaic and mosaic_images > 0 and mosaic_images == len(mosaic_ref_times)):
            # This needs to be cast as a string or else it throws an error about microseconds.
            if(is_mosaic):
                # Available python libraries for handling time are not working in my environment; it is easy enough to do it manually
                # Grab the volume acquisition time from the correct volume. It appears to be the time at which the volume began being recorded. The format is HHMMSS.
                volume_acq_time = float(dcm_by_instance[((instanceNumber - 1) // nVols) + 1].AcquisitionTime)
                vol_sec = 3600 * (volume_acq_time // 10000) + 60 * ((volume_acq_time % 10000) // 100) + (volume_acq_time % 100) # Convert to seconds

                # Grab the mosaic ref acquistion time. This appears to be the offset from the volume acq time at which the slice was recorded. This appears to be stored in milliseconds. Divide by 1000 to get seconds.
                ref_sec = mosaic_ref_times[(instanceNumber - 1) % mosaic_images] / 1000
                content_sec = vol_sec + ref_sec # Add volume time and mosaic ref time

                content_hour = str(int(content_sec // 3600))
                if(len(content_hour) == 1): content_hour = '0' + content_hour # Pad for formatting

                content_min = str(int((content_sec % 3600) // 60))
                if(len(content_min) == 1): content_min = '0' + content_min

                content_sec = content_sec % 60
                if(content_sec < 10): content_sec = '0' + str(content_sec)
                else: content_sec = str(content_sec)

                metaOut.ContentTime = content_hour + content_min + content_sec
                metaOut.AcquisitionTime = dcm_by_instance[((instanceNumber - 1) // nVols) + 1].AcquisitionTime
        dcmname = outDir + os.path.sep + nam + '_' + '%04d'% (instanceNumber) + '.dcm'

        staticI = nii.shape[0]-1 if(flips[0]==1) else 0
        staticJ = nii.shape[1]-1 if(flips[1]==1) else 0
        staticK = nii.shape[2]-1 if(flips[2]==1) else 0
        if(slice_dim == 0):
            slice = imgMat3D[-i-1,:,:,v] if len(nii.shape)>3 else imgMat3D[-i-1,:,:]
            slice = np.rot90(slice)
            slice = np.fliplr(slice)
            coords = np.matmul(M,np.transpose([ nii.shape[0]-i-1 , staticJ, staticK,1])).transpose()
        elif(slice_dim == 1):
            slice = imgMat3D[:,-i-1,:,v] if len(nii.shape)>3 else imgMat3D[:,-i-1,:]
            slice = np.rot90(slice)
            if(flips[0] == 1): slice = np.fliplr(slice)
            coords = np.matmul(M,np.transpose([ staticI, nii.shape[1]-i-1, staticK,1])).transpose()
        else: #slice_dim == 2 and catch any invalid value passed by user
            slice = imgMat3D[:,:,-i-1,v] if len(nii.shape)>3 else imgMat3D[:,:,-i-1]
            slice = np.rot90(slice)
            coords = np.matmul(M,np.transpose([ staticI, staticJ, nii.shape[2]-i-1,1])).transpose()
            if(flips[0] == 1): slice = np.fliplr(slice)

        # Copy tags specific to this volume that only changed once per volume
        if(nVols > 1 and (dcm_tags is not None or input_dicoms)):
            for key,tag in dcm_tags.items():
                if(len(tag.values[v])==1 and tag.tag not in staticMetaOut and (tag.VR != 'SQ' or tag.description == 'RadionuclideCodeSequence') ):
                    if(not tag.sqTag): #tag is not in sequence
                        metaOut[tag.tag] = tag.DataElements[v][0]
                    else: #tag is in sequence
                        try: #if we try access metaOut[tag.sqTag] for a sequence not yet in metaOut it will throw an exception, if the sequence is in metaOut we can access and write to it by normal means.
                            #when using nii2dicom for converting a niifi to dicom with original dicoms present the sequence will already be in metaOut, we need this special handling when nii2dicom is
                            #called by another program, such as ADIR_PET_Convert where we are making the dicom tags from ECAT/Interfile info.
                            metaOut[tag.sqTag][tag.sqNum][tag.tag] = tag.DataElements[v][0]
                        except:
                            #pyDicom is picky when writing certain tags/sequences.  Using hex values in the 0x00540013 format instead of (0054, 0013) seems to override this and allow writing.
                            #Converting the tag to hex results in a trailing 'L', i.e. hex('(0054, 0013)') becomes '0x540013L', pyDicom doesn't like this so we rstrip the trailing 'L'
                            metaOut[hex(tag.sqTag).rstrip('L')] = DataElement(hex(tag.sqTag).rstrip('L'),'SQ',Sequence([Dataset()]))  
                            metaOut[tag.sqTag][tag.sqNum][tag.tag] = tag.DataElements[v][0] #after the tag/sequence exists in metaOut we can address using the normal tag format.

        metaOut.Rows = slice.shape[0]
        metaOut.Columns = slice.shape[1]
        metaOut.ImagePositionPatient = [ '%f' % -coords[0],'%f' % -coords[1],'%f' % coords[2] ]

        # The short integer format ranges from -32,767 to 32,767 for the signed version and from 0 to 65,535 for the unsigned.
        # PyDicom states signed OR unsigned short for Smallest/LargestImagePixelValue but seems to treat them as needing
        # values to be valid for signed AND unsigned.  We manually enter the tags, allowing us to force US or SS with US prefered.
        # We are using the slope and intercept from the nifti file for all slices, which causes some of the largest values to be out
        # of range for Unsigned. Since these tags are optional, we exclude if value found is not valid.
        smallest_pixel = int(np.amin(slice*slope + inter))
        largest_pixel = int(np.amax(slice*slope + inter))
        if(smallest_pixel != largest_pixel):
            if(smallest_pixel < 0 and smallest_pixel > -32767):   # Signed Short needed
                metaOut[0x0028, 0x0106] = dicom.DataElement(0x00280106, 'SS', smallest_pixel)
            elif(smallest_pixel >= 0 and smallest_pixel < 65535): # Unsigned can be used
                metaOut[0x0028, 0x0106] = dicom.DataElement(0x00280106, 'US', smallest_pixel)
            if(largest_pixel < 0 and largest_pixel > -32767):     # Signed Short needed
                metaOut[0x0028, 0x0107] = dicom.DataElement(0x00280107, 'SS', largest_pixel)
            elif(largest_pixel >= 0 and largest_pixel < 65535):   # Unsigned can be used
                metaOut[0x0028, 0x0107] = dicom.DataElement(0x00280107, 'US', largest_pixel)

        if(preferDicomGeom and input_dicoms): 
            for j in range(3):
                # Tolerance is intentionally looser here because small deviations get compounded over space.
                metaOut.ImagePositionPatient[j] = replaceIfClose(metaOut.ImagePositionPatient[j].original_string,[dcm_by_instance[(instanceNumber - 1) % nSlices + 1].ImagePositionPatient[j].original_string],0.1,'ImagePositionPatient[%d]' % j)

            if(dcm_tags is not None and input_dicoms == 0 and metaOut.Manufacturer == 'SIEMENS'): #only run when we call from outside (ADIR_PET_Convert) and is Siemens scan
                metaOut.SliceLocation = metaOut.ImagePositionPatient[slice_dim].original_string

        metaOut.PixelData = slice.tobytes()
        if(not dry_run):
            metaOut.save_as(dcmname, write_like_original=False)
        if(verbose):
            print("Saving: %s (file %d of %d)" % (dcmname, instanceNumber, nVols*nSlices))

    # Main loop
    if(nVols>1 and slice_order_4d == 'xyzt'):
        for v in range(nVols):
            # We must redefine the range every time we wish to restart the counter
            range_slices = reversed(range(nSlices)) if reversed_slices else range(nSlices)
            for i in range_slices:
                instanceNumber = instanceNumber + 1
                save_dcm()
    else: # xytz, the more common case
        # We only traverse the slice range one
        range_slices = reversed(range(nSlices)) if reversed_slices else range(nSlices)
        for i in range_slices:
            for v in range(nVols):
                instanceNumber = instanceNumber + 1
                save_dcm()

    print("Finished writing " + str(nVols * nSlices) + " dcm files to outDir: " + outDir)
    if(dry_run):
        print('Dry-run completed, no data were written.')