Add initial files.

image/downsample_tiffs.py

@@ -0,0 +1,75 @@
+"""
+Downsample the tiffs for easier postprocessing
+"""
+
+import os
+import logging
+import math
+from argparse import ArgumentParser
+from functools import partial
+import multiprocessing as mp
+import skimage.external.tifffile
+import skimage.transform
+from skimage import exposure
+import utility
+
+def do_downsample(input_file, working_directory, downsample_ratios):
+    
+    im = skimage.external.tifffile.imread(input_file)
+    downsampled = skimage.transform.downscale_local_mean(im, downsample_ratios)
+    out_path = working_directory + "/" + os.path.splitext(os.path.basename(input_file))[0] + "_downsampled.tiff"
+    scaled_image = exposure.rescale_intensity(downsampled, out_range=(0.0, 1.0))
+    skimage.external.tifffile.imsave(out_path, skimage.img_as_ubyte(scaled_image))
+    
+def downsample_tiffs(input_folder, downsample_ratios, channel = 0, series = 0):
+    
+    file_names = []
+    dir_name = os.path.dirname(input_folder)
+    for eachFile in os.listdir(input_folder):
+        if eachFile.endswith(".tiff"):
+            file_names.append(dir_name + "/" + eachFile)
+    
+    # Make a partial function for multiple arg mapping
+    working_directory = os.getcwd() + "/Downsampled_X"+str(downsample_ratios[2]) + "_Y"+str(downsample_ratios[1])
+    working_directory += "_Z"+str(downsample_ratios[0]) + "/"
+    if not os.path.exists(working_directory):
+        os.makedirs(working_directory)
+    partial_extract_pixel_values_only = partial(do_downsample, 
+                                                working_directory=working_directory, 
+                                                downsample_ratios=downsample_ratios)
+     
+    # Limit to max 6 cpus
+    num_available_cpus = mp.cpu_count()
+    if num_available_cpus > 6:
+        num_available_cpus = 6
+        
+    # Distribute across processes - Be conservative with available CPUs
+    run_cpus = int(math.floor(num_available_cpus/2))
+    if run_cpus == 0:
+        run_cpus = 1
+    #pool = mp.Pool(processes = run_cpus)
+    pool = mp.Pool(processes = 3)
+    pool.map(partial_extract_pixel_values_only, file_names)
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "zeiss_to_tiff"
+    utility.do_setup(tool_name)
+    logger1 = logging.getLogger('downsample_tiffs.'+tool_name)
+    
+    # Suppress XML Parse warnings
+    pyxb_logger = logging.getLogger('pyxb')
+    pyxb_logger.setLevel(logging.CRITICAL)
+    
+    parser = ArgumentParser()
+    parser.add_argument("-i", "--input_folder", type=str, help='Folder with tiffs to be downsampled.')
+    parser.add_argument("-x", "--x_reduction", type=int, help='Factor reduction in x.', default=1)
+    parser.add_argument("-y", "--y_reduction", type=int, help='Factor reduction in y.', default=1)
+    parser.add_argument("-z", "--z_reduction", type=int, help='Factor reduction in z.', default=1)
+    args = parser.parse_args()
+    
+    logger1.info('Start Downsampling Images at: ' + args.input_folder)
+    downsample_ratios = (args.z_reduction, args.y_reduction, args.x_reduction)
+    downsample_tiffs(args.input_folder, downsample_ratios)
+    logger1.info('Completed Downsampling')

image/processing.py

