mcdlee · October 21, 2021 06:22 · Oct 21, 2021 · Jan 20, 2014 · Jan 8, 2014 · Jan 8, 2014
diff --git a/PET_parsing.py b/PET_parsing.py
@@ -9,7 +9,7 @@
 # collect pixel_volume and resolution
 
 def getSec(s):
-    b =int(s[0:2]) *3600 + int(s[2:4])*60 + int(s[5:6])
+    b =int(s[0:2]) *3600 + int(s[2:4])*60 + int(s[4:6])
     return b
 
 for record in ds.DirectoryRecordSequence:

diff --git a/PET_parsing.py b/PET_parsing.py
@@ -8,6 +8,10 @@
 
 # collect pixel_volume and resolution
 
+def getSec(s):
+    b =int(s[0:2]) *3600 + int(s[2:4])*60 + int(s[5:6])
+    return b
+
 for record in ds.DirectoryRecordSequence:
     if record.DirectoryRecordType == "IMAGE":
         if record.ImageID == "PET AC 2D":
@@ -18,7 +22,13 @@
             row = dcm.Rows
             column = dcm.Columns
             slice_number = dcm.InstanceNumber
+            body_weight = dcm.PatientWeight *1000
+            injection_time = getSec(dcm.RadiopharmaceuticalInformationSequence[0].RadiopharmaceuticalStartTime)
+            acqusition_time = getSec(dcm.AcquisitionTime)
+            dose = dcm.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose
+            half_life = dcm.RadiopharmaceuticalInformationSequence[0].RadionuclideHalfLife
             break
+actural_dose = dose * 0.5**((acqusition_time-injection_time)/half_life)
 
 pixel_volume = pixel_spacing ** 2 * slice_thickness/1000
 print("pixel size =", pixel_spacing, "x", pixel_spacing, "x", slice_thickness, "mm")
@@ -35,9 +45,7 @@
             # Extract the relative path to the DICOM file
             path = os.path.join(*record.ReferencedFileID)
             dcm = dicom.read_file(path)
-            body_weight = dcm.PatientWeight *1000
-            dose = dcm.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose
-            SUV = dcm.pixel_array * body_weight / dose * record[0x0028, 0x1053].value # Rescale Slope
+            SUV = dcm.pixel_array * body_weight / actural_dose * record[0x0028, 0x1053].value # Rescale Slope
             # Now get your image data
             pixel_data[i,:,:] = SUV 
             i = i-1
diff --git a/PET_parsing.py b/PET_parsing.py
@@ -29,7 +29,7 @@
 
 #combine the arrays into 3D dataset
 i = slice_number-1 #Because our scan is from pelvis to head
-for record in ds.DirectoryRecordSequence:
+for record in ds.DirectoryRecordSequence: #https://groups.google.com/forum/#!msg/pydicom/TR2nl6_JRtM/aDODx-eNUTcJ
     if record.DirectoryRecordType == "IMAGE":
         if record.ImageID == "PET AC 2D":
             # Extract the relative path to the DICOM file

diff --git a/PET_parsing.py b/PET_parsing.py
@@ -4,7 +4,7 @@
 import pylab
 
 # readind data
-ds = dicom.read_file("DICOMDIR") # read a dicom of pre-reconstructed SPECT
+ds = dicom.read_file("DICOMDIR") # read a dicomdir of PET/CT
 
 # collect pixel_volume and resolution
 

diff --git a/PET_parsing.py b/PET_parsing.py
@@ -28,7 +28,7 @@
 pixel_data = numpy.ones((slice_number,row,column))
 
 #combine the arrays into 3D dataset
-i = slice_number-1
+i = slice_number-1 #Because our scan is from pelvis to head
 for record in ds.DirectoryRecordSequence:
     if record.DirectoryRecordType == "IMAGE":
         if record.ImageID == "PET AC 2D":

