AIQ verified version of NPZ generation

pythonScripts/generateNPZ.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 11,
    "id": "1da1e1c8-1f5a-4c9e-aa4a-1baaa8c71df1",
    "metadata": {},
    "outputs": [
@@ -29,6 +29,8 @@
     "import importlib\n",
     "import numpy\n",
     "import random\n",
+    "import SimpleITK\n",
+    "import datetime\n",
     "\n",
     "nixSuite=os.path.join(os.path.expanduser('~'),'software','src','nixSuite')\n",
     "sys.path.append(os.path.join(nixSuite,'wrapper'))\n",
@@ -47,7 +49,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 2,
    "id": "40fce609-1142-4746-a912-70481ff809be",
    "metadata": {},
    "outputs": [],
@@ -73,7 +75,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 3,
    "id": "e7386af6-7b00-4607-a05e-5581db804869",
    "metadata": {},
    "outputs": [
@@ -81,7 +83,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "User: andrej studen CSRF: 9b381ce234d4a838c53c6f93dab2f5a7\n"
+      "User: andrej studen CSRF: f7abf32866c67a6f503e050311032978\n"
      ]
     }
    ],
@@ -92,7 +94,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 69,
+   "execution_count": 35,
    "id": "5910e7dd-de77-4260-a23a-78e2a6baabb4",
    "metadata": {},
    "outputs": [
@@ -100,8 +102,16 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[-592.9801940917969, -867.9606628417969, -2823.0]\n",
-      "[-1644.0222664799999, -2395.3331437499996, -4710.0]\n"
+      "Returning 1\n",
+      "0.0003862617854116126\n",
+      "sliceThickness=2.0 +- 0.0\n",
+      "928\n",
+      "[-389.23828125, -569.73828125, -941.0]\n",
+      "[[1.5234375, 0.0, 0.0], [0.0, 1.5234375, 0.0], [0.0, 0.0, 2.0]]\n",
+      "sliceThickness=3.0 +- 0.0\n",
+      "0\n",
+      "[-403.656, -588.125, -942.0]\n",
+      "[[4.07283, 0.0, 0.0], [0.0, 4.07283, 0.0], [0.0, 0.0, 3.0]]\n"
      ]
     }
    ],
@@ -125,36 +135,37 @@
     "    sd=odb.getData('series',seriesId)\n",
     "    firstId=sd['Instances'][0] \n",
     "    \n",
-    "    #check to see if orientation changes (it doesn't)\n",
-    "    #ior=[odb.getDicomField(y,'0020-0037') for y in sd[x]['Instances']] \n",
-    "    #same=all(y==ior[0] for y in ior)\n",
-    "    #print(same)\n",
-    "    \n",
-    "    #check to see if pixel spacing changes (it doesn't)\n",
-    "    #ips=[odb.getDicomField(y,'0028-0030') for y in sd['Instances']] \n",
-    "    #same=all(y==ips[0] for y in ips)\n",
-    "    #print(same)\n",
-    "    \n",
-    "    #check to see if sliceThickness changes (it doesn't)\n",
-    "    #ist=[odb.getDicomField(y,'0018-0050') for y in sd['Instances']]\n",
-    "    #same=all(y==ist[0] for y in ist)\n",
-    "    #print(same)\n",
-    "        \n",
     "    ior=parseOrientation(odb.getDicomField(firstId,'0020-0037'))\n",
+    "    #spacing\n",
     "    ips=parsePos(odb.getDicomField(firstId,'0028-0030'))\n",
     "    ist=parsePos(odb.getDicomField(firstId,'0018-0050')) \n",
-    "    ips.append(ist[0])\n",
-    "        \n",
-    "    for y in range(3):\n",
-    "        ior[y]=[ips[y]*w for w in ior[y]]\n",
-    "    \n",
+    "    #don't use slice thickness for pixel spacing\n",
+    "    #ips.append(ist[0])\n",
     "    \n",
     "    ipos=[parsePos(odb.getDicomField(y,'0020-0032')) for y in sd['Instances']] \n",
