Merge pull request dave3d#10 from dave3d/circleci-project-setup

Circleci project setup

dicom2stl.py

@@ -244,30 +244,13 @@ def usage():
 
 from utils import dicomutils
 
-#  Unzip a zipfile of dicom images into a temp directory, then
-#  load the series that has the most slices
-#
-
-
-def loadZipDicom(name):
-    print("Reading Dicom zip file:", name)
-    myzip = zipfile.ZipFile(name, 'r')
-
-    try:
-        myzip.extractall(tempDir)
-    except:
-        print("Zip extract failed")
-
-    return loadLargestSeries(tempDir)
-
-
 #  Load our Dicom data
 #
 if zipFlag:
     # Case for a zip file of images
     if verbose:
         print("zip")
-    img, modality = loadZipDicom(fname[0])
+    img, modality = dicomutils.loadZipDicom(fname[0], tempDir)
 
 
 else:

tests/create_data.py

@@ -0,0 +1,108 @@
+#! /usr/bin/env python
+
+import sys
+import SimpleITK as sitk
+import getopt
+
+def make_tetra(dim=128, pixel_type=sitk.sitkUInt8):
+    dim4 = dim>>2
+    dim2 = dim>>1
+    h = .86602540378443864676*dim2+dim4
+    size =[ dim,dim,dim ]
+    p1 = [dim4,dim4,dim4]
+    p2 = [dim2+dim4,dim4,dim4]
+    p3 = [dim2, h,dim4]
+    p4 = [dim2,dim2,h]
+    points = [p1,p2,p3,p4]
+
+    sigma=[dim/6,dim/6,dim/6]
+    v1 = sitk.GaussianSource(pixel_type, size, sigma=sigma, mean=p1, scale=200)
+    v2 = sitk.GaussianSource(pixel_type, size, sigma=sigma, mean=p2, scale=200)
+    v3 = sitk.GaussianSource(pixel_type, size, sigma=sigma, mean=p3, scale=200)
+    v4 = sitk.GaussianSource(pixel_type, size, sigma=sigma, mean=p4, scale=200)
+
+    vsum = v1+v2+v3+v4
+
+    return vsum
+
+def make_cylinder(dim=64, pixel_type=sitk.sitkUInt8):
+    mean=[dim/2,dim/2]
+    sigma=[dim/4,dim/4]
+    img = sitk.GaussianSource(pixel_type, [dim,dim], sigma=sigma, mean=mean, scale=200)
+
+    series = []
+    for i in range(dim):
+        series.append(img)
+
+    vol = sitk.JoinSeries(series)
+    return vol
+
+def vol2dicom(vol, root_name):
+    """ this kinda sucks.  doesn't do keep the UID tag properly. """
+    for z in range(vol.GetDepth()):
+       img = vol[:,:,z]
+       img.SetMetaData("0020|000e", "1.2.3.4.5.6.7.8")
+       name = "%s.%d.dcm" % (root_name, z)
+       print(name)
+       sitk.WriteImage(img, name)
+
+
+def usage():
+    print()
+    print(" create_data.py [options] output_image")
+    print()
+    print(" -h, --help          This message")
+    print(" -d int, --dim int   Output image dimensions (default=32)")
+    print(" -c , --cylinder     Cylinder volume (default)")
+    print(" -t , --tetra        Tetrahedral volume")
+    print(" -p type , --pixel type        Pixel type by name (default=UInt8)")
+    print()
+
+if __name__ == "__main__":
+
+    dim = 32
+    vtype = "cylinder"
+    ptype = sitk.sitkUInt8
+
+    try:
+        opts, args = getopt.getopt(sys.argv[1:], "hd:ctp:",
+            [ "help", "dim", "cylinder", "tetra", "pixel=" ] )
+    except getopt.GetoptError as err:
+        print (str(err))
+        usage()
+        sys.exit(2)
+
+
+    for o, a in opts:
+        if o in ("-h", "--help"):
+            usage()
+            sys.exit()
+        elif o in ("-d", "--dim"):
+            dim = int(a)
+        elif o in ("-c", "--cylinder"):
+            vtype = "cylinder"
+        elif o in ("-t", "--tetra"):
+            vtype = "tetra"
+        elif o in ("-p", "--pixel"):
+            if a.lower() == "uint16":
+                ptype = sitk.sitkUInt16
+            elif a.lower() == "int16":
+                ptype = sitk.sitkInt16
+            elif a.lower() == "int32":
+                ptype = sitk.sitkInt32
+            elif a.lower() == "float32":
+                ptype = sitk.sitkFloat32
+            elif a.lower() == "float64":
+                ptype = sitk.sitkFloat64
+        else:
+            assert False, "unhandled options"
+
+    if vtype == "tetra":
+        print("Making tetra")
+        vol = make_tetra(dim, ptype)
+    else:
+        print("Making cylinder")
+        vol = make_cylinder(dim, ptype)
+
+    print("Writing", args[0])
+    sitk.WriteImage(vol, args[0])

