<filename>vespa/interfaces/inline/vespa_inline_engine.py
﻿
# Python modules
from __future__ import division, print_function, absolute_import

import os
import io
import sys
import datetime
import traceback

# 3rd party modules
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
from matplotlib.backends.backend_pdf import PdfPages

# Our modules
import vespa.analysis.util_import      as util_import
import vespa.analysis.util_file_import as util_file_import
import vespa.analysis.figure_layouts   as figure_layouts
import vespa.common.util.export as util_export
import vespa.common.util.time_  as util_time


# TODO - bjs Future work and thoughts:
# ---------------------------------------------------------
# DICOM Tag for lots of text ...
#
# Need to stack Exception messages to create call stack info as it propagates up the program
#
# Can use any input/output directories, maybe an ini file for users to modify?
#  or std locs if ini not available
#
# can ditch STD mix=1 data, but should modify parser in future to keep in keyword optional
#
# Show work flow for


# Change to True to enable the assert() statements sprinkled through the code
ASSERTIONS_ENABLED = False



DESC =  \
"""
 Command line interface to process MRS data inline with an MR scanner 
 using the Vespa-Analysis package.
  
 Note. You may have to enclose data/preset/output strings in double 
 quotation marks for them to process properly if they have  
 spaces or other special characters embedded in them.
"""


#==============================================================================
# Helper Classes and Methods

class VespaInlineError(Exception):
    """Exception for errors that occur in the Vespa Inline Engine"""
    pass


def _repackage_exception():
    """
    Returns a 3-tuple of (exception type, exception value, trace).

    The type and value are as described in the doc for sys.exc_info().

    However, the traceback has been manipulated a little. It's been turned
    from a traceback object into a list. At some point, some of these
    execptions (and therefore tracebacks) could be generated by code in a
    process other than the main one and returned to the main process as
    results. That means they must be pickleable. Traceback objects can't be
    pickled; doing so raises "TypeError: can't pickle traceback objects".
    Therefore, if we try to return a raw traeback object, pickling will
    fail. Repackaging them as a list retains the traceback info without
    offending the pickle gods.

    """
    lines = []

    type_, value, trace = sys.exc_info()

    # Turn the trace into a list.
    trace = traceback.extract_tb(trace)

    lines = [str(type_),]
    lines.append(str(value))

    # Each entry in the trace list is a 4-tuple --
    #    (filename, line number, function name, text)
    for r in trace:
        lines.append(str(r.filename)+' line:'+str(r.lineno)+', method:'+str(r.name)+', '+str(r.line))

    #trace = [entry for entry in trace]

    return lines



class VespaInlineSettings:

    def __init__(self):
        '''
        This class embodies a list of input/output settings that control the
        behaviour of the Vespa Inline Engine (VIE), and to some extent the scanner
        dependent code module that calls the VIE. It has default values that
        should work for a standard setup, but can be tailored by the user to
        perform for a variety of other workflow locations.

        '''
        self.vespa_version = 'Vespa-CLI'                            # auto-set
        self.base_path = os.path.dirname(os.path.abspath(__file__)) # auto-set
        self.import_class  = 'import_philips_dicom'                 # auto-set
        self.dataformat    = 'philips_press28_dicom'   # user must set

        # Input/Output (absolute) directory paths where data or results reside
        # - coding trick, self.base_path is auto-set on instantiation to
        #    directory in which this module resides. This may also be where
        #    sub-directories reside for presets, output, debug, data, etc.
        # - Ex. self.preset_dir = os.path.join(self.base_path, 'presets')

        self.data_dir   = os.path.join(self.base_path, 'datadir')
        self.preset_dir = os.path.join(self.base_path, 'presets')
        self.output_dir = os.path.join(self.base_path, 'output')
        self.debug_dir  = os.path.join(self.base_path, 'debug')

        # Output flags - control output produced, True means output is created
        #  err – text output of filenames and traceback due to fitting error
        #  xml - full fitting provenance can load into Vespa-Analysis GUI
        #  pdf – hires image of plot/table fitting results
        #  png – fitting result image of plots/table – same as in DICOM output
        #  dcm – ‘screenshot’ fitting result image of plots/table
        #  dcm_pdf - DICOM encapsulated PDF from pdf file created above

        self.save_err = True
        self.save_xml = False
        self.save_pdf = False
        self.save_png = False
        self.save_dcm = True
        self.save_dcm_pdf = False

        # Flags for unique file naming
        #  - if True, append timestamp to filename just before the extension
        #  - timestamp includes: yymmdd.hhmmss._usec_
        #  - Ex. dicom_out_20200603.163207.478473.dcm

        self.err_fname_unique = True
        self.xml_fname_unique = True
        self.pdf_fname_unique = True
        self.png_fname_unique = True
        self.dcm_fname_unique = True
        self.dcm_pdf_fname_unique = True

        # Output filenames
        # - should be absolute path/filename.extn
        # - can use self.debug_dir, self.output_dir to create these
        # - Ex. self.png_fname = os.path.join(self.output_dir,'output_png.png')

        self.err_fname     = os.path.join(self.debug_dir,'default_err.png')
        self.xml_fname     = os.path.join(self.debug_dir,'default.xml')
        self.pdf_fname     = os.path.join(self.debug_dir,'default.pdf')
        self.png_fname     = os.path.join(self.output_dir,'default_output_png.png')
        self.dcm_fname     = os.path.join(self.output_dir,'default_output_dicom.dcm')
        self.dcm_pfd_fname = os.path.join(self.output_dir,'default_output_dicom_pdf.dcm')

        # The following sections set pdf, png and dcm creation parameters.
        # Many correspond to matplotlib keywords (used to create the images)
        #
        # PDF results plot settings

        self.pdf_plotstyle   = 'lcm_multi'
        self.pdf_file_label  = 'default_pdf'
        self.pdf_inbuf       = False
        self.pdf_minppm      = 0.5
        self.pdf_maxppm      = 4.2
        self.pdf_apply_phase = False
        self.pdf_remove_base = False
        self.pdf_fontname    = 'Courier New'
        self.pdf_dpi         = 300
        self.pdf_pad_inches  = 0.5

        # PNG and DICOM results file settings

        self.png_plotstyle   = 'lcm_square'
        self.png_file_label  = 'default_png'
        self.png_inbuf       = False
        self.png_minppm      = 0.5
        self.png_maxppm      = 4.2
        self.png_apply_phase = False
        self.png_remove_base = False
        self.png_fontname    = 'Courier New'
        self.png_dpi         = 100
        self.png_pad_inches  = 0.5

        # Note. The PNG settings also determine the DICOM output. The PNG
        #   'lcm_square' layout Figure is defined to be 10.24 x 10.24 inches
        #   so a png_dpi=100 would give a 1024x1024 DICOM image resolution.

        # Error reporting output settings
        # - These also affect the size of a DICOM image, just a debug readout.
        # - The 'inline_error' layout is define to be 10.24 x 10.24 inches
        # - A err_dpi = 100 would result in a 1024x1024 DICOM image.
        # - If possible, the text is much clearer at err_dpi = 200, but ...

        self.err_dpi         = 300
        self.err_pad_inches  = 0.5

        # Generic flag for other debug actions. Setting to True will output
        # the PNG buffer, that becomes the DICOM result, into a numpy save file
        #  using ndarray.tofile() method in same directory as self.png_fname

        self.debug = False     # generic flag for other debuggery