@@ -0,0 +1,135 @@
+import itk
+
+def default_image_type():
+    return itk.Image[itk.UC, 3]
+
+def image_type_2d():
+    return itk.Image[itk.UC, 2]
+
+def read_image(path, image_type):
+    
+    """
+    Read the image from file
+    """
+    
+    reader = itk.ImageFileReader[image_type].New()
+    reader.SetFileName(path)
+    reader.Update()
+    return reader
+
+def write_image(path, image_type, image_container):
+    
+    """
+    Write the image to file
+    """
+    
+    writer = itk.ImageFileWriter[image_type].New()
+    writer.SetInput(image_container.GetOutput())
+    writer.SetFileName(path)
+    writer.Update()
+    
+def convert_vtk(image_container, image_type):
+    
+    itk_vtk_converter = itk.ImageToVTKImageFilter[image_type].New()
+    itk_vtk_converter.SetInput(image_container.GetOutput())
+    itk_vtk_converter.Update()
+    return itk_vtk_converter
+
+def merge_tiles(image_containers, image_type, num_x):
+    
+    tiler = itk.TileImageFilter[image_type,image_type].New()
+    layout = [num_x, num_x, 0]
+    tiler.SetLayout(layout)
+    for idx, eachContainer in enumerate(image_containers):
+        tiler.SetInput(idx, eachContainer.GetOutput())
+        
+    tiler.Update()
+    return tiler
+
+def median_filter(image_container, image_type, radius):
+    
+    median_filter = itk.MedianImageFilter[image_type, image_type].New()
+    median_filter.SetInput(image_container.GetOutput())
+    median_filter.SetRadius(radius)
+    median_filter.Update()
+    
+    return median_filter
+    
+def correct_spacing(image_container, image_type, z_spacing):
+    
+    correct = itk.ChangeInformationImageFilter[image_type].New()
+    correct.SetInput(image_container.GetOutput())
+    spacing = image_container.GetOutput().GetSpacing()
+    spacing[2] = z_spacing
+    correct.SetOutputSpacing(spacing)
+    correct.ChangeSpacingOn()
+    correct.Update()
+    
+    return correct
+
+def sum_projection(image_container, image_type):
+    output_type = itk.Image[itk.UC, 2]
+    
+    sum_projection = itk.SumProjectionImageFilter[image_type, output_type].New()
+    sum_projection.SetInput(image_container.GetOutput())
+    
+    sum_projection.Update()
+    return sum_projection
+
+def shrink_image(image_type, image_container, factors = [2,2,2]):
+    
+    """
+    Reduce the size of the image in each dimension by the specified factors
+    """
+    
+    shrink = itk.ShrinkImageFilter[image_type, image_type].New()
+    shrink.SetInput(image_container.GetOutput())
+    shrink.SetShrinkFactor(0, factors[0])
+    shrink.SetShrinkFactor(1, factors[1])
+    shrink.SetShrinkFactor(2, factors[2])
+    shrink.Update()
+    return shrink
+
+def get_roi(image_type, image_container, start = [0,0,0], size = [200, 200, 10]):
+    
+    """
+    Select a region of interest in the image
+    """
+    
+    roi = itk.RegionOfInterestImageFilter[image_type, image_type].New()
+    roi.SetInput(image_container.GetOutput())
+
+    region = itk.ImageRegion[3]()
+    region.SetIndex(start)
+    region.SetSize(size)
+
+    roi.SetRegionOfInterest(region)
+    roi.Update()
+    return roi
+
+def threshold_image(image_type, image_container, lower = 195, upper = 200):
+    
+    """
+    Do a binary thresholding
+    """
+    
+    threshold = itk.BinaryThresholdImageFilter[image_type,image_type].New()
+    threshold.SetInput(image_container.GetOutput())
+    threshold.SetLowerThreshold(lower)
+    threshold.SetUpperThreshold(upper)
+    threshold.Update()
+    return threshold
+    
+def get_intensity_histogram(image_type, image_container):
+    
+    """
+    Write an intensity histogram
+    """
+    
+    histogram = itk.ImageToHistogramFilter[image_type].New()
+    histogram.SetInput(image_container.GetOutput())
+    histogram.Update()
+    
+    for idx in range(histogram.GetOutput().GetSize()[0]):
+        print "Hist Intensity: ", idx, " Frequency: ", histogram.GetOutput().GetFrequency(idx)
+    

image/processing/median.py

@@ -0,0 +1,30 @@
+"""
+Median filter on chosen image
+"""
+
+import logging
+from code.image.tools import it_wrapper
+import code.image.settings
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "median"
+    work_dir, image_data = code.image.settings.do_setup(tool_name)
+    logger1 = logging.getLogger('processing.'+tool_name)
+
+    output_path = work_dir+"/median"
+    
+    logger1.info('Start Filtering Image at: ' + work_dir)
+    image_type = it_wrapper.default_image_type()
+    
+    input_path = work_dir + "/convert/itk_tile_0_2_6.tiff"
+    logger1.info("Reading the image: " +input_path + " into ITK")
+    reader = it_wrapper.read_image(input_path, image_type)
+    correct = it_wrapper.correct_spacing(reader, image_type, image_data["VoxelSizeZ"])
+            
+    logger1.info('Doing Median Filter')
+    median_filter = it_wrapper.median_filter(correct, image_type, 1)
+    
+    logger1.info('Writing TIFF Stack')
+    it_wrapper.write_image(output_path+"/itk_tile_0_2_6.tiff", image_type, median_filter)

