Хороший алгоритм вычисления объема или площади поверхности в python

Я пытаюсь вычислить объем (или площадь поверхности) массива 3D numpy. Воксели во многих случаях анизотропны, и у меня есть коэффициент преобразования пикселей в пиксель в каждом направлении.Хороший алгоритм вычисления объема или площади поверхности в python

Кто-нибудь знает хорошее место, чтобы найти набор инструментов или пакет, чтобы сделать выше?

Прямо сейчас, у меня есть код внутри дома, но я хочу перейти на более прочную прочность с точки зрения точности.

Редактировать 1: Вот некоторые (плохие) образцы data. Это намного меньше, чем типичная сфера. Я добавлю лучшие данные, когда смогу его создать! Он находится в (я). TumorBrain.tumor.

источник

2012-01-18 tylerthemiler

Вы должны определить, что вы имеете в виду немного более четко. Если ваши данные на самом деле представляют собой трехмерный массив, то объем, который занимает вся сетка, это 'nx * dx * ny * dy * nz * dz', но я уверен, что вы не имели в виду это ... Вы хотите, чтобы объем внутри изоповерхности? –

Я ДУМАЮ, что ты прав. Это двоичный массив X x Y x Z, и мне нужен объем всего, что содержится в периметре его двоичной части. Это обычно (но не всегда) сфера. – tylerthemiler

это звучит забавно, можете ли вы опубликовать ссылку на некоторые примеры данных? просто сохраните массив numpy с помощью 'pickle.dump' – wim

Один из вариантов - использовать VTK. (Я собираюсь использовать привязки python tvtk для этого здесь ...)

По крайней мере, в некоторых случаях получение области внутри изоповерхности будет более точным.

Кроме того, насколько площадь поверхности идет, tvtk.MassProperties также вычисляет площадь поверхности. Это mass.surface_area (с объектом mass в коде ниже).

import numpy as np 
from tvtk.api import tvtk 

def main(): 
    # Generate some data with anisotropic cells... 
    # x,y,and z will range from -2 to 2, but with a 
    # different (20, 15, and 5 for x, y, and z) number of steps 
    x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j] 
    r = np.sqrt(x**2 + y**2 + z**2) 

    dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))] 

    # Your actual data is a binary (logical) array 
    max_radius = 1.5 
    data = (r <= max_radius).astype(np.int8) 

    ideal_volume = 4.0/3 * max_radius**3 * np.pi 
    coarse_volume = data.sum() * dx * dy * dz 
    est_volume = vtk_volume(data, (dx, dy, dz), (x.min(), y.min(), z.min())) 

    coarse_error = 100 * (coarse_volume - ideal_volume)/ideal_volume 
    vtk_error = 100 * (est_volume - ideal_volume)/ideal_volume 

    print 'Ideal volume', ideal_volume 
    print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%' 
    print 'VTK approximation', est_volume, 'Error', vtk_error, '%' 

def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)): 
    data[data == 0] = -1 
    grid = tvtk.ImageData(spacing=spacing, origin=origin) 
    grid.point_data.scalars = data.T.ravel() # It wants fortran order??? 
    grid.point_data.scalars.name = 'scalars' 
    grid.dimensions = data.shape 

    iso = tvtk.ImageMarchingCubes(input=grid) 
    mass = tvtk.MassProperties(input=iso.output) 
    return mass.volume 

main()

Это дает:

Ideal volume 14.1371669412 
Coarse approximation 14.7969924812 Error 4.66731094565 % 
VTK approximation 14.1954890878 Error 0.412544796894 %

источник

2012-01-18 06:42:43

Это выглядит потрясающе! Я попробую, когда мой босс вернется (у меня даже нет админа, поэтому я не могу установить vtk :(). – tylerthemiler

Voxels будут довольно простыми, правильными многогранниками, нет? Вычислите объем каждого из них и суммируйте их.

источник

2012-01-18 00:52:23

Основная проблема заключается в том, что направление z слишком грубое. Мы делаем своего рода метод shells для усреднения между срезами, но это недостаточно. – tylerthemiler

Хороший алгоритм вычисления объема или площади поверхности в python

ответ

Смежные вопросы