def is_dicom(filename):
    """
    Returns True if the file in question is a DICOM file, else False.

    Per the DICOM specs, a DICOM file starts with 128 reserved
    bytes followed by "DICM".

    ref: DICOM spec, Part 10: Media Storage and File Format for Media
         Interchange, 7.1 DICOM FILE META INFORMATION

    """
    if os.path.isfile(filename):
        f = open(filename, "rb")
        s = f.read(132)
        f.close()
        pattern = "DICM"
        binary_pattern = pattern.encode()
        return s.endswith(binary_pattern)
    else:
        return False


def get_time():

    now = datetime.datetime.now()
    current_time = now.strftime("%H:%M:%S")
    return current_time


#==============================================================================
# CLI Interface Methods


def analysis_cli_chain(datasets, presets,
                                 basis_mmol,
                                 verbose=False,
                                 process_id='single'):
    
    # Sort datasets into variables --------------------------------------------
    
    data_metab,   data_water,   data_coil,    data_ecc   = datasets
    preset_metab, preset_water, preset_coil,  preset_ecc = presets

    # -------------------------------------------------------------------------
    # Preset/Process Coil Combine Dataset

    if data_coil is not None:
        if verbose: print(process_id + " - apply preset - coil combine")
        # update dataset object with preset blocks and chains
        data_coil.apply_preset(preset_coil, voxel=(0, 0, 0))

        msg = process_id + " - run chain - coil combine"
        if verbose: print(msg)
        _process_all_blocks(data_coil)

    # ----------------------------------------------------------------------
    # Apply presets to ecc, water and metab datasets

    if data_ecc is not None and preset_ecc is not None:
        if verbose: print(process_id+" - apply preset - ecc")
        data_ecc.apply_preset(preset_ecc, voxel=(0,0,0))

    if data_water is not None and preset_water is not None:
        if verbose: print(process_id+" - apply preset - water")
        data_water.apply_preset(preset_water, voxel=(0,0,0))

    if data_metab is not None and preset_metab is not None:
        if verbose: print(process_id+" - apply preset - metab")
        data_metab.apply_preset(preset_metab, voxel=(0,0,0))

    #----------------------------------------------------------------------
    # Attach coil combine to ecc, water and metab datasets - run chain ecc

    if data_coil is not None:
        if verbose: print(process_id+" - attach coil combine to - ecc, water and metab")
        for dset in [data_ecc, data_water, data_metab]:
            if dset is not None:
                dset.set_associated_dataset_combine(data_coil)

    if data_ecc is not None:
        if verbose: print(process_id+" - run chain - ecc")
        _process_all_blocks(data_ecc)       # get combined FID for next steps

    #----------------------------------------------------------------------
    # Attach ecc to water and metab datasets - run chain water

    if data_ecc is not None:
        if verbose: print(process_id+" - attach ecc to - water and metab")
        for dset in [data_water, data_metab]:
            if dset is not None:
                    dset.set_associated_dataset_ecc(data_ecc)

    if data_water is not None:
        if verbose: print(process_id+" - run chain - water")
        _process_all_blocks(data_water)

    #----------------------------------------------------------------------
    # Attach mmol_basis and water to metab dataset - run chain metab

    for dset in [data_metab,]:
        if basis_mmol is not None:
            if dset is not None:
                if verbose: print(process_id + " - attach mmol_basis to - metab")
                dset.set_associated_dataset_mmol(basis_mmol)

        if data_water is not None:
            if verbose: print(process_id + " - attaching water to - metab")
            dset.set_associated_dataset_quant(data_water)

    if verbose: print(process_id+" - run chain - metab")
    _process_all_blocks(data_metab)

    return data_metab
    
    
    