-    "    m=[[numpy.dot(numpy.array(v),numpy.array(x)) for x in ipos] for v in ior]\n",
-    "    m=[numpy.min(numpy.array(x)) for x in m]\n",
-    "    print(m)    \n",
+    "    #guess spacing from ipos\n",
+    "    m=[numpy.dot(numpy.array(ior[2]),numpy.array(x)) for x in ipos] \n",
+    "    #sort array\n",
+    "    m1=numpy.sort(m)\n",
+    "    dm=m1[1:]-m1[:-1]\n",
+    "    zSpc=numpy.average(dm)\n",
+    "    zDev=numpy.std(dm)\n",
+    "    print(f'sliceThickness={zSpc} +- {zDev}')\n",
+    "    ips.append(zSpc)\n",
+    "    \n",
+    "    #direction, ie. scaled orientation\n",
+    "    #copy orientation to direction\n",
+    "    xdir=[y for y in ior]\n",
+    "    for i in range(3):\n",
+    "        xdir[i]=[ips[i]*w for w in xdir[i]]\n",
     "    \n",
-    "    return {f'{code}origin':m,f'{code}orientation':ior}\n",
+    "    #find minimum from original array to idenfify lower slice\n",
+    "    im=numpy.argmin(numpy.array(m))\n",
+    "    print(im)\n",
+    "    org=ipos[im]\n",
+    "    print(org)\n",
+    "    print(xdir)\n",
+    "    return {f'{code}origin':org,f'{code}orientation':xdir}\n",
     "    \n",
     "        \n",
     "        \n",
@@ -162,7 +173,8 @@
     "    project='/iPNUMMretro/Study'\n",
     "    schema='study'\n",
     "    query='Imaging1'\n",
-    "    coreDir=os.path.join(os.path.expanduser('~'),'temp')\n",
+    "    #coreDir=os.path.join(os.path.expanduser('~'),'temp')\n",
+    "    coreDir=os.path.join('Z:')\n",
     "    ds=db.selectRows(project,schema,query,[])\n",
     "    rows=ds['rows'][0:1]\n",
     "    imgs={'CT':'CT_orthancId','PET':'PETWB_orthancId'}\n",
@@ -176,8 +188,13 @@
     "            \n",
     "        fname=os.path.join(coreDir,f'{code}.npz')\n",
     "        arr={}\n",
-    "        arr.update({x:numpy.load(odb.getNumpy('series',r[imgs[x]])) for x in imgs})\n",
+    "        #debug  don't load array\n",
+    "        iSUV={x:getSUVfactor(odb,r[imgs[x]]) for x in imgs}\n",
+    "        arr.update({x:iSUV[x]*numpy.load(odb.getNumpy('series',r[imgs[x]])) for x in imgs})\n",
+    "        #debug don't calculate geometry\n",
     "        igeo={x:getGeometry(x,odb,r[imgs[x]]) for x in imgs}\n",
+    "        #debug skip saving\n",
+    "        #continue\n",
     "        _={arr.update(igeo[x]) for x in igeo}\n",
     "        #print(arr)\n",
     "        \n",
@@ -186,15 +203,105 @@
     "        if doUpdate:\n",
     "            db.modifyRows('update',project,schema,query,[r])\n",
     "        \n",
