<s> r"""<STR_LIT>""" <EOL> import scipy as sp <EOL> from transforms3d import _gohlketransforms as tr <EOL> from OpenPNM . Base import logging <EOL> logger = logging . getLogger ( __name__ ) <EOL> def distance_transform ( network , geometry , offset , ** kwargs ) : <EOL> r"""<STR_LIT>""" <EOL> import math <EOL> import numpy as np <EOL> from skimage . morphology import convex_hull_image <EOL> from skimage . measure import regionprops <EOL> from scipy import ndimage <EOL> Nt = geometry . num_throats ( ) <EOL> area = sp . zeros ( Nt ) <EOL> perimeter = sp . zeros ( Nt ) <EOL> centroid = sp . zeros ( [ Nt , <NUM_LIT:3> ] ) <EOL> incentre = sp . zeros ( [ Nt , <NUM_LIT:3> ] ) <EOL> inradius = sp . zeros ( Nt ) <EOL> equiv_diameter = sp . zeros ( Nt ) <EOL> eroded_verts = sp . ndarray ( Nt , dtype = object ) <EOL> res = <NUM_LIT:200> <EOL> vertices = geometry [ '<STR_LIT>' ] <EOL> normals = geometry [ '<STR_LIT>' ] <EOL> z_axis = [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ] <EOL> for i in range ( Nt ) : <EOL> logger . info ( "<STR_LIT>" + str ( i + <NUM_LIT:1> ) + "<STR_LIT>" + str ( Nt ) ) <EOL> angle = tr . angle_between_vectors ( normals [ i ] , z_axis ) <EOL> if angle == <NUM_LIT:0.0> or angle == np . pi : <EOL> rotate_facet = False <EOL> facet = vertices [ i ] <EOL> else : <EOL> rotate_facet = True <EOL> M = tr . rotation_matrix ( tr . angle_between_vectors ( normals [ i ] , z_axis ) , <EOL> tr . vector_product ( normals [ i ] , z_axis ) ) <EOL> facet = np . dot ( vertices [ i ] , M [ : <NUM_LIT:3> , : <NUM_LIT:3> ] . T ) <EOL> x = facet [ : , <NUM_LIT:0> ] <EOL> y = facet [ : , <NUM_LIT:1> ] <EOL> z = facet [ : , <NUM_LIT:2> ] <EOL> pts = np . column_stack ( ( x , y ) ) <EOL> translation = [ pts [ : , <NUM_LIT:0> ] . min ( ) , pts [ : , <NUM_LIT:1> ] . min ( ) ] <EOL> pts -= translation <EOL> order = np . int ( math . ceil ( - np . log10 ( np . max ( pts ) ) ) ) <EOL> max_factor = np . max ( [ pts [ : , <NUM_LIT:0> ] . max ( ) , pts [ : , <NUM_LIT:1> ] . max ( ) ] ) <EOL> f = res / max_factor <EOL> r = f * offset <EOL> if r <= res / <NUM_LIT:2> : <EOL> pts *= f <EOL> minp1 = pts [ : , <NUM_LIT:0> ] . min ( ) <EOL> minp2 = pts [ : , <NUM_LIT:1> ] . min ( ) <EOL> maxp1 = pts [ : , <NUM_LIT:0> ] . max ( ) <EOL> maxp2 = pts [ : , <NUM_LIT:1> ] . max ( ) <EOL> img = np . zeros ( [ np . int ( math . ceil ( maxp1 - minp1 ) + <NUM_LIT:1> ) , <EOL> np . int ( math . ceil ( maxp2 - minp2 ) + <NUM_LIT:1> ) ] ) <EOL> int_pts = np . around ( pts , <NUM_LIT:0> ) . astype ( int ) <EOL> for pt in int_pts : <EOL> img [ pt [ <NUM_LIT:0> ] ] [ pt [ <NUM_LIT:1> ] ] = <NUM_LIT:1> <EOL> img_pad = np . zeros ( [ np . shape ( img ) [ <NUM_LIT:0> ] + <NUM_LIT:2> , np . shape ( img ) [ <NUM_LIT:1> ] + <NUM_LIT:2> ] ) <EOL> img_pad [ <NUM_LIT:1> : np . shape ( img ) [ <NUM_LIT:0> ] + <NUM_LIT:1> , <NUM_LIT:1> : np . shape ( img ) [ <NUM_LIT:1> ] + <NUM_LIT:1> ] = img <EOL> z_plane = sp . unique ( np . around ( z , order + <NUM_LIT:2> ) ) <EOL> if len ( z_plane ) > <NUM_LIT:1> : <EOL> logger . error ( '<STR_LIT>' ) <EOL> temp_arr = np . ones ( <NUM_LIT:1> ) <EOL> temp_arr . fill ( np . mean ( z_plane ) ) <EOL> z_plane = temp_arr <EOL> "<STR_LIT>" <EOL> convhullimg = convex_hull_image ( img_pad ) <EOL> eroded = ndimage . distance_transform_edt ( convhullimg ) <EOL> eroded [ eroded <= r ] = <NUM_LIT:0> <EOL> eroded [ eroded > r ] = <NUM_LIT:1> <EOL> if np . sum ( eroded ) >= <NUM_LIT:3> : <EOL> regions = regionprops ( eroded [ <NUM_LIT:1> : np . shape ( img ) [ <NUM_LIT:0> ] + <NUM_LIT:1> , <EOL> <NUM_LIT:1> : np . shape ( img ) [ <NUM_LIT:1> ] + <NUM_LIT:1> ] . astype ( int ) ) <EOL> if len ( regions ) == <NUM_LIT:1> : <EOL> for props in regions : <EOL> x0 , y0 = props . centroid <EOL> equiv_diameter [ i ] = props . equivalent_diameter <EOL> area [ i ] = props . area <EOL> perimeter [ i ] = props . perimeter <EOL> coords = props . coords <EOL> centroid2d = [ x0 , y0 ] / f <EOL> centroid2d += ( translation ) <EOL> centroid3d = np . concatenate ( ( centroid2d , z_plane ) ) <EOL> dt = ndimage . distance_transform_edt ( eroded ) <EOL> inx0 , iny0 = np . asarray ( np . unravel_index ( dt . argmax ( ) , dt . shape ) ) . astype ( float ) <EOL> incentre2d = [ inx0 , iny0 ] <EOL> incentre2d /= f <EOL> incentre2d += ( translation ) <EOL> incentre3d = np . concatenate ( ( incentre2d , z_plane ) ) <EOL> offset_verts = [ ] <EOL> for pt in int_pts : <EOL> vert = np . argmin ( np . sum ( np . square ( coords - pt ) , axis = <NUM_LIT:1> ) ) <EOL> if vert not in offset_verts : <EOL> offset_verts . append ( vert ) <EOL> if len ( offset_verts ) >= <NUM_LIT:3> : <EOL> offset_coords = coords [ offset_verts ] . astype ( float ) <EOL> offset_coords /= f <EOL> offset_coords_3d = np . vstack ( ( offset_coords [ : , <NUM_LIT:0> ] + translation [ <NUM_LIT:0> ] , <EOL> offset_coords [ : , <NUM_LIT:1> ] + translation [ <NUM_LIT:1> ] , <EOL> np . ones ( len ( offset_verts ) ) * z_plane ) ) . T <EOL> if rotate_facet : <EOL> MI = tr . inverse_matrix ( M ) <EOL> incentre [ i ] = np . dot ( incentre3d , MI [ : <NUM_LIT:3> , : <NUM_LIT:3> ] . T ) <EOL> centroid [ i ] = np . dot ( centroid3d , MI [ : <NUM_LIT:3> , : <NUM_LIT:3> ] . T ) <EOL> eroded_verts [ i ] = np . dot ( offset_coords_3d , MI [ : <NUM_LIT:3> , : <NUM_LIT:3> ] . T ) <EOL> else : <EOL> incentre [ i ] = incentre3d <EOL> centroid [ i ] = centroid3d <EOL> eroded_verts [ i ] = offset_coords_3d <EOL> inradius [ i ] = dt . max ( ) <EOL> area [ i ] /= f * f <EOL> perimeter [ i ] /= f <EOL> equiv_diameter [ i ] /= f <EOL> inradius [ i ] /= f <EOL> else : <EOL> area [ i ] = <NUM_LIT:0> <EOL> perimeter [ i ] = <NUM_LIT:0> <EOL> equiv_diameter [ i ] = <NUM_LIT:0> <EOL> if kwargs [ '<STR_LIT>' ] is True : <EOL> geometry [ '<STR_LIT>' ] = area <EOL> geometry [ '<STR_LIT>' ] = perimeter <EOL> geometry [ '<STR_LIT>' ] = centroid <EOL> geometry [ '<STR_LIT>' ] = equiv_diameter <EOL> geometry [ '<STR_LIT>' ] = inradius * <NUM_LIT:2> <EOL> geometry [ '<STR_LIT>' ] = incentre <EOL> return eroded_verts </s>