def analysis_cli_output(data_metab, settings, verbose=False, process_id='single'):

    # Test keyword values -----------------------------------------------------

    if settings is None:
        settings = VespaInlineSettings()

    tstamp       = util_time.now(util_time.DISPLAY_TIMESTAMP_FORMAT)
    fname_tstamp = util_time.filename_timestamp()  # yyyymmdd.hhmmss.usecs

    # determine a default output base name, just in case

    if data_metab.dataset_filename != '':
        fdefault = data_metab.dataset_filename
    elif data_metab.blocks['raw'].data_sources[0] != '':
        fdefault = data_metab.blocks['raw'].data_sources[0]
    else:
        fdefault = '.' + os.sep + 'default_out'

    #--------------------------------------------------------------------------
    # Save Dataset to XML for provenance

    if settings.save_xml:

        # Determine if a (unique) output filename exists or can be created ----

        if settings.xml_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            xml_fname = fpath + os.sep + fbase + '.xml'
        else:
            xml_fname, _ = os.path.splitext(settings.xml_fname)
        if settings.xml_fname_unique:
            xml_fname += '_'+fname_tstamp
        xml_fname += '.xml'

        if verbose: print(process_id+" - Saving dataset to XML file %s. " % xml_fname)

        data_metab.dataset_filename = xml_fname
        util_export.export(xml_fname, [data_metab,], None, None, False)

    #--------------------------------------------------------------------------
    # Save results to PDF

    pdf_buf = None
    if settings.save_pdf or settings.save_dcm_pdf:

        fig_call = figure_layouts.null_call # default

        # Determine if a (unique) output filename exists or can be created ----

        if settings.pdf_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            pdf_fname = fpath + os.sep + fbase
        else:
            pdf_fname, _ = os.path.splitext(settings.pdf_fname)
        if settings.pdf_fname_unique:
            pdf_fname += '_'+fname_tstamp

        if settings.dcm_pdf_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext = os.path.splitext(fname)
            dcm_pdf_fname = fpath + os.sep + fbase
        else:
            dcm_pdf_fname, _ = os.path.splitext(settings.dcm_pdf_fname)
        if settings.dcm_pdf_fname_unique:
            dcm_pdf_fname += '_'+fname_tstamp       # make it unique
        settings.dcm_pdf_fname = dcm_pdf_fname+'.dcm'   # store for later use

        # Decide what picture to create ------------------------------------

        if settings.pdf_plotstyle == 'lcm':
            pdf_fname  += '_lcm.pdf'
            fig_call = figure_layouts.lcm_like
        elif settings.pdf_plotstyle == 'lcm_multi':
            pdf_fname += '_lcm_multi.pdf'
            fig_call = figure_layouts.lcm_multipage_pdf

        if verbose: print(process_id+" - Saving Results to PDF %s " % pdf_fname)

        figs = fig_call(data_metab,
                        viffpath=settings.pdf_file_label,
                        vespa_version=settings.vespa_version,
                        timestamp='',
                        fontname=settings.pdf_fontname,
                        minplot=settings.pdf_minppm,
                        maxplot=settings.pdf_maxppm,
                        nobase=settings.pdf_remove_base,
                        extfig=None,
                        fixphase=settings.pdf_apply_phase,
                        verbose=False,
                        debug=False,
                        quantvals=True)

        # Create the PdfPages object to which we will save the pages:
        #  - with statement endsures object closed at end of block, even if Exception

        pdf_buf = io.BytesIO()
        with PdfPages(pdf_buf) as pdf:
            for fig in figs:
                pdf.savefig(fig, dpi=settings.pdf_dpi,
                                 pad_inches=settings.pdf_pad_inches,
                                 facecolor=fig.get_facecolor(),
                                 edgecolor='none')

            # We can also set the file's metadata via the PdfPages object:
            today = datetime.date.today()
            d = pdf.infodict()
            d['Title']        = u'Vespa Output - PDF'
            d['Author']       = u'<NAME>'
            d['Subject']      = u'Vespa results output'
            d['Keywords']     = u'PdfPages Vespa output lcm multi-page'
            d['CreationDate'] = datetime.datetime(today.year, today.month, today.day)
            d['ModDate']      = datetime.datetime.today()

        if settings.save_pdf:
            with open(pdf_fname, "wb") as f:
                f.write(pdf_buf.getvalue())

    # Create PNG output -------------------------------------------

    png_buf = None             # default setting
    if settings.save_png or settings.save_dcm:

        # Determine if (unique) output filenames exist or can be created ------

        if settings.dcm_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            dcm_fname = fpath + os.sep + fbase
        else:
            dcm_fname, _ = os.path.splitext(settings.dcm_fname)
        if settings.dcm_fname_unique:
            dcm_fname += '_'+fname_tstamp       # make it unique
        settings.dcm_fname = dcm_fname+'.dcm'   # store for later use

        if settings.png_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            png_fname = fpath + os.sep + fbase
        else:
            png_fname, _ = os.path.splitext(settings.png_fname)
        if settings.png_fname_unique:
            png_fname += '_'+fname_tstamp       # make it unique

        # Decide what picture to create ---------------------------------------

        if settings.png_plotstyle == 'lcm':
            fig_call = figure_layouts.lcm_like
        elif settings.png_plotstyle == 'lcm_square':
            fig_call = figure_layouts.lcm_square
        elif settings.png_plotstyle == 'lcm_multi':
            fig_call = figure_layouts.lcm_multipage_pdf
        elif settings.png_plotstyle == 'brp_generic':
            fig_call = figure_layouts.analysis_brp_generic

        fig = fig_call( data_metab,
                        viffpath=settings.png_file_label,
                        vespa_version=settings.vespa_version,
                        timestamp=tstamp,
                        fontname=settings.png_fontname,
                        minplot=settings.png_minppm,
                        maxplot=settings.png_maxppm,
                        nobase=settings.png_remove_base,
                        extfig=None,
                        fixphase=settings.png_apply_phase,
                        verbose=False, 
                        debug=False,
                        dpi=settings.png_dpi)

        # convert figure plot to string to byte arry for dicom output ---------

        buf1 = fig[0].canvas.tostring_rgb()
        png_buf = np.fromstring(buf1, dtype=np.uint8)

        if settings.save_png:
            if verbose: print(process_id+" - saving debug PNG file to - "+str(png_fname+'.png'))
            fig[0].savefig( png_fname+'.png',
                            dpi=settings.png_dpi,
                            pad_inches=settings.png_pad_inches)
        if settings.debug:
            if verbose: print(process_id+" - saving debug NPY file to - "+str(png_fname+'_png_buf.npy'))
            png_buf.tofile(png_fname+'_png_buf.npy')

    if verbose: print(process_id+' - finished analysis_cli_output()')

    return png_buf, pdf_buf
            
            
               
def _process_all_blocks(dataset):
    """ for all voxels, run chain in all blocks to update """
    
    chain_outputs = {}
    
    voxel = dataset.all_voxels
    for key in dataset.blocks.keys():
        if key == 'spectral':
            key = 'spectral'
            block = dataset.blocks[key]
            tmp = block.chain.run(voxel, entry='all')
            chain_outputs[key] = tmp
            if 'fit' in dataset.blocks.keys():
                key = 'fit'
                block = dataset.blocks[key]
                block.chain.run(voxel, entry='initial_only')
                key = 'spectral'
                block = dataset.blocks[key]
                block.set_do_fit(True, voxel[0])
                tmp = block.chain.run(voxel, entry='all')
                chain_outputs[key] = tmp
        else:
            block = dataset.blocks[key]
            tmp = block.chain.run(voxel, entry='all')
            chain_outputs[key] = tmp

    return chain_outputs