tests/test_dicomutils.py

@@ -0,0 +1,80 @@
+#! /usr/bin/env python
+
+import unittest
+import os, shutil, zipfile
+import vtk
+import SimpleITK as sitk
+
+from utils import dicomutils
+from tests import create_data
+from tests import write_series
+
+class TestDicomUtils(unittest.TestCase):
+
+    TMPDIR = "testtmp"
+    SIZE = 32
+
+    @classmethod
+    def setUpClass(cls):
+        print("\nBuildin\' it up!")
+        cyl = create_data.make_cylinder(TestDicomUtils.SIZE, sitk.sitkUInt16)
+        try:
+            os.mkdir(TestDicomUtils.TMPDIR)
+        except:
+            print("Oopsie")
+        write_series.write_series( cyl, TestDicomUtils.TMPDIR )
+
+    @classmethod
+    def tearDownClass(cls):
+        print("\nTearin\' it down!")
+        shutil.rmtree(TestDicomUtils.TMPDIR)
+
+    def test_scanDirForDicom(self):
+        print("\nTesting DicomUtils.scanDirForDicom")
+        matches, dirs = dicomutils.scanDirForDicom(TestDicomUtils.TMPDIR)
+        print(matches, dirs)
+        self.assertEqual(len(matches), TestDicomUtils.SIZE)
+
+    def test_getAllSeries(self):
+        print("\nTesting DicomUtils.getAllSeries")
+        seriessets = dicomutils.getAllSeries([TestDicomUtils.TMPDIR])
+        print(seriessets)
+        self.assertEqual(len(seriessets),1)
+        series_id = seriessets[0][0]
+        if series_id.startswith('1.2.826.0.1.3680043'):
+           print("    Series looks good")
+        else:
+           self.fail("    Bad series: "+series_id)
+
+    def test_getModality(self):
+        print("\nTesting DicomUtils.getModality")
+        img = sitk.Image(10,10,sitk.sitkUInt16)
+        m1 = dicomutils.getModality(img)
+        self.assertEqual(m1, "")
+        img.SetMetaData("0008|0060", "dude")
+        m2 = dicomutils.getModality(img)
+        self.assertEqual(m2, "dude")
+
+    def test_loadLargestSeries(self):
+        print("\nTesting DicomUtils.loadLargestSeries")
+        img, mod = dicomutils.loadLargestSeries(TestDicomUtils.TMPDIR)
+        self.assertEqual(img.GetSize(), (TestDicomUtils.SIZE,TestDicomUtils.SIZE,TestDicomUtils.SIZE))
+        self.assertEqual(mod, "CT")
+
+    def test_loadZipDicom(self):
+        print("\nTesting DicomUtils.loadZipDicom")
+        zf = zipfile.ZipFile('tests/testzip.zip', 'w')
+        for z in range(TestDicomUtils.SIZE):
+            zf.write(TestDicomUtils.TMPDIR+'/'+str(z)+'.dcm')
+        zf.close()
+        img, mod = dicomutils.loadZipDicom('tests/testzip.zip', 'tests/ziptmp')
+        print(img.GetSize())
+        print(mod)
+        os.unlink('tests/testzip.zip')
+        shutil.rmtree('tests/ziptmp')
+
+
+
+if __name__ == "__main__":
+    unittest.main()
+