image/processing/projection.py

@@ -0,0 +1,36 @@
+"""
+SUM projection of tiled images
+"""
+
+import logging
+from code.image.tools import it_wrapper
+import code.image.settings
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "projection"
+    work_dir, image_data = code.image.settings.do_setup(tool_name)
+    logger1 = logging.getLogger('processing.'+tool_name)
+
+    output_path = work_dir+"/projection"
+    
+    logger1.info('Start Projecting Image at: ' + work_dir)
+    image_type = it_wrapper.default_image_type()
+    
+    num_tiles_y = 3
+    num_tiles_x = 3
+    for jdx in range(num_tiles_y):
+        for idx in range(num_tiles_x):
+            label = str(idx) + "_" + str(jdx) +"_"+ str(idx+ num_tiles_x*jdx)
+            input_path = work_dir + "/convert/itk_tile_" + label + ".tiff"
+            logger1.info("Reading the image: " +input_path + " into ITK")
+            reader = it_wrapper.read_image(input_path, image_type)
+            correct = it_wrapper.correct_spacing(reader, image_type, image_data["VoxelSizeZ"])
+            
+            logger1.info('Doing SUM projection')
+            project = it_wrapper.sum_projection(correct, image_type)
+    
+            output_image_type = it_wrapper.image_type_2d()
+            logger1.info('Writing TIFF Stack')
+            it_wrapper.write_image(output_path+"/sum_tile_" + label + ".tiff", output_image_type, project)

image/processing/reduce.py

@@ -0,0 +1,30 @@
+"""
+Median filter on chosen image
+"""
+
+import logging
+from code.image.tools import it_wrapper
+import code.image.settings
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "reduce"
+    work_dir, image_data = code.image.settings.do_setup(tool_name)
+    logger1 = logging.getLogger('processing.'+tool_name)
+
+    output_path = work_dir+"/reduce"
+    
+    logger1.info('Start Reducing Image at: ' + work_dir)
+    image_type = it_wrapper.default_image_type()
+    
+    input_path = work_dir + "/median/itk_tile_0_2_6.tiff"
+    logger1.info("Reading the image: " +input_path + " into ITK")
+    reader = it_wrapper.read_image(input_path, image_type)
+    correct = it_wrapper.correct_spacing(reader, image_type, image_data["VoxelSizeZ"])
+            
+    logger1.info('Doing Reduction')
+    im_reduce = it_wrapper.shrink_image(image_type, correct, [2,2,1])
+    
+    logger1.info('Writing TIFF Stack')
+    it_wrapper.write_image(output_path+"/itk_tile_0_2_6.tiff", image_type, im_reduce)

image/processing/show.py

@@ -0,0 +1,17 @@
+"""
+Show a tiff using vtk
+"""
+
+import code.image.settings
+from code.image.tools import simple_render
+from code.image.tools import it_wrapper
+
+work_dir, input_data = code.image.settings.do_setup("show")
+input_path = work_dir + "/threshold/itk_tile_0_2_6.tiff"
+
+image_type = it_wrapper.default_image_type()
+reader = it_wrapper.read_image(input_path, image_type)
+correct = it_wrapper.correct_spacing(reader, image_type, input_data["VoxelSizeZ"])
+converter = it_wrapper.convert_vtk(reader, image_type)
+    
+simple_render.render_image(converter, threshold = 1.0)

image/processing/threshold.py

@@ -0,0 +1,30 @@
+"""
+Threshold on chosen image
+"""
+
+import logging
+from code.image.tools import it_wrapper
+import code.image.settings
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "threshold"
+    work_dir, image_data = code.image.settings.do_setup(tool_name)
+    logger1 = logging.getLogger('processing.'+tool_name)
+
+    output_path = work_dir+"/threshold"
+    
+    logger1.info('Start Thresholding Image at: ' + work_dir)
+    image_type = it_wrapper.default_image_type()
+    
+    input_path = work_dir + "/reduce/itk_tile_0_2_6.tiff"
+    logger1.info("Reading the image: " +input_path + " into ITK")
+    reader = it_wrapper.read_image(input_path, image_type)
+    correct = it_wrapper.correct_spacing(reader, image_type, image_data["VoxelSizeZ"])
+            
+    logger1.info('Doing Thresholding')
+    threshold = it_wrapper.threshold_image(image_type, correct, lower = 20, upper = 255)
+    
+    logger1.info('Writing TIFF Stack')
+    it_wrapper.write_image(output_path+"/itk_tile_0_2_6.tiff", image_type, threshold)