def _load_preset(presetfile, verbose=False, process_id='single'):
    """ Load Vespa-Analsys preset file """

    if not presetfile:      # required by 'analysis_kernel()'
        return None

    msg = process_id+" - load_preset - fname = %s"  % presetfile
    if verbose: print(msg)
     
    importer = util_import.DatasetImporter(presetfile)

    # Time to rock and roll!
    presets = importer.go()
    preset  = presets[0]

    return preset


def analysis_kernel(params, verbose=False, process_id='single', exception_as_list=False):

    png_buf, pdf_buf = None, None   # defaults
    try:
        fdatasets, fpresets, fbasis_mmol, settings = params

        if verbose:
            print(process_id + ' - Begin - metab filename = ' + fdatasets['metab'])

        msg = 'begin loading presets'
        presets = []
        for key in ['metab','water','ecc','coil']:
            fname = fpresets[key]
            if fname is not None:
                presets.append(_load_preset(fname, verbose=True, process_id=process_id))
            else:
                presets.append(None)

        msg = 'begin loading datasets'
        datasets = []
        if settings.dataformat in ['philips_press28_dicom', 'philips_slaser30_cmrr_spar', 'vasf']:
            for key in ['metab', 'water', 'ecc', 'coil']:
                fname = fdatasets[key]
                if fname is not None:
                    dataset = util_file_import.get_datasets_cli(fname, settings.import_class, None)
                    datasets.append(dataset[0])
                else:
                    datasets.append(None)

        elif settings.dataformat in ['slaser_cmrr',]:
            fname = fdatasets['metab']
            if fname is not None:
                dataset = util_file_import.get_datasets_cli(fname, settings.import_class, None)
                datasets.append(dataset[0],dataset[1],dataset[2],dataset[3])            # all files in one

        else:
            msg = "VIE.analysis_kernel: Unknown 'settings.dataformat' name, can not load dataset, returning."
            raise ValueError(msg)

        msg = 'begin loading mmol basis'
        if fbasis_mmol is not None:
            basis_mmol, msg = util_file_import.open_viff_dataset_file([fbasis_mmol,])
            basis_mmol = basis_mmol[0]
        else:
            basis_mmol = None

        msg = 'begin run cli_chain'
        data_metab = analysis_cli_chain(datasets, presets, basis_mmol,
                                        verbose=verbose,
                                        process_id=process_id)
        msg = 'begin run cli_output'
        png_buf, pdf_buf = analysis_cli_output(data_metab, settings,
                                               verbose=verbose,
                                               process_id=process_id)
        if verbose:
            print(process_id+' - Finish - metab filename = ' + fdatasets['metab'])

    except Exception as e:
        if exception_as_list:
            return png_buf, pdf_buf, _repackage_exception()  # report from multiprocess
        else:
            raise VespaInlineError(_repackage_exception()) # single thread

    return png_buf, pdf_buf, ''
    



def do_main():
    """
    no unit test implemented yet
    
    """
    pass




if __name__ == '__main__':
    
    do_main()



<filename>vespa/interfaces/inline/vespa_inline_engine.py
﻿
# Python modules
from __future__ import division, print_function, absolute_import

import os
import io
import sys
import datetime
import traceback

# 3rd party modules
import numpy as np
import matplotlib as mpl
mpl.use('Agg')
from matplotlib.backends.backend_pdf import PdfPages

# Our modules
import vespa.analysis.util_import      as util_import
import vespa.analysis.util_file_import as util_file_import
import vespa.analysis.figure_layouts   as figure_layouts
import vespa.common.util.export as util_export
import vespa.common.util.time_  as util_time


# TODO - bjs Future work and thoughts:
# ---------------------------------------------------------
# DICOM Tag for lots of text ...
#
# Need to stack Exception messages to create call stack info as it propagates up the program
#
# Can use any input/output directories, maybe an ini file for users to modify?
#  or std locs if ini not available
#
# can ditch STD mix=1 data, but should modify parser in future to keep in keyword optional
#
# Show work flow for


# Change to True to enable the assert() statements sprinkled through the code
ASSERTIONS_ENABLED = False



DESC =  \
"""
 Command line interface to process MRS data inline with an MR scanner 
 using the Vespa-Analysis package.
  
 Note. You may have to enclose data/preset/output strings in double 
 quotation marks for them to process properly if they have  
 spaces or other special characters embedded in them.
"""


#==============================================================================
# Helper Classes and Methods

class VespaInlineError(Exception):
    """Exception for errors that occur in the Vespa Inline Engine"""
    pass


def _repackage_exception():
    """
    Returns a 3-tuple of (exception type, exception value, trace).

    The type and value are as described in the doc for sys.exc_info().

    However, the traceback has been manipulated a little. It's been turned
    from a traceback object into a list. At some point, some of these
    execptions (and therefore tracebacks) could be generated by code in a
    process other than the main one and returned to the main process as
    results. That means they must be pickleable. Traceback objects can't be
    pickled; doing so raises "TypeError: can't pickle traceback objects".
    Therefore, if we try to return a raw traeback object, pickling will
    fail. Repackaging them as a list retains the traceback info without
    offending the pickle gods.

    """
    lines = []

    type_, value, trace = sys.exc_info()

    # Turn the trace into a list.
    trace = traceback.extract_tb(trace)

    lines = [str(type_),]
    lines.append(str(value))

    # Each entry in the trace list is a 4-tuple --
    #    (filename, line number, function name, text)
    for r in trace:
        lines.append(str(r.filename)+' line:'+str(r.lineno)+', method:'+str(r.name)+', '+str(r.line))

    #trace = [entry for entry in trace]

    return lines



