Thresholded Region Growing (3D Flood Fill)

region_growing.py

import numpy as np

# start_point = (1,2,3)  # Z, X, Y coordinate values

def region_grow(vol, mask, start_point, epsilon=5, HU_mid=0, HU_range=0, fill_with=1):

    sizez = vol.shape[0] - 1
    sizex = vol.shape[1] - 1
    sizey = vol.shape[2] - 1

    items = []
    visited = []

    def enqueue(item):
        items.insert(0,item)

    def dequeue():
        s = items.pop()
        visited.append(s)
        return s

    enqueue((start_point[0], start_point[1], start_point[2]))

    while not items==[]:

        z,x,y = dequeue()
        
        voxel = vol[z,x,y]
        mask[z,x,y] = fill_with
        
        if x<sizex:
            tvoxel = vol[z,x+1,y]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z,x+1,y))

        if x>0:
            tvoxel = vol[z,x-1,y]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z,x-1,y))

        if y<sizey:
            tvoxel = vol[z,x,y+1]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z,x,y+1))

        if y>0:
            tvoxel = vol[z,x,y-1]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z,x,y-1))

        if z<sizez:
            tvoxel = vol[z+1,x,y]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z+1,x,y))

        if z>0:
            tvoxel = vol[z-1,x,y]
            if abs(tvoxel-voxel)<epsilon  and  abs(tvoxel-HU_mid)<HU_range:  enqueue((z-1,x,y))