fixed serious bug regarding wrong reshaping order (was 'C' now 'F') of 3dim to 1dim and back.

This commit is contained in:
Philip Eisenlohr 2016-04-28 13:24:08 -04:00
parent b9021f3b3f
commit d3a7ceff15
1 changed files with 21 additions and 13 deletions

View File

@ -159,7 +159,9 @@ for name in filenames:
if table.label_dimension(options.id) != 1: errors.append('grain identifier {} not found.'.format(options.id)) if table.label_dimension(options.id) != 1: errors.append('grain identifier {} not found.'.format(options.id))
else: idCol = table.label_index(options.id) else: idCol = table.label_index(options.id)
if remarks != []: damask.util.croak(remarks) if remarks != []:
damask.util.croak(remarks)
remarks = []
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
@ -184,6 +186,8 @@ for name in filenames:
N = grid.prod() N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}.'.format(N,'x'.join(map(str,grid)))) if N != len(table.data): errors.append('data count {} does not match grid {}.'.format(N,'x'.join(map(str,grid))))
else: remarks.append('grid: {}x{}x{}'.format(*grid))
if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
@ -194,33 +198,37 @@ for name in filenames:
stack = [table.data] stack = [table.data]
neighborhood = neighborhoods[options.neighborhood] neighborhood = neighborhoods[options.neighborhood]
convoluted = np.empty([len(neighborhood)]+list(grid+2),'i') diffToNeighbor = np.empty(list(grid+2)+[len(neighborhood)],'i')
microstructure = periodic_3Dpad(np.array(table.data[:,idCol].reshape(grid),'i')) microstructure = periodic_3Dpad(table.data[:,idCol].astype('i').reshape(grid,order='F'))
for i,p in enumerate(neighborhood): for i,p in enumerate(neighborhood):
stencil = np.zeros((3,3,3),'i') stencil = np.zeros((3,3,3),'i')
stencil[1,1,1] = -1 stencil[1,1,1] = -1
stencil[p[0]+1, stencil[p[0]+1,
p[1]+1, p[1]+1,
p[2]+1] = 1 p[2]+1] = 1
convoluted[i,:,:,:] = ndimage.convolve(microstructure,stencil) diffToNeighbor[:,:,:,i] = ndimage.convolve(microstructure,stencil) # compare ID at each point...
# ...to every one in the specified neighborhood
# for same IDs at both locations ==> 0
distance = np.ones((len(feature_list),grid[0],grid[1],grid[2]),'d') diffToNeighbor = np.sort(diffToNeighbor) # sort diff such that number of changes in diff (steps)...
# ...reflects number of unique neighbors
convoluted = np.sort(convoluted,axis = 0) uniques = np.where(diffToNeighbor[1:-1,1:-1,1:-1,0] != 0, 1,0) # initialize unique value counter (exclude myself [= 0])
uniques = np.where(convoluted[0,1:-1,1:-1,1:-1] != 0, 1,0) # initialize unique value counter (exclude myself [= 0])
for i in xrange(1,len(neighborhood)): # check remaining points in neighborhood for i in xrange(1,len(neighborhood)): # check remaining points in neighborhood
uniques += np.where(np.logical_and( uniques += np.where(np.logical_and(
convoluted[i,1:-1,1:-1,1:-1] != convoluted[i-1,1:-1,1:-1,1:-1], # flip of ID difference detected? diffToNeighbor[1:-1,1:-1,1:-1,i] != 0, # not myself?
convoluted[i,1:-1,1:-1,1:-1] != 0), # not myself? diffToNeighbor[1:-1,1:-1,1:-1,i] != diffToNeighbor[1:-1,1:-1,1:-1,i-1],
1,0) # count flip ), # flip of ID difference detected?
1,0) # count that flip
distance = np.ones((len(feature_list),grid[0],grid[1],grid[2]),'d')
for i,feature_id in enumerate(feature_list): for i,feature_id in enumerate(feature_list):
distance[i,:,:,:] = np.where(uniques >= features[feature_id]['aliens'],0.0,1.0) # seed with 0.0 when enough unique neighbor IDs are present distance[i,:,:,:] = np.where(uniques >= features[feature_id]['aliens'],0.0,1.0) # seed with 0.0 when enough unique neighbor IDs are present
distance[i,:,:,:] = ndimage.morphology.distance_transform_edt(distance[i,:,:,:])*[options.scale]*3 distance[i,:,:,:] = ndimage.morphology.distance_transform_edt(distance[i,:,:,:])*[options.scale]*3
distance.shape = ([len(feature_list),grid.prod(),1]) distance = distance.reshape([len(feature_list),grid.prod(),1],order='F')
for i in xrange(len(feature_list)): for i in xrange(len(feature_list)):
stack.append(distance[i,:]) stack.append(distance[i,:])