class VespaInlineSettings:

    def __init__(self):
        '''
        This class embodies a list of input/output settings that control the
        behaviour of the Vespa Inline Engine (VIE), and to some extent the scanner
        dependent code module that calls the VIE. It has default values that
        should work for a standard setup, but can be tailored by the user to
        perform for a variety of other workflow locations.

        '''
        self.vespa_version = 'Vespa-CLI'                            # auto-set
        self.base_path = os.path.dirname(os.path.abspath(__file__)) # auto-set
        self.import_class  = 'import_philips_dicom'                 # auto-set
        self.dataformat    = 'philips_press28_dicom'   # user must set

        # Input/Output (absolute) directory paths where data or results reside
        # - coding trick, self.base_path is auto-set on instantiation to
        #    directory in which this module resides. This may also be where
        #    sub-directories reside for presets, output, debug, data, etc.
        # - Ex. self.preset_dir = os.path.join(self.base_path, 'presets')

        self.data_dir   = os.path.join(self.base_path, 'datadir')
        self.preset_dir = os.path.join(self.base_path, 'presets')
        self.output_dir = os.path.join(self.base_path, 'output')
        self.debug_dir  = os.path.join(self.base_path, 'debug')

        # Output flags - control output produced, True means output is created
        #  err – text output of filenames and traceback due to fitting error
        #  xml - full fitting provenance can load into Vespa-Analysis GUI
        #  pdf – hires image of plot/table fitting results
        #  png – fitting result image of plots/table – same as in DICOM output
        #  dcm – ‘screenshot’ fitting result image of plots/table
        #  dcm_pdf - DICOM encapsulated PDF from pdf file created above

        self.save_err = True
        self.save_xml = False
        self.save_pdf = False
        self.save_png = False
        self.save_dcm = True
        self.save_dcm_pdf = False

        # Flags for unique file naming
        #  - if True, append timestamp to filename just before the extension
        #  - timestamp includes: yymmdd.hhmmss._usec_
        #  - Ex. dicom_out_20200603.163207.478473.dcm

        self.err_fname_unique = True
        self.xml_fname_unique = True
        self.pdf_fname_unique = True
        self.png_fname_unique = True
        self.dcm_fname_unique = True
        self.dcm_pdf_fname_unique = True

        # Output filenames
        # - should be absolute path/filename.extn
        # - can use self.debug_dir, self.output_dir to create these
        # - Ex. self.png_fname = os.path.join(self.output_dir,'output_png.png')

        self.err_fname     = os.path.join(self.debug_dir,'default_err.png')
        self.xml_fname     = os.path.join(self.debug_dir,'default.xml')
        self.pdf_fname     = os.path.join(self.debug_dir,'default.pdf')
        self.png_fname     = os.path.join(self.output_dir,'default_output_png.png')
        self.dcm_fname     = os.path.join(self.output_dir,'default_output_dicom.dcm')
        self.dcm_pfd_fname = os.path.join(self.output_dir,'default_output_dicom_pdf.dcm')

        # The following sections set pdf, png and dcm creation parameters.
        # Many correspond to matplotlib keywords (used to create the images)
        #
        # PDF results plot settings

        self.pdf_plotstyle   = 'lcm_multi'
        self.pdf_file_label  = 'default_pdf'
        self.pdf_inbuf       = False
        self.pdf_minppm      = 0.5
        self.pdf_maxppm      = 4.2
        self.pdf_apply_phase = False
        self.pdf_remove_base = False
        self.pdf_fontname    = 'Courier New'
        self.pdf_dpi         = 300
        self.pdf_pad_inches  = 0.5

        # PNG and DICOM results file settings

        self.png_plotstyle   = 'lcm_square'
        self.png_file_label  = 'default_png'
        self.png_inbuf       = False
        self.png_minppm      = 0.5
        self.png_maxppm      = 4.2
        self.png_apply_phase = False
        self.png_remove_base = False
        self.png_fontname    = 'Courier New'
        self.png_dpi         = 100
        self.png_pad_inches  = 0.5

        # Note. The PNG settings also determine the DICOM output. The PNG
        #   'lcm_square' layout Figure is defined to be 10.24 x 10.24 inches
        #   so a png_dpi=100 would give a 1024x1024 DICOM image resolution.

        # Error reporting output settings
        # - These also affect the size of a DICOM image, just a debug readout.
        # - The 'inline_error' layout is define to be 10.24 x 10.24 inches
        # - A err_dpi = 100 would result in a 1024x1024 DICOM image.
        # - If possible, the text is much clearer at err_dpi = 200, but ...

        self.err_dpi         = 300
        self.err_pad_inches  = 0.5

        # Generic flag for other debug actions. Setting to True will output
        # the PNG buffer, that becomes the DICOM result, into a numpy save file
        #  using ndarray.tofile() method in same directory as self.png_fname

        self.debug = False     # generic flag for other debuggery


def is_dicom(filename):
    """
    Returns True if the file in question is a DICOM file, else False.

    Per the DICOM specs, a DICOM file starts with 128 reserved
    bytes followed by "DICM".

    ref: DICOM spec, Part 10: Media Storage and File Format for Media
         Interchange, 7.1 DICOM FILE META INFORMATION

    """
    if os.path.isfile(filename):
        f = open(filename, "rb")
        s = f.read(132)
        f.close()
        pattern = "DICM"
        binary_pattern = pattern.encode()
        return s.endswith(binary_pattern)
    else:
        return False


def get_time():

    now = datetime.datetime.now()
    current_time = now.strftime("%H:%M:%S")
    return current_time


#==============================================================================
# CLI Interface Methods


