3d скелет из сегментации

Я хочу создать скелет на основе существующей сегментации, аналогично тому, как это сделано здесь (из sk-образ):

моя цель

Однако я хочу сделать это на 3D-данных. Есть ли код для этого где-нибудь? Желательно на питоне, но помогает любой язык.

Я знаю об этот замечательный сайт, однако я думаю, они не предлагают никакого кода.

Я планирую использовать это на объемах около 500x500x500 пикселей, поэтому он должен хорошо масштабироваться...


person kyra    schedule 07.08.2015    source источник
comment
Там не так много кода для 3D-данных... Я не знаю, найдете ли вы его, но вы можете попробовать посмотреть код scikits-image (у него есть общедоступный репозиторий на github) и попробуйте запрограммировать его 3D-версию . Если вы это сделаете, scikits-image с радостью примет код для расширения библиотеки.   -  person Imanol Luengo    schedule 07.08.2015
comment
Приготовьтесь испачкать руки, это будет непросто! Удачи в кодировании вашего 3D-скелетирования и, пожалуйста, отправьте свой прогресс в scikits-image!!   -  person Ander Biguri    schedule 07.08.2015
comment
FWIW, scikit-image реализовал 3D-скелетонизацию в 2016 году. См., Например, github.com/scikit-image/scikit-image/pull/   -  person Stefan van der Walt    schedule 12.12.2019


Ответы (1)


Я разрабатываю эти инструменты по этой ссылке ниже. Функция getSkeletonize3D в программе с именем convOptimize.py позволяет уменьшить объем ваших 3D-данных. Выдача результата для имеющегося у меня 512 куба заняла около 30 минут. Дайте мне знать, если у вас возникнут проблемы. https://github.com/3Scan/3scan-skeleton. Документ, который я использовал для реализации, находится в комментариях в коде ниже.

По сути, этот алгоритм 3D-скелетирования работает следующим образом: в каждом проходе он имеет 12 подитераций, в которых он итеративно удаляет границы в определенных направлениях, пока вы не получите скелет в центре.

Основной код Python, необходимый для скелетирования ваших данных, приведен ниже. Поскольку для этого требуется импорт из других программ, ротационные операторы имеют импорт из другого файла с именем Thin3dtemplates. Я рекомендую вам загрузить rotatealOperators, Thin3dtemplates, преобразовать файлы сценариев Python, а также загрузить файл lookuparray.npy, который используется в качестве таблицы поиска в формате массива numpy, предварительно рассчитанного для проверки вокселя для маркировки для удаления или нет. Для запуска этих кодов вам потребуется версия python > 3, модули scipy, numpy и pyeda.

import numpy as np
import time
from scipy import ndimage
from scipy.ndimage.filters import convolve

"""
   the following subiteration functions are how each image is rotated to the next direction for removing
   boundary voxels in the order described in the reference paper
   us, ne, wd,..
"""
from rotationalOperators import firstSubiteration, secondSubiteration, thirdSubiteration, fourthSubiteration, fifthSubiteration, sixthSubiteration, seventhSubiteration, eighthSubiteration, ninthSubiteration, tenthSubiteration, eleventhSubiteration, twelvethSubiteration

"""
   reference paper
   http://web.inf.u-szeged.hu/ipcg/publications/papers/PalagyiKuba_GMIP1999.pdf
   input should be a binary image/ already segmented
"""


"""
   array that has calculated the validity of the 14 templates beforehand and stored each index which is
   decimal number of the binary string of 26 values (sqrt(3) connectivity) that are around a single voxel 
"""

lookUpTablearray = np.load('lookupTablearray.npy')


def _convolveImage(arr, flippedKernel):
    arr = np.ascontiguousarray(arr, dtype=np.uint64)
    result = convolve(arr, flippedKernel, mode='constant', cval=0)
    result[arr == 0] = 0
    return result


"""
each of the 12 iterations corresponds to each of the following
directions - us, ne, wd, es, uw, nd, sw, un, ed, nw, ue, sd
imported from template expressions
evaluated in advance using pyeda
https://pyeda.readthedocs.org/en/latest/expr.html
"""

sElement = ndimage.generate_binary_structure(3, 1)


def _getBouondariesOfimage(image):
    """
       function to find boundaries/border/edges of the array/image
    """

    erode_im = ndimage.morphology.binary_erosion(image, sElement)
    boundaryIm = image - erode_im
    return boundaryIm

"""
each of the 12 iterations corresponds to each of the following
directions - us, ne, wd, es, uw, nd, sw, un, ed, nw, ue, sd
imported from template expressions
evaluated in advance using pyeda
https://pyeda.readthedocs.org/en/latest/expr.html
"""

directionList = [firstSubiteration, secondSubiteration, thirdSubiteration, fourthSubiteration,
                 fifthSubiteration, sixthSubiteration, seventhSubiteration, eighthSubiteration,
                 ninthSubiteration, tenthSubiteration, eleventhSubiteration, twelvethSubiteration]


def _skeletonPass(image):
    """
        each pass consists of 12 serial subiterations and finding the
        boundaries of the padded image/array
    """
    boundaryIm = _getBouondariesOfimage(image)
    numPixelsremovedList = [] * 12
    boundaryIndices = list(set(map(tuple, list(np.transpose(np.nonzero(boundaryIm))))))
    for i in range(0, 12):
        convImage = _convolveImage(image, directionList[i])
        totalPixels, image = _applySubiter(image, boundaryIndices, convImage)
        print("number of pixels removed in the {} direction is {}". format(i, totalPixels))
        numPixelsremovedList.append(totalPixels)
    numPixelsremoved = sum(numPixelsremovedList)
    return numPixelsremoved, image


def _applySubiter(image, boundaryIndices, convImage):
    """
       each subiteration paralleley reduces the border voxels in 12 directions
       going through each voxel and marking if it can be deleted or not in a
       different image named temp_del and finally multiply it with the original
       image to delete the voxels so marked
    """
    temp_del = np.zeros_like(image)
    # boundaryIndicesCopy = copy.deepcopy(boundaryIndices)
    lenB = len(boundaryIndices)
    for k in range(0, lenB):
        temp_del[boundaryIndices[k]] = lookUpTablearray[convImage[boundaryIndices[k]]]
    numpixel_removed = np.einsum('ijk->', image * temp_del, dtype=int)
    image[temp_del == 1] = 0
    return numpixel_removed, image


def getSkeletonize3D(image):
    """
    function to skeletonize a 3D binary image with object in brighter contrast than background.
    In other words, 1 = object, 0 = background
    """
    assert np.max(image) in [0, 1]
    zOrig, yOrig, xOrig = np.shape(image)
    padImage = np.lib.pad(image, 1, 'constant', constant_values=0)
    start_skeleton = time.time()
    pass_no = 0
    numpixel_removed = 0
    while pass_no == 0 or numpixel_removed > 0:
        numpixel_removed, padImage = _skeletonPass(padImage)
        print("number of pixels removed in pass {} is {}".format(pass_no, numpixel_removed))
        pass_no += 1
    print("done %i number of pixels in %f seconds" % (np.sum(image), time.time() - start_skeleton))
    return padImage[1: zOrig + 1, 1: yOrig + 1, 1: xOrig + 1]

if __name__ == '__main__':
    sample = np.ones((5, 5, 5), dtype=np.uint8)
    resultSkel = getSkeletonize3D(sample)
    # gives a single voxel at the center
    print("resultSkel", resultSkel)
person Pranathi    schedule 18.12.2015