-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmaxProjVol.py
101 lines (100 loc) · 6.06 KB
/
maxProjVol.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
import numpy as np
import sys, random, nrrd
if (len(sys.argv) < 2):
print('Error: missing arguments!')
print('e.g. python maxProjVol.py file.nrrd [outfile.nrrd]')
else:
filein = str(sys.argv[1]).replace('.obj', '.nrrd')
print('Loading %s...' % (filein))
data1, header1 = nrrd.read(filein)
out = str(sys.argv[1])
shape = np.shape(data1)
if len(sys.argv) > 2:
out = str(sys.argv[2])
data2 = np.zeros(shape, dtype=np.uint8)
max = np.max(data1)
th = 255
vertices = '####\n#\n# OBJ File Generated by VirtualFlyBrain.org\n#\n####\n'
print('Creating 3D max projection...')
c = 0
p = 1.0
binCount = np.nonzero(np.bincount(data1.flatten()))[0].size
print(binCount)
while ((c < 15000 and th > 10) or (p < 10 and c < 250000)) and (binCount > 3 or c < 10):
th = max / p
for x in [j for (i, j) in zip(np.max(np.max(data1, axis=2), axis=1), range(0, shape[0])) if i >= th]:
# print(str('{0:.1f}'.format((np.float(x) / shape[0]) * 100.0)) + '%')
for y in [j for (i, j) in zip(np.max(data1[x, :, :], axis=1), range(0, shape[1])) if i >= th]:
l = [data1[x, y, :].argmax()]
if l[0] < shape[2]-2:
l += [shape[2]-1-(data1[x, y, :l[0]+1:-1].argmax())]
for z in l:
try:
if data1[x, y, z] > th:
c += 1
data2[x, y, z] = data1[x, y, z]
radius = data1[x, y, z] / 256.0
vertices += 'v ' + str((x * np.float(header1['space directions'][0][0]))+(random.uniform(-1, 1)*np.float(header1['space directions'][0][0]))) + ' '
vertices += str((y * np.float(header1['space directions'][1][1]))+(random.uniform(-1, 1)*np.float(header1['space directions'][1][1]))) + ' '
vertices += str((z * np.float(header1['space directions'][2][2]))+(random.uniform(-1, 1)*np.float(header1['space directions'][2][2]))) + ' ' + str(radius) + '\n'
if binCount < 4:
data1[x:x+1, y:y+1, z:z+1] = np.uint(0)
else:
data2[x, y, z] = np.uint8(0)
except (RuntimeError, ValueError):
data2[x, y, z] = np.uint8(0)
print([x, y, z])
except IndexError:
print([x, y, z])
print(str(c) + ' points found with threshold of ' + str(th))
p += 1.0
# adding full surface for binary images (3 due to potential boundary markers)
if binCount < 4:
for z in [j for (i, j) in zip(np.max(np.max(data1, axis=1), axis=0), range(0, shape[2])) if i >= th]:
for y in [j for (i, j) in zip(np.max(data1[:, :, z], axis=0), range(0, shape[1])) if i >= th]:
l = [data1[:, y, z].argmax()]
if l[0] < shape[0] - 2:
l += [shape[0] - 1 - (data1[:l[0] + 1:-1, y, z].argmax())]
for x in l:
try:
if data1[x, y, z] > th:
c += 1
data2[x, y, z] = data1[x, y, z]
radius = data1[x, y, z] / 256.0
vertices += 'v ' + str((x * np.float(header1['space directions'][0][0]))+(random.uniform(-1, 1)*np.float(header1['space directions'][0][0]))) + ' '
vertices += str((y * np.float(header1['space directions'][1][1]))+(random.uniform(-1, 1)*np.float(header1['space directions'][1][1]))) + ' '
vertices += str((z * np.float(header1['space directions'][2][2]))+(random.uniform(-1, 1)*np.float(header1['space directions'][2][2]))) + ' ' + str(radius) + '\n'
data1[x:x+1, y:y+1, z:z+1] = np.uint(0)
except (RuntimeError, ValueError):
data2[x, y, z] = np.uint8(0)
print([x, y, z])
except IndexError:
print([x, y, z])
print(str(c) + ' points found with threshold of ' + str(th))
for z in [j for (i, j) in zip(np.max(np.max(data1, axis=1), axis=0), range(0, shape[2])) if i >= th]:
for x in [j for (i, j) in zip(np.max(data1[:, :, z], axis=1), range(0, shape[0])) if i >= th]:
l = [data1[x, :, z].argmax()]
if l[0] < shape[1] - 2:
l += [shape[1] - 1 - (data1[x, :l[0] + 1:-1, z].argmax())]
for y in l:
try:
if data1[x, y, z] > th:
c += 1
data2[x, y, z] = data1[x, y, z]
radius = data1[x, y, z] / 256.0
vertices += 'v ' + str((x * np.float(header1['space directions'][0][0]))+(random.uniform(-1, 1)*np.float(header1['space directions'][0][0]))) + ' '
vertices += str((y * np.float(header1['space directions'][1][1]))+(random.uniform(-1, 1)*np.float(header1['space directions'][1][1]))) + ' '
vertices += str((z * np.float(header1['space directions'][2][2]))+(random.uniform(-1, 1)*np.float(header1['space directions'][2][2]))) + ' ' + str(radius) + '\n'
data1[x:x+1, y:y+1, z:z+1] = np.uint(0)
except (RuntimeError, ValueError):
data2[x, y, z] = np.uint8(0)
print([x, y, z])
except IndexError:
print([x, y, z])
print(str(c) + ' points found with threshold of ' + str(th))
if len(sys.argv) > 2:
print("Saving result to nrrd " + out.replace('.nrrd', '').replace('.obj', ''))
nrrd.write(out.replace('.obj', '.nrrd').replace('.nrrd', '_max.nrrd'), np.uint8(data2), header=header1)
print("Saving result to obj " + out.replace('.nrrd', '').replace('.obj', ''))
with open(out.replace('.nrrd', '.obj'), 'w+') as objfile:
objfile.write(vertices)