def analysis_cli_chain(datasets, presets,
                                 basis_mmol,
                                 verbose=False,
                                 process_id='single'):
    
    # Sort datasets into variables --------------------------------------------
    
    data_metab,   data_water,   data_coil,    data_ecc   = datasets
    preset_metab, preset_water, preset_coil,  preset_ecc = presets

    # -------------------------------------------------------------------------
    # Preset/Process Coil Combine Dataset

    if data_coil is not None:
        if verbose: print(process_id + " - apply preset - coil combine")
        # update dataset object with preset blocks and chains
        data_coil.apply_preset(preset_coil, voxel=(0, 0, 0))

        msg = process_id + " - run chain - coil combine"
        if verbose: print(msg)
        _process_all_blocks(data_coil)

    # ----------------------------------------------------------------------
    # Apply presets to ecc, water and metab datasets

    if data_ecc is not None and preset_ecc is not None:
        if verbose: print(process_id+" - apply preset - ecc")
        data_ecc.apply_preset(preset_ecc, voxel=(0,0,0))

    if data_water is not None and preset_water is not None:
        if verbose: print(process_id+" - apply preset - water")
        data_water.apply_preset(preset_water, voxel=(0,0,0))

    if data_metab is not None and preset_metab is not None:
        if verbose: print(process_id+" - apply preset - metab")
        data_metab.apply_preset(preset_metab, voxel=(0,0,0))

    #----------------------------------------------------------------------
    # Attach coil combine to ecc, water and metab datasets - run chain ecc

    if data_coil is not None:
        if verbose: print(process_id+" - attach coil combine to - ecc, water and metab")
        for dset in [data_ecc, data_water, data_metab]:
            if dset is not None:
                dset.set_associated_dataset_combine(data_coil)

    if data_ecc is not None:
        if verbose: print(process_id+" - run chain - ecc")
        _process_all_blocks(data_ecc)       # get combined FID for next steps

    #----------------------------------------------------------------------
    # Attach ecc to water and metab datasets - run chain water

    if data_ecc is not None:
        if verbose: print(process_id+" - attach ecc to - water and metab")
        for dset in [data_water, data_metab]:
            if dset is not None:
                    dset.set_associated_dataset_ecc(data_ecc)

    if data_water is not None:
        if verbose: print(process_id+" - run chain - water")
        _process_all_blocks(data_water)

    #----------------------------------------------------------------------
    # Attach mmol_basis and water to metab dataset - run chain metab

    for dset in [data_metab,]:
        if basis_mmol is not None:
            if dset is not None:
                if verbose: print(process_id + " - attach mmol_basis to - metab")
                dset.set_associated_dataset_mmol(basis_mmol)

        if data_water is not None:
            if verbose: print(process_id + " - attaching water to - metab")
            dset.set_associated_dataset_quant(data_water)

    if verbose: print(process_id+" - run chain - metab")
    _process_all_blocks(data_metab)

    return data_metab
    
    
    
def analysis_cli_output(data_metab, settings, verbose=False, process_id='single'):

    # Test keyword values -----------------------------------------------------

    if settings is None:
        settings = VespaInlineSettings()

    tstamp       = util_time.now(util_time.DISPLAY_TIMESTAMP_FORMAT)
    fname_tstamp = util_time.filename_timestamp()  # yyyymmdd.hhmmss.usecs

    # determine a default output base name, just in case

    if data_metab.dataset_filename != '':
        fdefault = data_metab.dataset_filename
    elif data_metab.blocks['raw'].data_sources[0] != '':
        fdefault = data_metab.blocks['raw'].data_sources[0]
    else:
        fdefault = '.' + os.sep + 'default_out'

    #--------------------------------------------------------------------------
    # Save Dataset to XML for provenance

    if settings.save_xml:

        # Determine if a (unique) output filename exists or can be created ----

        if settings.xml_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            xml_fname = fpath + os.sep + fbase + '.xml'
        else:
            xml_fname, _ = os.path.splitext(settings.xml_fname)
        if settings.xml_fname_unique:
            xml_fname += '_'+fname_tstamp
        xml_fname += '.xml'

        if verbose: print(process_id+" - Saving dataset to XML file %s. " % xml_fname)

        data_metab.dataset_filename = xml_fname
        util_export.export(xml_fname, [data_metab,], None, None, False)

    #--------------------------------------------------------------------------
    # Save results to PDF

    pdf_buf = None
    if settings.save_pdf or settings.save_dcm_pdf:

        fig_call = figure_layouts.null_call # default

        # Determine if a (unique) output filename exists or can be created ----

        if settings.pdf_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            pdf_fname = fpath + os.sep + fbase
        else:
            pdf_fname, _ = os.path.splitext(settings.pdf_fname)
        if settings.pdf_fname_unique:
            pdf_fname += '_'+fname_tstamp

        if settings.dcm_pdf_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext = os.path.splitext(fname)
            dcm_pdf_fname = fpath + os.sep + fbase
        else:
            dcm_pdf_fname, _ = os.path.splitext(settings.dcm_pdf_fname)
        if settings.dcm_pdf_fname_unique:
            dcm_pdf_fname += '_'+fname_tstamp       # make it unique
        settings.dcm_pdf_fname = dcm_pdf_fname+'.dcm'   # store for later use

        # Decide what picture to create ------------------------------------

        if settings.pdf_plotstyle == 'lcm':
            pdf_fname  += '_lcm.pdf'
            fig_call = figure_layouts.lcm_like
        elif settings.pdf_plotstyle == 'lcm_multi':
            pdf_fname += '_lcm_multi.pdf'
            fig_call = figure_layouts.lcm_multipage_pdf

        if verbose: print(process_id+" - Saving Results to PDF %s " % pdf_fname)

        figs = fig_call(data_metab,
                        viffpath=settings.pdf_file_label,
                        vespa_version=settings.vespa_version,
                        timestamp='',
                        fontname=settings.pdf_fontname,
                        minplot=settings.pdf_minppm,
                        maxplot=settings.pdf_maxppm,
                        nobase=settings.pdf_remove_base,
                        extfig=None,
                        fixphase=settings.pdf_apply_phase,
                        verbose=False,
                        debug=False,
                        quantvals=True)

        # Create the PdfPages object to which we will save the pages:
        #  - with statement endsures object closed at end of block, even if Exception

        pdf_buf = io.BytesIO()
        with PdfPages(pdf_buf) as pdf:
            for fig in figs:
                pdf.savefig(fig, dpi=settings.pdf_dpi,
                                 pad_inches=settings.pdf_pad_inches,
                                 facecolor=fig.get_facecolor(),
                                 edgecolor='none')

            # We can also set the file's metadata via the PdfPages object:
            today = datetime.date.today()
            d = pdf.infodict()
            d['Title']        = u'Vespa Output - PDF'
            d['Author']       = u'<NAME>'
            d['Subject']      = u'Vespa results output'
            d['Keywords']     = u'PdfPages Vespa output lcm multi-page'
            d['CreationDate'] = datetime.datetime(today.year, today.month, today.day)
            d['ModDate']      = datetime.datetime.today()

        if settings.save_pdf:
            with open(pdf_fname, "wb") as f:
                f.write(pdf_buf.getvalue())

    # Create PNG output -------------------------------------------

    png_buf = None             # default setting
    if settings.save_png or settings.save_dcm:

        # Determine if (unique) output filenames exist or can be created ------

        if settings.dcm_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            dcm_fname = fpath + os.sep + fbase
        else:
            dcm_fname, _ = os.path.splitext(settings.dcm_fname)
        if settings.dcm_fname_unique:
            dcm_fname += '_'+fname_tstamp       # make it unique
        settings.dcm_fname = dcm_fname+'.dcm'   # store for later use

        if settings.png_fname == '':
            fpath, fname = os.path.split(os.path.abspath(fdefault))
            fbase, fext  = os.path.splitext(fname)
            png_fname = fpath + os.sep + fbase
        else:
            png_fname, _ = os.path.splitext(settings.png_fname)
        if settings.png_fname_unique:
            png_fname += '_'+fname_tstamp       # make it unique

        # Decide what picture to create ---------------------------------------

        if settings.png_plotstyle == 'lcm':
            fig_call = figure_layouts.lcm_like
        elif settings.png_plotstyle == 'lcm_square':
            fig_call = figure_layouts.lcm_square
        elif settings.png_plotstyle == 'lcm_multi':
            fig_call = figure_layouts.lcm_multipage_pdf
        elif settings.png_plotstyle == 'brp_generic':
            fig_call = figure_layouts.analysis_brp_generic

        fig = fig_call( data_metab,
                        viffpath=settings.png_file_label,
                        vespa_version=settings.vespa_version,
                        timestamp=tstamp,
                        fontname=settings.png_fontname,
                        minplot=settings.png_minppm,
                        maxplot=settings.png_maxppm,
                        nobase=settings.png_remove_base,
                        extfig=None,
                        fixphase=settings.png_apply_phase,
                        verbose=False, 
                        debug=False,
                        dpi=settings.png_dpi)

        # convert figure plot to string to byte arry for dicom output ---------

        buf1 = fig[0].canvas.tostring_rgb()
        png_buf = np.fromstring(buf1, dtype=np.uint8)

        if settings.save_png:
            if verbose: print(process_id+" - saving debug PNG file to - "+str(png_fname+'.png'))
            fig[0].savefig( png_fname+'.png',
                            dpi=settings.png_dpi,
                            pad_inches=settings.png_pad_inches)
        if settings.debug:
            if verbose: print(process_id+" - saving debug NPY file to - "+str(png_fname+'_png_buf.npy'))
            png_buf.tofile(png_fname+'_png_buf.npy')

    if verbose: print(process_id+' - finished analysis_cli_output()')

    return png_buf, pdf_buf
            
            
               
