aux for signa

install_hdf5.md

@@ -0,0 +1,17 @@
+Download source from:
+https://www.hdfgroup.org/downloads/hdf5/source-code/
+
+e.g.:
+https://www.hdfgroup.org/package/hdf5-1-10-6-tar-bz2/?wpdmdl=14134&refresh=5e34a41db4c8c1580508189
+
+decompress to a folder, e.g., /home/user
+
+./configure --prefix=/usr/local/hdf5 <more configure_flags>
+
+make
+
+make check                # run test suite.
+
+make install
+
+make check-install        # verify installation.

niftypet/nipet/CMakeLists.txt

@@ -34,3 +34,12 @@ install(EXPORT ${PROJECT_NAME}Targets FILE NiftyPET${PROJECT_NAME}Targets.cmake
 add_subdirectory(lm)
 add_subdirectory(prj)
 add_subdirectory(sct)
+# HDF5 needed for GE data
+set(HDF5_ROOT "/usr/local/hdf5/")
+find_package(HDF5)
+if(HDF5_FOUND)
+  message(STATUS "HDF5 include dir: ${HDF5_INCLUDE_DIRS}")
+  add_subdirectory(lm_sig)
+else()
+  message(WARNING "Could not find HDF5. Not building lm_sig")
+endif()

niftypet/nipet/__init__.py

@@ -32,6 +32,8 @@
     'back_prj', 'frwd_prj', 'simulate_recon', 'simulate_sino',
     # sct
     'vsm',
+    # signa
+    'aux_sig', 'lm_sig', 'lminfo_sig',
     # optional
     'video_dyn', 'video_frm', 'xnat']  # yapf: disable
 from pkg_resources import resource_filename
@@ -72,10 +74,8 @@
     xnat = None
 
 # > GE Signa
-# from . import aux_sig
-
-# from . import lm_sig
-# from .lm_sig.hst_sig import lminfo_sig
+from . import aux_sig, lm_sig
+from .lm_sig.hst_sig import lminfo_sig
 
 # for use in `cmake -DCMAKE_PREFIX_PATH=...`
 cmake_prefix = resource_filename(__name__, "cmake")

niftypet/nipet/aux_sig.py