+    "def mergeDateAndTime(date,time):\n",
+    "    time=time.replace(' ','')\n",
+    "    dateTimeStr=' '.join([date,time])\n",
+    "    #print(f'[{date}] [{time}] [{dateTimeStr}]')\n",
+    "    dt=datetime.datetime.strptime(dateTimeStr,'%Y%m%d %H%M%S.%f')\n",
+    "    #print(dt)\n",
+    "    return dt\n",
+    "        \n",
+    "def getSUVfactor(odb,seriesId):\n",
+    "    sd=odb.getData('series',seriesId)\n",
+    "    orthancId=sd['Instances'][0] \n",
+    "    \n",
+    "    MODALITY='0008-0060'\n",
+    "    PATIENT_WEIGHT='0010-1030'\n",
+    "    RADIOPHARMACEUTICAL_STARTTIME='0054-0016/0/0018-1072'\n",
+    "    RADIONUCLIDE_TOTALDOSE='0054-0016/0/0018-1074'\n",
+    "    RADIONUCLIDE_HALFLIFE='0054-0016/0/0018-1075'\n",
+    "    RADIOPHARMACEUTICAL_STARTDATETIME='0054-0016/0/0018-1078'\n",
+    "    SERIES_TIME='0008-0031'\n",
+    "    SERIES_DATE='0008-0021'\n",
+    "    modality=odb.getDicomField(orthancId,MODALITY)\n",
+    "    if modality!='PT':\n",
+    "        print('Returning 1')\n",
+    "        return 1\n",
+    "    weight = float(odb.getDicomField(orthancId,PATIENT_WEIGHT))\n",
+    "    #print(weight)\n",
+    "    halfLife = float(odb.getDicomField(orthancId,RADIONUCLIDE_HALFLIFE))#seconds\n",
+    "    #print(halfLife)\n",
+    "    injDose = float(odb.getDicomField(orthancId,RADIONUCLIDE_TOTALDOSE))#in Bq\n",
+    "    #print(injDose)\n",
+    "    #two options for time of activity measurement (full date time or time only, then match it with studyDate)\n",
+    "    #injDateTime = odb.getDicomField(orthancId,RADIOPHARMACEUTICAL_STARTDATETIME)\n",
+    "    #if troubles:\n",
+    "    injDateTime = mergeDateAndTime(odb.getDicomField(orthancId,SERIES_DATE),\n",
+    "            odb.getDicomField(orthancId,RADIOPHARMACEUTICAL_STARTTIME))\n",
+    "    scanDateTime=mergeDateAndTime(odb.getDicomField(orthancId,SERIES_DATE),\n",
+    "                                  odb.getDicomField(orthancId,SERIES_TIME))\n",
+    "    t=(scanDateTime-injDateTime).total_seconds()\n",
+    "    scanDose = injDose*numpy.exp(-(numpy.log(2)*t)/halfLife);\n",
+    "    SUV_factor = 1/((scanDose/weight)*0.001);\n",
+    "    print(SUV_factor)\n",
+    "    return SUV_factor\n",
+    "\n",
+    "    \n",
     "getFile(db,odb)\n",
     "    "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 36,
    "id": "05cb025f-1c44-4ddf-b859-8b8e9d2584a6",
    "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[1.5234375, 1.5234375, 2.0]\n",
+      "[1. 0. 0. 0. 1. 0. 0. 0. 1.]\n",
+      "[4.07283, 4.07283, 3.0]\n",
+      "[1. 0. 0. 0. 1. 0. 0. 0. 1.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "def convertToITK(dt,label):\n",
+    "    imgArray=dt[label]\n",
+    "    org=dt['{}origin'.format(label)]\n",
+    "    orient=dt['{}orientation'.format(label)]\n",
+    "    spacing=[numpy.linalg.norm(x) for x in orient]\n",
+    "    print(spacing)\n",
+    "    orient=[orient[i]/spacing[i] for i in range(3)]\n",
+    "    direction=numpy.reshape(orient,9)\n",
+    "    print(direction)\n",
+    "    img = SimpleITK.GetImageFromArray(imgArray)\n",
+    "    img.SetOrigin(org)\n",
+    "    img.SetDirection(direction)\n",
+    "    img.SetSpacing(spacing)\n",
+    "    return img\n",
+    "\n",
+    "#convert to SimpleITKimage to display it\n",
+    "infile=os.path.join('Z:','6ebf549213805ee0.npz')\n",
+    "dt=numpy.load(infile)\n",
+    "ct=convertToITK(dt,'CT')\n",
+    "pet=convertToITK(dt,'PET')\n",
+    "outdir=os.path.join('Z:')\n",
+    "SimpleITK.WriteImage(ct,os.path.join(outdir,'CT.nii.gz'))\n",
+    "SimpleITK.WriteImage(pet,os.path.join(outdir,'PET.nii.gz'))\n",
+    "del ct\n",
+    "del pet\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6d523d18-ede0-41aa-a93c-96123fa6abbc",
+   "metadata": {},
    "outputs": [],
    "source": []
   }