def _process_all_blocks(dataset):
    """ for all voxels, run chain in all blocks to update """
    
    chain_outputs = {}
    
    voxel = dataset.all_voxels
    for key in dataset.blocks.keys():
        if key == 'spectral':
            key = 'spectral'
            block = dataset.blocks[key]
            tmp = block.chain.run(voxel, entry='all')
            chain_outputs[key] = tmp
            if 'fit' in dataset.blocks.keys():
                key = 'fit'
                block = dataset.blocks[key]
                block.chain.run(voxel, entry='initial_only')
                key = 'spectral'
                block = dataset.blocks[key]
                block.set_do_fit(True, voxel[0])
                tmp = block.chain.run(voxel, entry='all')
                chain_outputs[key] = tmp
        else:
            block = dataset.blocks[key]
            tmp = block.chain.run(voxel, entry='all')
            chain_outputs[key] = tmp

    return chain_outputs



def _load_preset(presetfile, verbose=False, process_id='single'):
    """ Load Vespa-Analsys preset file """

    if not presetfile:      # required by 'analysis_kernel()'
        return None

    msg = process_id+" - load_preset - fname = %s"  % presetfile
    if verbose: print(msg)
     
    importer = util_import.DatasetImporter(presetfile)

    # Time to rock and roll!
    presets = importer.go()
    preset  = presets[0]

    return preset


def analysis_kernel(params, verbose=False, process_id='single', exception_as_list=False):

    png_buf, pdf_buf = None, None   # defaults
    try:
        fdatasets, fpresets, fbasis_mmol, settings = params

        if verbose:
            print(process_id + ' - Begin - metab filename = ' + fdatasets['metab'])

        msg = 'begin loading presets'
        presets = []
        for key in ['metab','water','ecc','coil']:
            fname = fpresets[key]
            if fname is not None:
                presets.append(_load_preset(fname, verbose=True, process_id=process_id))
            else:
                presets.append(None)

        msg = 'begin loading datasets'
        datasets = []
        if settings.dataformat in ['philips_press28_dicom', 'philips_slaser30_cmrr_spar', 'vasf']:
            for key in ['metab', 'water', 'ecc', 'coil']:
                fname = fdatasets[key]
                if fname is not None:
                    dataset = util_file_import.get_datasets_cli(fname, settings.import_class, None)
                    datasets.append(dataset[0])
                else:
                    datasets.append(None)

        elif settings.dataformat in ['slaser_cmrr',]:
            fname = fdatasets['metab']
            if fname is not None:
                dataset = util_file_import.get_datasets_cli(fname, settings.import_class, None)
                datasets.append(dataset[0],dataset[1],dataset[2],dataset[3])            # all files in one

        else:
            msg = "VIE.analysis_kernel: Unknown 'settings.dataformat' name, can not load dataset, returning."
            raise ValueError(msg)

        msg = 'begin loading mmol basis'
        if fbasis_mmol is not None:
            basis_mmol, msg = util_file_import.open_viff_dataset_file([fbasis_mmol,])
            basis_mmol = basis_mmol[0]
        else:
            basis_mmol = None

        msg = 'begin run cli_chain'
        data_metab = analysis_cli_chain(datasets, presets, basis_mmol,
                                        verbose=verbose,
                                        process_id=process_id)
        msg = 'begin run cli_output'
        png_buf, pdf_buf = analysis_cli_output(data_metab, settings,
                                               verbose=verbose,
                                               process_id=process_id)
        if verbose:
            print(process_id+' - Finish - metab filename = ' + fdatasets['metab'])

    except Exception as e:
        if exception_as_list:
            return png_buf, pdf_buf, _repackage_exception()  # report from multiprocess
        else:
            raise VespaInlineError(_repackage_exception()) # single thread

    return png_buf, pdf_buf, ''
    