image/processing/tile.py

@@ -0,0 +1,39 @@
+"""
+Merging of tiled images
+"""
+
+import logging
+from code.image.tools import it_wrapper
+import code.image.settings
+
+if __name__ == "__main__":
+    
+    # Do setup
+    tool_name = "tile"
+    work_dir, image_data = code.image.settings.do_setup(tool_name)
+    logger1 = logging.getLogger('processing.'+tool_name)
+
+    output_path = work_dir+"/projection"
+    
+    logger1.info('Start Merging Images at: ' + work_dir)
+    image_type = it_wrapper.default_image_type()
+    
+    num_tiles_y = 3
+    num_tiles_x = 3
+    container = []
+    for jdx in range(num_tiles_y):
+        for idx in range(num_tiles_x):
+            
+            label = str(idx) + "_" + str(jdx) +"_"+ str(idx+ num_tiles_x*jdx)
+            input_path = work_dir + "/projection/sum_tile_" + label + ".tiff"
+            logger1.info("Reading the image: " +input_path + " into ITK")
+            reader = it_wrapper.read_image(input_path, image_type)
+            correct = it_wrapper.correct_spacing(reader, image_type, image_data["VoxelSizeZ"])
+            container.append(correct)
+    
+    logger1.info('Merging Tiles')    
+    merged = it_wrapper.merge_tiles(container, image_type, num_tiles_x)
+      
+    output_image_type = it_wrapper.default_image_type()
+    logger1.info('Writing TIFF Stack')
+    it_wrapper.write_image(output_path+"/sum_merge.tiff", output_image_type, merged)

image/surface.py