tests/write_series.py

@@ -0,0 +1,146 @@
+#! /usr/bin/env python
+
+import SimpleITK as sitk
+
+import sys, time, os
+import numpy as np
+
+pixel_dtypes = {"int16" : np.int16,
+                "float64" : np.float64}
+
+
+def writeSlices(series_tag_values, new_img, out_dir, writer, i):
+    image_slice = new_img[:,:,i]
+
+    # Tags shared by the series.
+    list(map(lambda tag_value: image_slice.SetMetaData(tag_value[0], tag_value[1]), series_tag_values))
+
+    # Slice specific tags.
+    image_slice.SetMetaData("0008|0012", time.strftime("%Y%m%d")) # Instance Creation Date
+    image_slice.SetMetaData("0008|0013", time.strftime("%H%M%S")) # Instance Creation Time
+
+    # Setting the type to CT preserves the slice location.
+    image_slice.SetMetaData("0008|0060", "CT")  # set the type to CT so the thickness is carried over
+
+    # (0020, 0032) image position patient determines the 3D spacing between slices.
+    image_slice.SetMetaData("0020|0032", '\\'.join(map(str,new_img.TransformIndexToPhysicalPoint((0,0,i))))) # Image Position (Patient)
+    image_slice.SetMetaData("0020,0013", str(i)) # Instance Number
+
+    # Write to the output directory and add the extension dcm, to force writing in DICOM format.
+    writer.SetFileName(os.path.join(out_dir, str(i)+'.dcm'))
+    writer.Execute(image_slice)
+
+
+
+# Write the 3D image as a series
+# IMPORTANT: There are many DICOM tags that need to be updated when you modify an
+#            original image. This is a delicate opration and requires knowlege of
+#            the DICOM standard. This example only modifies some. For a more complete
+#            list of tags that need to be modified see:
+#                           http://gdcm.sourceforge.net/wiki/index.php/Writing_DICOM
+#            If it is critical for your work to generate valid DICOM files,
+#            It is recommended to use David Clunie's Dicom3tools to validate the files
+#                           (http://www.dclunie.com/dicom3tools.html).
+
+def write_series(new_img, data_directory, pixel_dtype=np.int16):
+    writer = sitk.ImageFileWriter()
+    # Use the study/series/frame of reference information given in the meta-data
+    # dictionary and not the automatically generated information from the file IO
+    writer.KeepOriginalImageUIDOn()
+
+    modification_time = time.strftime("%H%M%S")
+    modification_date = time.strftime("%Y%m%d")
+
+    # Copy some of the tags and add the relevant tags indicating the change.
+    # For the series instance UID (0020|000e), each of the components is a number, cannot start
+    # with zero, and separated by a '.' We create a unique series ID using the date and time.
+    # tags of interest:
+    direction = new_img.GetDirection()
+    series_tag_values = [("0008|0031",modification_time), # Series Time
+                      ("0008|0021",modification_date), # Series Date
+                      ("0008|0008","DERIVED\\SECONDARY"), # Image Type
+                      ("0020|000e", "1.2.826.0.1.3680043.2.1125."+modification_date+".1"+modification_time), # Series Instance UID
+                      ("0020|0037", '\\'.join(map(str, (direction[0], direction[3], direction[6],# Image Orientation (Patient)
+                                                        direction[1],direction[4],direction[7])))),
+                      ("0008|103e", "Created-SimpleITK")] # Series Description
+
+    if pixel_dtype == np.float64:
+        # If we want to write floating point values, we need to use the rescale slope, "0028|1053", to select the
+        # number of digits we want to keep. We also need to specify additional pixel storage and representation
+        # information.
+        rescale_slope = 0.001 #keep three digits after the decimal point
+        series_tag_values = series_tag_values + \
+                            [('0028|1053', str(rescale_slope)), #rescale slope
+                             ('0028|1052','0'),   #rescale intercept
+                             ('0028|0100', '16'), #bits allocated
+                             ('0028|0101', '16'), #bits stored
+                             ('0028|0102', '15'), #high bit
+                             ('0028|0103','1')] #pixel representation
+
+    # Write slices to output directory
+    list(map(lambda i: writeSlices(series_tag_values, new_img, data_directory, writer, i), range(new_img.GetDepth())))
+
+
+def do_test(data_directory):
+    # Re-read the series
+    # Read the original series. First obtain the series file names using the
+    # image series reader.
+    series_IDs = sitk.ImageSeriesReader.GetGDCMSeriesIDs(data_directory)
+    if not series_IDs:
+        print("ERROR: given directory \""+data_directory+"\" does not contain a DICOM series.")
+        sys.exit(1)
+    series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(data_directory, series_IDs[0])
+
+    series_reader = sitk.ImageSeriesReader()
+    series_reader.SetFileNames(series_file_names)
+
+    # Configure the reader to load all of the DICOM tags (public+private):
+    # By default tags are not loaded (saves time).
+    # By default if tags are loaded, the private tags are not loaded.
+    # We explicitly configure the reader to load tags, including the
+    # private ones.
+    series_reader.LoadPrivateTagsOn()
+    image3D = series_reader.Execute()
+    print(image3D.GetSpacing(),'vs',new_img.GetSpacing())
+    sys.exit( 0 )
+
+if __name__ == "__main__":
+
+    if len( sys.argv ) < 3:
+        print( "Usage: python " + __file__ + " <output_directory> [" + ", ".join(pixel_dtypes) + "]" )
+        print( "  or " )
+        print( "       python " + __file__ + " input_volume <output_directory> [" + ", ".join(pixel_dtypes) + "]" )
+        sys.exit ( 1 )
+
+    inname = ""
+
+    if len(sys.argv)>3:
+        inname = sys.argv.pop(1)
+        print("Reading volume", inname)
+
+    # Create a new series from a numpy array
+    try:
+        pixel_dtype = pixel_dtypes[sys.argv[2]]
+    except KeyError:
+        pixel_dtype = pixel_dtypes["int16"]
+
+
+    data_directory = sys.argv[1]
+
+    if len(inname)>0:
+        new_img = sitk.ReadImage(inname)
+        if pixel_dtype == "float64":
+           new_img = sitk.Cast(new_img, sitk.sitkFloat64)
+        else:
+           new_img = sitk.Cast(new_img, sitk.sitkInt16)
+
+    else:
+        new_arr = np.random.uniform(-10, 10, size = (3,4,5)).astype(pixel_dtype)
+        new_img = sitk.GetImageFromArray(new_arr)
+        new_img.SetSpacing([2.5,3.5,4.5])
+
+    write_series(new_img, pixel_dtype, data_directory)
+
+    if len(inname) == 0:
+        do_test(data_directory)
+

utils/dicomutils.py

@@ -95,6 +95,22 @@ def loadLargestSeries(dicomdir):
     return img, modality
 
 
+def loadZipDicom(name, tempDir):
+    """ Unzip a zipfile of dicom images into a temp directory, then
+        load the series that has the most slices.
+    """
+
+    print("Reading Dicom zip file:", name)
+    myzip = zipfile.ZipFile(name, 'r')
+
+    try:
+        myzip.extractall(tempDir)
+    except:
+        print("Zip extract failed")
+
+    return loadLargestSeries(tempDir)
+
+
 #
 #   Main (test code)
 #