def do_main():
    """
    no unit test implemented yet
    
    """
    pass




if __name__ == '__main__':
    
    do_main()



# Python modules # 3rd party modules # Our modules # TODO - bjs Future work and thoughts: # --------------------------------------------------------- # DICOM Tag for lots of text ... # # Need to stack Exception messages to create call stack info as it propagates up the program # # Can use any input/output directories, maybe an ini file for users to modify? #  or std locs if ini not available # # can ditch STD mix=1 data, but should modify parser in future to keep in keyword optional # # Show work flow for # Change to True to enable the assert() statements sprinkled through the code Command line interface to process MRS data inline with an MR scanner   using the Vespa-Analysis package.     Note. You may have to enclose data/preset/output strings in double   quotation marks for them to process properly if they have    spaces or other special characters embedded in them. #============================================================================== # Helper Classes and Methods Exception for errors that occur in the Vespa Inline Engine Returns a 3-tuple of (exception type, exception value, trace).      The type and value are as described in the doc for sys.exc_info().      However, the traceback has been manipulated a little. It's been turned     from a traceback object into a list. At some point, some of these     execptions (and therefore tracebacks) could be generated by code in a     process other than the main one and returned to the main process as     results. That means they must be pickleable. Traceback objects can't be     pickled; doing so raises "TypeError: can't pickle traceback objects".     Therefore, if we try to return a raw traeback object, pickling will     fail. Repackaging them as a list retains the traceback info without     offending the pickle gods. # Turn the trace into a list. # Each entry in the trace list is a 4-tuple -- #    (filename, line number, function name, text) #trace = [entry for entry in trace] This class embodies a list of input/output settings that control the         behaviour of the Vespa Inline Engine (VIE), and to some extent the scanner         dependent code module that calls the VIE. It has default values that         should work for a standard setup, but can be tailored by the user to         perform for a variety of other workflow locations. # auto-set # auto-set # auto-set # user must set # Input/Output (absolute) directory paths where data or results reside # - coding trick, self.base_path is auto-set on instantiation to #    directory in which this module resides. This may also be where #    sub-directories reside for presets, output, debug, data, etc. # - Ex. self.preset_dir = os.path.join(self.base_path, 'presets') # Output flags - control output produced, True means output is created #  err – text output of filenames and traceback due to fitting error #  xml - full fitting provenance can load into Vespa-Analysis GUI #  pdf – hires image of plot/table fitting results #  png – fitting result image of plots/table – same as in DICOM output #  dcm – ‘screenshot’ fitting result image of plots/table #  dcm_pdf - DICOM encapsulated PDF from pdf file created above # Flags for unique file naming #  - if True, append timestamp to filename just before the extension #  - timestamp includes: yymmdd.hhmmss._usec_ #  - Ex. dicom_out_20200603.163207.478473.dcm # Output filenames # - should be absolute path/filename.extn # - can use self.debug_dir, self.output_dir to create these # - Ex. self.png_fname = os.path.join(self.output_dir,'output_png.png') # The following sections set pdf, png and dcm creation parameters. # Many correspond to matplotlib keywords (used to create the images) # # PDF results plot settings # PNG and DICOM results file settings # Note. The PNG settings also determine the DICOM output. The PNG #   'lcm_square' layout Figure is defined to be 10.24 x 10.24 inches #   so a png_dpi=100 would give a 1024x1024 DICOM image resolution. # Error reporting output settings # - These also affect the size of a DICOM image, just a debug readout. # - The 'inline_error' layout is define to be 10.24 x 10.24 inches # - A err_dpi = 100 would result in a 1024x1024 DICOM image. # - If possible, the text is much clearer at err_dpi = 200, but ... # Generic flag for other debug actions. Setting to True will output # the PNG buffer, that becomes the DICOM result, into a numpy save file #  using ndarray.tofile() method in same directory as self.png_fname # generic flag for other debuggery Returns True if the file in question is a DICOM file, else False.      Per the DICOM specs, a DICOM file starts with 128 reserved     bytes followed by "DICM".      ref: DICOM spec, Part 10: Media Storage and File Format for Media          Interchange, 7.1 DICOM FILE META INFORMATION #============================================================================== # CLI Interface Methods # Sort datasets into variables -------------------------------------------- # ------------------------------------------------------------------------- # Preset/Process Coil Combine Dataset # update dataset object with preset blocks and chains # ---------------------------------------------------------------------- # Apply presets to ecc, water and metab datasets #---------------------------------------------------------------------- # Attach coil combine to ecc, water and metab datasets - run chain ecc # get combined FID for next steps #---------------------------------------------------------------------- # Attach ecc to water and metab datasets - run chain water #---------------------------------------------------------------------- # Attach mmol_basis and water to metab dataset - run chain metab # Test keyword values ----------------------------------------------------- # yyyymmdd.hhmmss.usecs # determine a default output base name, just in case #-------------------------------------------------------------------------- # Save Dataset to XML for provenance # Determine if a (unique) output filename exists or can be created ---- #-------------------------------------------------------------------------- # Save results to PDF # default # Determine if a (unique) output filename exists or can be created ---- # make it unique # store for later use # Decide what picture to create ------------------------------------ # Create the PdfPages object to which we will save the pages: #  - with statement endsures object closed at end of block, even if Exception # We can also set the file's metadata via the PdfPages object: # Create PNG output ------------------------------------------- # default setting # Determine if (unique) output filenames exist or can be created ------ # make it unique # store for later use # make it unique # Decide what picture to create --------------------------------------- # convert figure plot to string to byte arry for dicom output --------- for all voxels, run chain in all blocks to update Load Vespa-Analsys preset file # required by 'analysis_kernel()' # Time to rock and roll! # defaults # all files in one # report from multiprocess # single thread no unit test implemented yet