@@ -0,0 +1,171 @@
+import vtk
+from vmtk import pypes
+
+class Surface2d():
+    
+    def __init__(self, input_path):
+        self.input_path = input_path
+        self.surface = None
+        self.boundaries = None
+    
+    def generate(self, target_reduction = 0.9994, num_subdivisions = 3, resize_factors = None, clipping_factor = 0.0):
+        
+        # Convert to VTI
+        myArguments = 'vmtkimagereader -ifile ' + self.input_path
+        my_pype = pypes.PypeRun(myArguments)
+        
+        vtk_image = my_pype.GetScriptObject('vmtkimagereader','0').Image
+        
+        # Size the image if needed
+        if resize_factors is not None:
+            reduction_filter = vtk.vtkImageResize()
+            reduction_filter.SetInput(vtk_image)
+            reduction_filter.SetMagnificationFactors(resize_factors[0], resize_factors[1], resize_factors[2])
+            reduction_filter.SetResizeMethodToMagnificationFactors()
+            reduction_filter.Update()
+        
+        # Threshold
+        threshold = vtk.vtkThreshold()
+        if resize_factors is not None:
+            threshold.SetInputConnection(reduction_filter.GetOutputPort())
+        else:
+            threshold.SetInput(vtk_image)
+        threshold.ThresholdByUpper(1.0)
+        threshold.Update()
+        
+        # Convert to polydata
+        surface = vtk.vtkGeometryFilter()
+        surface.SetInputConnection(threshold.GetOutputPort())
+        surface.Update()
+        
+        # Triangulate
+        triangle = vtk.vtkTriangleFilter()
+        triangle.SetInputConnection(surface.GetOutputPort())
+        triangle.Update()
+    
+        # Decimate
+        decimate = vtk.vtkDecimatePro()
+        decimate.SetInputConnection(triangle.GetOutputPort())
+        decimate.SetTargetReduction(target_reduction)
+        decimate.SetFeatureAngle(15.0)
+        decimate.Update()
+        
+        # Do loop subdivision
+        su = vtk.vtkLinearSubdivisionFilter()
+        su.SetInputConnection(decimate.GetOutputPort())
+        su.SetNumberOfSubdivisions(num_subdivisions)
+        su.Update()
+        
+        # Clip the boundaries, recommended to ensure straight inlets and outlets
+        if clipping_factor > 0.0:
+            bounds = su.GetOutput().GetBounds()
+            
+            width = bounds[1] - bounds[0]
+            height = bounds[3] - bounds[2]
+            print bounds[2], bounds[3], height
+            
+            p1 = vtk.vtkPlane()
+            p1.SetOrigin(bounds[0] + clipping_factor*width, 0.0, 0)
+            p1.SetNormal(1, 0, 0)
+            
+            p2 = vtk.vtkPlane()
+            p2.SetOrigin(0.0, bounds[2] + clipping_factor*height, 0)
+            p2.SetNormal(0, 1, 0)
+            
+            p3 = vtk.vtkPlane()
+            p3.SetOrigin(bounds[1] - clipping_factor*width, 0.0, 0)
+            p3.SetNormal(-1, 0, 0)
+        
+            p4 = vtk.vtkPlane()
+            p4.SetOrigin(0.0, bounds[3] - clipping_factor*height, 0)
+            p4.SetNormal(0, -1, 0)        
+        
+            c1 = vtk.vtkClipPolyData()
+            c1.SetInputConnection(su.GetOutputPort())
+            c1.SetClipFunction(p1)
+            c1.SetValue(0.0)
+            
+            c2 = vtk.vtkClipPolyData()
+            c2.SetInputConnection(c1.GetOutputPort())
+            c2.SetClipFunction(p2)
+            c2.SetValue(0.0)
+# 
+            c3 = vtk.vtkClipPolyData()
+            c3.SetInputConnection(c2.GetOutputPort())
+            c3.SetClipFunction(p3)
+            c3.SetValue(0.0)
+             
+            c4 = vtk.vtkClipPolyData()
+            c4.SetInputConnection(c3.GetOutputPort())
+            c4.SetClipFunction(p4)
+            c4.SetValue(0.0)
+            
+            su = vtk.vtkGeometryFilter()
+            su.SetInputConnection(c4.GetOutputPort())
+            su.Update()
+        
+        feature_edges = vtk.vtkFeatureEdges()
+        feature_edges.SetInputConnection(su.GetOutputPort())  
+        feature_edges.Update()   
+        
+        clean = vtk.vtkCleanPolyData()
+        clean.SetInputConnection(feature_edges.GetOutputPort())
+        clean.Update()
+         
+        triangle2 = vtk.vtkTriangleFilter()
+        triangle2.SetInputConnection(clean.GetOutputPort())
+        triangle2.Update()
+        
+        self.surface = su.GetOutput()
+        self.boundaries = triangle2.GetOutput()
+        return su.GetOutput(), triangle2.GetOutput()
+    
+    def write(self, output_directory):
+        
+        surface_writer = vtk.vtkXMLPolyDataWriter()
+        surface_writer.SetFileName(output_directory + "/surface.vtp")
+        surface_writer.SetInput(self.surface)
+        surface_writer.Write() 
+        
+        surface_writer.SetFileName(output_directory + "/boundaries.vtp")
+        surface_writer.SetInput(self.boundaries)
+        surface_writer.Write() 
+        
+class Surface3d():
+    
+    def __init__(self):
+        pass
+    
+    def generate(self, input_path, output_directory, part_name):
+        
+        input_path = work_dir + "/vtk_image/itk_tile_0_2_6.vti"
+        
+        reader = vtk.vtkXMLImageDataReader()
+        reader.SetFileName(input_path)
+        reader.Update()
+         
+        pad_filter = vtk.vtkImageConstantPad()
+        pad_filter.SetInputConnection(reader.GetOutputPort())
+        pad_filter.SetConstant(0)
+        pad_filter.SetOutputWholeExtent(reader.GetOutput().GetExtent()[0]-5, 
+                                        reader.GetOutput().GetExtent()[1]+5, 
+                                        reader.GetOutput().GetExtent()[2]-5, 
+                                        reader.GetOutput().GetExtent()[3]+5, 
+                                        reader.GetOutput().GetExtent()[4]-5, 
+                                        reader.GetOutput().GetExtent()[5]+5)
+        pad_filter.Update()
+         
+        writer = vtk.vtkXMLImageDataWriter()
+        writer.SetFileName(work_dir + "/surface/itk_tile_0_2_6_pad.vti")
+        writer.SetInputConnection(pad_filter.GetOutputPort())
+        writer.Update()    
+        
+        
+        myArguments = 'vmtklevelsetsegmentation -ifile ' + work_dir + "/surface/itk_tile_0_2_6_pad.vti" + ' -ofile ' + output_path+"/itk_tile_0_2_6_ls.vti"
+     
+        myPype = pypes.PypeRun(myArguments)
+         
+        myArguments = 'vmtkmarchingcubes -ifile ' + output_path+"/itk_tile_0_2_6_ls.vti" + ' -ofile ' + output_path+"/itk_tile_0_2_6_model.vtp"
+     
+        myPype = pypes.PypeRun(myArguments)
+