@@ -0,0 +1,253 @@
+import logging
+import os
+from math import pi
+
+import h5py
+import numpy as np
+
+log = logging.getLogger(__name__)
+
+
+def constants_h5(pthfn):
+    with h5py.File(pthfn, 'r') as f:
+        # coincidence event mode
+        cncdmd = f['HeaderData']['AcqParameters']['EDCATParameters']['coinOutputMode'][0]
+        if cncdmd == 802:
+            # bytes per event in this mode:
+            bpe = 6
+            log.info("list-mode data in NOMINAL mode (6 bytes per event)")
+        elif cncdmd == 803:
+            bpe = 16
+            log.error(
+                "list-mode data in CALIBRATION mode (16 bytes per event) not currently supported")
+        elif cncdmd == 805:
+            bpe = 8
+            log.error(
+                "the ist-mode data in ENERGY mode (8 bytes per event) not currently supported")
+        else:
+            bpe = 0
+            log.error("list-mode data in UNKNOWN mode")
+
+        # toff: scan start time marker (used as offset)
+        CntH5 = {
+            'toff': f['HeaderData']['AcqStats']['frameStartCoincTStamp'][0],
+            'Deff': f['HeaderData']['SystemGeometry']['effectiveRingDiameter'][0],
+            'TFOV': f['HeaderData']['AcqParameters']['EDCATParameters']['transAxialFOV'][0],
+            'cpitch': f['HeaderData']['SystemGeometry']['interCrystalPitch'][0],
+            'bpitch': f['HeaderData']['SystemGeometry']['interBlockPitch'][0],
+            'exLOR': f['HeaderData']['AcqParameters']['RxScanParameters']['extraRsForTFOV'][0],
+            'axCB': f['HeaderData']['SystemGeometry']['axialCrystalsPerBlock'][0],
+            'axBU': f['HeaderData']['SystemGeometry']['axialBlocksPerUnit'][0],
+            'axUM': f['HeaderData']['SystemGeometry']['axialUnitsPerModule'][0],
+            'axMno': f['HeaderData']['SystemGeometry']['axialModulesPerSystem'][0],
+            'txCB': f['HeaderData']['SystemGeometry']['radialCrystalsPerBlock'][0],
+            'txBU': f['HeaderData']['SystemGeometry']['radialBlocksPerUnit'][0],
+            'txUM': f['HeaderData']['SystemGeometry']['radialUnitsPerModule'][0],
+            'txMno': f['HeaderData']['SystemGeometry']['radialModulesPerSystem'][0],
+            'MRD': f['HeaderData']['AcqParameters']['BackEndAcqFilters']['maxRingDiff'][0],
+            'tau0': f['HeaderData']['AcqParameters']['EDCATParameters']['negCoincidenceWindow'][0],
+            'tau1': f['HeaderData']['AcqParameters']['EDCATParameters']['posCoincidenceWindow'][0],
+            'tauP': f['HeaderData']['AcqParameters']['EDCATParameters']['coincTimingPrecision'][0],
+            'TOFC': f['HeaderData']['AcqParameters']['RxScanParameters']['tofCompressionFactor']
+            [0],
+            'LLD': f['HeaderData']['AcqParameters']['EDCATParameters']['lower_energy_limit'][0],
+            'ULD': f['HeaderData']['AcqParameters']['EDCATParameters']['upper_energy_limit'][0],
+            'BPE': bpe}
+    return CntH5
+
+
+def get_nbins(Cnt):
+    txUno = Cnt['txUM'] * Cnt['txMno']
+    txCU = Cnt['txCB'] * Cnt['txBU']
+    cpitch = Cnt['cpitch']
+    minValue = np.ceil(2 * (np.arcsin(
+        (Cnt['TFOV'] * 10.0) / Cnt['Deff']) - np.floor(txUno * np.arcsin(
+            (Cnt['TFOV'] * 10.0) / Cnt['Deff']) / pi) * pi / txUno) / cpitch)
+    if txCU < minValue:
+        minValue = txCU
+    halfFanLORs = np.floor(txUno * np.arcsin(
+        (Cnt['TFOV'] * 10.0) / Cnt['Deff']) / pi) * txCU + minValue
+    W = 2 * int(halfFanLORs) + 2 * Cnt['exLOR'] + 1
+    C = Cnt['txCB'] * Cnt['txBU'] * Cnt['txUM'] * Cnt['txMno'] * Cnt['axCB']
+    if W > C:
+        W = C
+    return W
+
+
+# ===================================================================================
+# SCANNER CONSTANTS
+def get_sig_constants(pthfn):
+
+    if not os.path.isfile(pthfn):
+        print('e> coult not open the file HDF5 to get SIGNA constants')
+        return
+
+    Cnt = constants_h5(pthfn)
+
+    # default logging set to WARNING only (30)
+    Cnt['LOG'] = 30
+
+    # number of sinogram angles
+    NSBINS = get_nbins(Cnt)
+    # number of transxial crystals
+    NCRS = Cnt['txCB'] * Cnt['txBU'] * Cnt['txUM'] * Cnt['txMno']
+    # number of rings
+    NRNG = Cnt['axCB'] * Cnt['axBU'] * Cnt['axUM'] * Cnt['axMno']
+    # number of 2D sinograms
+    NSN = NRNG**2 - (NRNG-1)
+    # bootstrapping of the list-mode data;
+    # 0: None, 1: Not used (was non-parametric for the mMR), 2: parametric
+    Cnt['BTP'] = 0
+    Cnt['NCRS'] = NCRS
+    Cnt['NRNG'] = NRNG
+    # defnie reduced detector ring
+    Cnt['RNG_END'] = NRNG
+    Cnt['RNG_STRT'] = 0
+    Cnt['NSN'] = NSN
+    Cnt['NSBINS'] = NSBINS
+    Cnt['NSANGLES'] = NCRS // 2
+
+    Cnt['NSEG0'] = 2 * Cnt['axUM'] * Cnt['axCB'] - 1
+
+    # LM processing
+    # integration time of 1 sec
+    Cnt['ITIME'] = 1000
+    # number of CUDA streams
+    Cnt['NSTREAMS'] = 32
+    # number of elements per data chunk
+    # 2^{28} = 268435456 elements (6Bytes) to make up 1.6GB
+    Cnt['ELECHNK'] = (268435456 // Cnt['NSTREAMS'])
+
+    # projection view integration time (length of the short time frames t = 2^VTIME)
+    Cnt['VTIME'] = 2
+    # the short time frame projection views are limited to 90 mins only
+    Cnt['MXNITAG'] = 5400
+
+    # transaxial and axial crystal sizes in cm
+    Cnt['TXCRS'] = 0.395
+    Cnt['AXCRS'] = 0.530
+    # gap between axial detector units
+    Cnt['AXGAP'] = 0.280
+    # axial FOV
+    Cnt['AXFOV'] = (Cnt['axUM'] - 1) * Cnt['AXGAP'] + Cnt['axUM'] * Cnt['axCB'] * Cnt['AXCRS']
+
+    return Cnt
+
+
+# ===================================================================================
+# AXIAL LUTS
+def get_axLUT(Cnt):
+
+    # calculated rings
+    NRNG_c = Cnt['RNG_END'] - Cnt['RNG_STRT']
+    NSN1_c = NRNG_c**2 - (NRNG_c-1)
+
+    # get the sino LUTs
+    M = np.zeros((NRNG_c, NRNG_c), dtype=np.int16)
+    # sino index
+    Msn = np.zeros((NRNG_c, NRNG_c), dtype=np.int16)
+    # sino index SIGNA native
+    Msig = np.zeros((NRNG_c, NRNG_c), dtype=np.int16)
+
+    sn_rno = np.zeros((NRNG_c**2, 2), dtype=np.int16)
+
+    # diagonal linear index (positive only)
+    dli = 0
+    # full diagonal linear index (positive and negative)
+    sni = 0
+
+    for ro in range(0, NRNG_c):
+        if ro == 0:
+            oblique = 1
+        else:
+            oblique = 2
+        for m in range(oblique):
+            # selects the sub-michelogram of the whole michelogram
+            strt = Cnt['NRNG'] * (ro + Cnt['RNG_STRT']) + Cnt['RNG_STRT']
+            stop = (Cnt['RNG_STRT'] + NRNG_c) * Cnt['NRNG']
+            step = Cnt['NRNG'] + 1
+
+            for li in range(strt, stop, step):
+                # goes along a diagonal started in the first row at r2o
+                # from the linear indices of michelogram get the subscript indices
+                if m == 0:
+                    r0 = li // Cnt['NRNG']
+                    r1 = li - r0 * Cnt['NRNG']
+                    M[r1, r0] = dli
+                    dli += 1
+                else:
+                    r1 = li // Cnt['NRNG']
+                    r0 = li - r1 * Cnt['NRNG']
+
+                Msn[r1, r0] = sni
+
+                sn_rno[sni, 0] = r0
+                sn_rno[sni, 1] = r1
+
+                sni += 1
+
+                # --- SIGNA native ---
+                rdiff = r0 - r1
+                rsum = r0 + r1
+                if (rdiff > 1):
+                    angle = rdiff // 2
+                    if (angle <= Cnt['MRD'] / 2):
+                        snis = rsum + (4*angle - 2) * Cnt['NRNG'] - (4*angle*angle - 1)
+                elif (rdiff < -1):
+                    angle = -rdiff // 2
+                    if (angle <= Cnt['MRD'] // 2):
+                        snis = rsum + (4*angle) * Cnt['NRNG'] - ((angle+1) * 4 * angle)
+                else:
+                    snis = rsum
+                Msig[r1, r0] = snis
+                # ------
+    axLUT = {'r2s': Msn, 'r2sig': Msig, 's2r': sn_rno}
+    return axLUT
+
+
+# ===================================================================================
+# TRANSIAXIAL LUTS
+def get_txLUT(Cnt):
+    # number of bins per sinogram angle
+    NSBINS = Cnt['NSBINS']
+    # number of sinogram angles
+    NSANGLES = Cnt['NSANGLES']
+    # number of rings
+    NCRS = Cnt['NCRS']
+    # crystal to sinogram index lookup table (LUT)
+    c2s = np.zeros((NCRS, NCRS), dtype=np.int32)
+    # sinogram to crystal index LUT
+    s2c = np.zeros((NSANGLES * NSBINS, 2), dtype=np.int16)
+    for c0 in range(NCRS):
+        for c1 in range(NCRS):
+            if ((NCRS // 2) <= (c0 + c1)) and ((c0 + c1) < (3 * NCRS // 2)):
+                iw = (NSBINS-1) // 2 + (c0 - c1 - NCRS//2)
+            else:
+                iw = (NSBINS-1) // 2 - (c0 - c1 - NCRS//2)
+            if (iw >= 0) and (iw <= (NSBINS - 1)):
+                ia = ((c0 + c1 + NCRS//2) % NCRS) // 2
+                aw = ia + NSANGLES*iw
+                c2s[c1, c0] = aw
+                c2s[c0, c1] = aw
+                s2c[aw, 0] = c0
+                s2c[aw, 1] = c1
+            else:
+                c2s[c1, c0] = -1
+                c2s[c0, c1] = -1
+
+    txLUT = {'c2s': c2s, 's2c': s2c}
+    return txLUT
+
+
+def init_sig(pthfn):
+
+    # get the constants for the mMR
+    Cnt = get_sig_constants(pthfn)
+
+    # transaxial look up tables
+    txLUT = get_txLUT(Cnt)
+
+    # axial look up tables
+    axLUT = get_axLUT(Cnt)
+
+    return Cnt, txLUT, axLUT

niftypet/nipet/img/mmrimg.py

@@ -155,9 +155,8 @@ def getinterfile_off(fmu, Cnt, Offst=OFFSET_DEFAULT):
     accounting for image offset (does slow interpolation).
     '''
     # pead the image file
-    f = open(fmu, 'rb')
-    mu = np.fromfile(f, np.float32)
-    f.close()
+    with open(fmu, 'rb') as f:
+        mu = np.fromfile(f, np.float32)
 
     # save_im(mur, Cnt, os.path.dirname(fmu) + '/mur.nii')
     # -------------------------------------------------------------------------
@@ -179,9 +178,8 @@ def getinterfile_off(fmu, Cnt, Offst=OFFSET_DEFAULT):
 def getinterfile(fim, Cnt):
     '''Return the floating point image file in an array from an Interfile file.'''
     # pead the image file
-    f = open(fim, 'rb')
-    im = np.fromfile(f, np.float32)
-    f.close()
+    with open(fim, 'rb') as f:
+        im = np.fromfile(f, np.float32)
 
     # pumber of voxels
     nvx = im.shape[0]
@@ -1016,19 +1014,17 @@ def hmu_offset(hdr):
 
 def rd_hmu(fh):
     # --read hdr file--
-    f = open(fh, 'r')
-    hdr = f.read()
-    f.close()
+    with open(fh, 'r') as f:
+        hdr = f.read()
     # -----------------
     # pegular expression to find the file name
     p = re.compile(r'(?<=:=)\s*\w*[.]\w*')
     i0 = hdr.find('!name of data file')
     i1 = i0 + hdr[i0:].find('\n')
     fbin = p.findall(hdr[i0:i1])[0]
     # --read img file--
-    f = open(os.path.join(os.path.dirname(fh), fbin.strip()), 'rb')
-    im = np.fromfile(f, np.float32)
-    f.close()
+    with open(os.path.join(os.path.dirname(fh), fbin.strip()), 'rb') as f:
+        im = np.fromfile(f, np.float32)
     # -----------------
     return hdr, im
 

niftypet/nipet/lm_sig/CMakeLists.txt

@@ -0,0 +1,28 @@
+project(lmproc_sig)
+
+file(GLOB SRC LIST_DIRECTORIES false "src/*.cu")
+include_directories(src)
+include_directories(${Python3_INCLUDE_DIRS})
+include_directories(${Python3_NumPy_INCLUDE_DIRS})
+include_directories(${HDF5_INCLUDE_DIRS})
+
+add_library(${PROJECT_NAME} ${LIB_TYPE} ${SRC})
+add_library(NiftyPET::${PROJECT_NAME} ALIAS ${PROJECT_NAME})
+target_include_directories(${PROJECT_NAME} PUBLIC
+  "$<BUILD_INTERFACE:${${CMAKE_PROJECT_NAME}_INCLUDE_DIRS}>"
+  "$<INSTALL_INTERFACE:niftypet/${CMAKE_PROJECT_NAME}/include>")
+target_link_libraries(${PROJECT_NAME} mmr_auxe ${Python3_LIBRARIES} ${HDF5_LIBRARIES} CUDA::cudart_static CUDA::curand_static)
+
+if(SKBUILD)
+python_extension_module(${PROJECT_NAME})
+endif()
+set_target_properties(${PROJECT_NAME} PROPERTIES
+  VERSION ${CMAKE_PROJECT_VERSION}
+  SOVERSION ${CMAKE_PROJECT_VERSION_MAJOR}
+  INTERFACE_${PROJECT_NAME}_MAJOR_VERSION ${CMAKE_PROJECT_VERSION_MAJOR})
+set_property(TARGET ${PROJECT_NAME} APPEND PROPERTY COMPATIBLE_INTERFACE_STRING ${PROJECT_NAME}_MAJOR_VERSION)
+install(TARGETS ${PROJECT_NAME} EXPORT ${PROJECT_NAME}Targets
+  INCLUDES DESTINATION niftypet/${CMAKE_PROJECT_NAME}/include
+  LIBRARY DESTINATION niftypet/${CMAKE_PROJECT_NAME}/lm_sig)
+install(EXPORT ${PROJECT_NAME}Targets FILE NiftyPET${PROJECT_NAME}Targets.cmake
+  NAMESPACE NiftyPET:: DESTINATION niftypet/${CMAKE_PROJECT_NAME}/cmake)

niftypet/nipet/lm_sig/__init__.py

@@ -0,0 +1,2 @@
+__all__ = ['hst_sig', 'lmproc_sig']
+from . import hst_sig, lmproc_sig