diff --git a/PET.py → PET_parsing.py b/PET.py → PET_parsing.py
@@ -28,7 +28,7 @@
 pixel_data = numpy.ones((slice_number,row,column))
 
 #combine the arrays into 3D dataset
-i = 256
+i = slice_number-1
 for record in ds.DirectoryRecordSequence:
     if record.DirectoryRecordType == "IMAGE":
         if record.ImageID == "PET AC 2D":

diff --git a/PET.py b/PET.py
@@ -6,24 +6,26 @@
 # readind data
 ds = dicom.read_file("DICOMDIR") # read a dicom of pre-reconstructed SPECT
 
-# collect pixel_volume
-i = 0
+# collect pixel_volume and resolution
+
 for record in ds.DirectoryRecordSequence:
     if record.DirectoryRecordType == "IMAGE":
         if record.ImageID == "PET AC 2D":
             path = os.path.join(*record.ReferencedFileID)
             dcm = dicom.read_file(path)
             pixel_spacing = dcm.PixelSpacing[1]
             slice_thickness = record.SliceThickness
-            if i ==0:
-                break
+            row = dcm.Rows
+            column = dcm.Columns
+            slice_number = dcm.InstanceNumber
+            break
 
 pixel_volume = pixel_spacing ** 2 * slice_thickness/1000
 print("pixel size =", pixel_spacing, "x", pixel_spacing, "x", slice_thickness, "mm")
 print("pixel volume=", pixel_volume, "cc")
+print("This image is", row, "x", column, "x", slice_number)
 
-pixel_data = numpy.ones((257,128,128))
-
+pixel_data = numpy.ones((slice_number,row,column))
 
 #combine the arrays into 3D dataset
 i = 256
@@ -33,9 +35,9 @@
             # Extract the relative path to the DICOM file
             path = os.path.join(*record.ReferencedFileID)
             dcm = dicom.read_file(path)
+            body_weight = dcm.PatientWeight *1000
+            dose = dcm.RadiopharmaceuticalInformationSequence[0].RadionuclideTotalDose
+            SUV = dcm.pixel_array * body_weight / dose * record[0x0028, 0x1053].value # Rescale Slope
             # Now get your image data
-            pixel_data[i,:,:] = dcm.pixel_array
-            i = i-1
-
-pylab.imshow(pixel_data[:,60,:])
-pylab.show()
+            pixel_data[i,:,:] = SUV 
+            i = i-1
diff --git a/PET.py b/PET.py
@@ -0,0 +1,41 @@
+import dicom
+import os
+import numpy
+import pylab
+
+# readind data
+ds = dicom.read_file("DICOMDIR") # read a dicom of pre-reconstructed SPECT
+
+# collect pixel_volume
+i = 0
+for record in ds.DirectoryRecordSequence:
+    if record.DirectoryRecordType == "IMAGE":
+        if record.ImageID == "PET AC 2D":
+            path = os.path.join(*record.ReferencedFileID)
+            dcm = dicom.read_file(path)
+            pixel_spacing = dcm.PixelSpacing[1]
+            slice_thickness = record.SliceThickness
+            if i ==0:
+                break
+
+pixel_volume = pixel_spacing ** 2 * slice_thickness/1000
+print("pixel size =", pixel_spacing, "x", pixel_spacing, "x", slice_thickness, "mm")
+print("pixel volume=", pixel_volume, "cc")
+
+pixel_data = numpy.ones((257,128,128))
+
+
+#combine the arrays into 3D dataset
+i = 256
+for record in ds.DirectoryRecordSequence:
+    if record.DirectoryRecordType == "IMAGE":
+        if record.ImageID == "PET AC 2D":
+            # Extract the relative path to the DICOM file
+            path = os.path.join(*record.ReferencedFileID)
+            dcm = dicom.read_file(path)
+            # Now get your image data
+            pixel_data[i,:,:] = dcm.pixel_array
+            i = i-1
+
+pylab.imshow(pixel_data[:,60,:])
+pylab.show()