forked from andreubs/mammoreplicator
-
Notifications
You must be signed in to change notification settings - Fork 0
/
mammoreplicator.py
256 lines (198 loc) · 12 KB
/
mammoreplicator.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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
#!/usr/bin/python
#
# ** Reproducing 2D mammograms with 3D printed surrogates **
# [Andreu Badal, 2015-07-10, 2015-02-03]
#
# Process a DICOM mammography to generate a printable triangle mesh
#
#
import numpy
from numpy import *
# Instal DICOM reader for Python: http://www.pydicom.org/ (just run ez_setup.py as root)
import dicom
def str4(n):
return "{0:.4f}".format(n) # Simple function to convert a float into a string with 4 decimal places
# Input data:
print '\n** "mammoreplicator": Replicating 2D mammograms with 3D printed surrogates **'
print ' [Andreu Badal, 2015-07-10]\n'
print ' Reading input data...'
image_file_name = "TG18-MM-2k-01.dcm"
compression_thickness = 5.0 # cm
average_breast_intensity = 700
output_binning = 8 # Rebin the pixels to reduce the print resolution (1=no-binning, 2=average_4pixels...)
subtract_layer = 2.0
output_dcm = 0 # 1==true, 0==false->do not output dicom file with thickness
output_mm = 1 # 1==true, 0==false->output triangles in cm
rows = [0,2048]
columns = [0,2048]
num_rows = int((rows[1]-rows[0])/output_binning)*output_binning # Processed sizwe will be a multiple of binning value
num_columns= int((columns[1]-columns[0])/output_binning)*output_binning
rows[1] = rows[0]+num_rows
columns[1] = columns[0]+num_columns
center_pixel = array([num_columns/2,num_rows]) # Location x-ray field of view center in pixel units [x,y]=[Row,Column]
center_coord = array([0.0,0.0, 0.0]) # Location x-ray field of view center in cm (printed phantom origin)
source_coord = array([0.0,0.0,80.0]) # Location x-ray source focal spot. Using NUMPY arrays for calculations.
pixel_size = 0.0070 # cm # Desired pixel size in the printed phantom
mfp_breast = 1.284 # cm # 16.5keV=1.284cm, 20keV=1.933cm; // Molybdenum target, 28 kVp -> average energy 16.5 keV
#mfp_plastic = 0.863 # cm # 16.5keV=0.863cm, 20keV=1.315cm
mfp_plastic = 1.0; # cm --> wild estimation assuming some combination of PMMA and TiO2 in the stereolithogaphy substrate !!STL!!
air_threshold = 0 #!!DeBuG!! Currently not removing air pixels: drawing 0 height columns instead # Any pixel below this threshold will be considered air and asigned 0 thickness
print " compression_thickness =", compression_thickness
print " average_breast_intensity =", average_breast_intensity
print " columns=["+str(columns[0])+","+str(columns[1])+"], num_columns = ", num_columns
print " rows=["+str(rows[0])+","+str(rows[1])+"], num_rows = ", num_rows
#print " max_pixel_value =", max_pixel_value
print " air_threshold =", air_threshold
print " mfp_breast =", mfp_breast
print " mfp_plastic =", mfp_plastic
print " center_pixel =", center_pixel
print " center_coord =", center_coord
print " source_coord =", source_coord
print " pixel_size =", pixel_size
print " subtract_layer =", subtract_layer
print " output_binning =", output_binning
print " binned columns = "+str(num_columns/output_binning)+", binned rows = "+str(num_rows/output_binning)
print " "
# Read image data:
print ' Reading input file: ',image_file_name,'...'
mammo = dicom.read_file(image_file_name)
# Process pixel values:
print '\n Processing pixel values; rows = ',num_rows,', columns = ',num_columns,'...'
disp( ' Current row: ', linefeed=False) # (disable new line at end of message)
pix = numpy.zeros((rows[1], columns[1])) # We will work on this copy of the data as floats, init to 0
num_air_pix = 0
conversion_factor = (compression_thickness/log(float(average_breast_intensity))) * (mfp_plastic/mfp_breast) # !!DeBuG!! Linearity assumption
for j in range(rows[0], rows[1], output_binning):
disp(j, linefeed=False); disp(', ', linefeed=False)
for i in range(columns[0], columns[1], output_binning):
# If requested, re-bin the output to reduce the print resolution (1=no-binning, 2=average_4pixels...):
if output_binning>1:
avg = 0.0
for k1 in range(0, output_binning):
for k2 in range(0, output_binning):
#print 'j='+str(j)+' i='+str(i)+': j+k1='+str(j+k1)+' i+k2='+str(i+k2) # !!Verbose!!
avg = avg + float(mammo.pixel_array[j+k1][i+k2])
mammo.pixel_array[j][i] = int(round(avg/(output_binning**2)))
# Skip air pixels (leave them as 0 thickness):
if mammo.pixel_array[j][i]<air_threshold:
num_air_pix = num_air_pix + 1
pix[j][i] = 0
else:
if (mammo.pixel_array[j][i]>1):
pix[j][i] = log(float(mammo.pixel_array[j][i])) * conversion_factor #!!DeBuG!! Calculate log of image?
else:
pix[j][i] = 0
if (output_dcm==1):
bin_value = int(round(pix[j][i]*100)) # Convert cm floats into tenths of mm integers to store in dicom image
for k1 in range(0, output_binning):
for k2 in range(0, output_binning):
mammo.pixel_array[j+k1][i+k2] = bin_value
if (output_dcm==1):
# Save final data as dicom: help(dicom.write_file)
print '\n\n Saving final data as DICOM image in integer units of tenths of mm of plastic...'
mammo.PixelData = mammo.pixel_array.tostring() # Save the data array into the actual pixel data string
#dicom.write_file(image_file_name+"_plastic_bin"+str(output_binning)+".dcm", mammo)
dicom.write_file(image_file_name+"_STL_bin"+str(output_binning)+".dcm", mammo) # !!STL!!
# Generate triangle mesh:
print '\n Generating the output triangle mesh focused to the x-ray source focal spot (num_air_pix='+str(num_air_pix)+')...'
num_binned_pixels = (num_columns/output_binning) * (num_rows/output_binning) - num_air_pix
num_vertices = num_binned_pixels * 8 # 8 vertices per printed column
num_triangles = num_binned_pixels * 8 + 2*(num_rows/output_binning) + 2*(num_columns/output_binning) #!!ReducedTriangles!!
print ' Writing ',num_triangles,' triangles in output file: '+image_file_name+"_bin"+str(output_binning)+'.ply'
image_file_name2 = image_file_name+"_STL_bin"+str(output_binning)+"_focusedZ"+str(source_coord[2]) # !!STL!!
if (subtract_layer>0.00001):
image_file_name2 = image_file_name2+"_-"+str(subtract_layer)+"cm"
if (output_mm==1):
image_file_name2 = image_file_name2+"_mm"
image_file_name2 = image_file_name2+".ply"
ply = open(image_file_name2, 'w') # Create PLY file and write header:
ply.write('ply\n')
ply.write('format ascii 1.0\n')
ply.write('comment ** Geometry created by mammoreplicator.py using the input file: '+image_file_name+'\n')
ply.write('comment [Andreu Badal, 2015-02-03]\n')
if (subtract_layer>0.00001):
ply.write('comment NOTE: subtracting '+str(subtract_layer)+' cm from each column for easier printing.\n')
if (output_mm==1):
ply.write('comment NOTE: vertex coordinates in mm.\n')
ply.write('element vertex '+str(num_vertices)+'\n')
ply.write('property float x\n')
ply.write('property float y\n')
ply.write('property float z\n')
ply.write('element face '+str(num_triangles)+'\n')
ply.write('property list uchar int vertex_index\n')
ply.write('end_header\n')
disp( ' Current row: ', linefeed=False) # (disable new line at end of message)
for j in range(rows[0], rows[1], output_binning):
disp(j, linefeed=False); disp(', ', linefeed=False)
for i in range(columns[0], columns[1], output_binning):
if pix[j][i]>-999990.01: # Skip air pixels
# Define vectors:
r0 = numpy.zeros(3)
#r1 = numpy.zeros(3)
d = numpy.zeros(3)
vertices = numpy.zeros((8,3))
offset = 0.5*pixel_size*output_binning
offset0 = array([+offset,+offset,0.0]) # 3 ----- 0 (vertex order clockwise)
offset1 = array([+offset,-offset,0.0]) # - c -
offset2 = array([-offset,-offset,0.0]) # 2 ----- 1
offset3 = array([-offset,+offset,0.0])
# Set pixel center coordinates in real space:
r0[0] = center_coord[0] + float(i-center_pixel[0]) * pixel_size
r0[1] = center_coord[1] + float(j-center_pixel[1]) * pixel_size
r0[2] = center_coord[2]
vertices[0] = r0 + offset0 # 4 vertices at the base
vertices[1] = r0 + offset1
vertices[2] = r0 + offset2
vertices[3] = r0 + offset3
# - Shorten object for printing (and make sure all pixels are 0 or more):
pix[j][i] = max(pix[j][i] - subtract_layer, 0.0)
d = source_coord-vertices[0] # 4 vertices at the top focused to the source: calculate mormalized direction vector for each vertex
vertices[4] = vertices[0] + pix[j][i]*d/numpy.linalg.norm(d)
d = source_coord-vertices[1]
vertices[5] = vertices[1] + pix[j][i]*d/numpy.linalg.norm(d)
d = source_coord-vertices[2]
vertices[6] = vertices[2] + pix[j][i]*d/numpy.linalg.norm(d)
d = source_coord-vertices[3]
vertices[7] = vertices[3] + pix[j][i]*d/numpy.linalg.norm(d)
# Write vertices coordinates:
for n in range(0, 8):
if (output_mm==1):
ply.write(str4(10.0*vertices[n][0])+' '+str4(10.0*vertices[n][1])+' '+str4(10.0*vertices[n][2])+'\n') # Outputing mesh in mm
else:
ply.write(str4(vertices[n][0])+' '+str4(vertices[n][1])+' '+str4(vertices[n][2])+'\n') # Outputing mesh in cm
# Write triangles in the 6 sides of each printed column giving groups of 3 vertices:
for n in range(0, num_binned_pixels):
# - Create a triangle mesh combining the walls of the columns of consecutive pixels:
ply.write('3 '+str(0+n*8)+' '+str(1+n*8)+' '+str(2+n*8)+'\n')
ply.write('3 '+str(2+n*8)+' '+str(3+n*8)+' '+str(0+n*8)+'\n') # 3 ----- 0 (vertex order clockwise)
ply.write('3 '+str(4+n*8)+' '+str(5+n*8)+' '+str(6+n*8)+'\n') # - c -
ply.write('3 '+str(6+n*8)+' '+str(7+n*8)+' '+str(4+n*8)+'\n') # 2 ----- 1
if (n<(num_columns/output_binning)): # 7 ----- 4
ply.write('3 '+str(1+n*8)+' '+str(5+n*8)+' '+str(6+n*8)+'\n') # -top-
ply.write('3 '+str(6+n*8)+' '+str(2+n*8)+' '+str(1+n*8)+'\n') # 6 ----- 5
elif (n>(num_binned_pixels-(num_columns/output_binning)-1)):
ply.write('3 '+str(1+n*8+3-(num_columns/output_binning)*8)+' '+str(5+n*8)+' '+str(6+n*8)+'\n')
ply.write('3 '+str(6+n*8)+' '+str(2+n*8+5-(num_columns/output_binning)*8)+' '+str(1+n*8+3-(num_columns/output_binning)*8)+'\n')
ply.write('3 '+str(3+n*8)+' '+str(7+n*8)+' '+str(4+n*8)+'\n')
ply.write('3 '+str(4+n*8)+' '+str(0+n*8)+' '+str(3+n*8)+'\n')
else:
ply.write('3 '+str(1+n*8+3-(num_columns/output_binning)*8)+' '+str(5+n*8)+' '+str(6+n*8)+'\n')
ply.write('3 '+str(6+n*8)+' '+str(2+n*8+5-(num_columns/output_binning)*8)+' '+str(1+n*8+3-(num_columns/output_binning)*8)+'\n')
#ply.write('3 '+str(3+n*8)+' '+str(7+n*8)+' '+str(4+n*8)+'\n') # No need to build this wall, already written in past column
#ply.write('3 '+str(4+n*8)+' '+str(0+n*8)+' '+str(3+n*8)+'\n')
if ((n+1)%(num_columns/output_binning))==0: #!!ReducedTriangles!!
ply.write('3 '+str(0+n*8)+' '+str(4+n*8)+' '+str(5+n*8)+'\n')
ply.write('3 '+str(5+n*8)+' '+str(1+n*8)+' '+str(0+n*8)+'\n')
elif (n%(num_columns/output_binning))==0: #!!ReducedTriangles!!
ply.write('3 '+str(0+n*8+7+8)+' '+str(4+n*8)+' '+str(5+n*8)+'\n')
ply.write('3 '+str(5+n*8)+' '+str(1+n*8+5+8)+' '+str(0+n*8+7+8)+'\n')
ply.write('3 '+str(2+n*8)+' '+str(6+n*8)+' '+str(7+n*8)+'\n')
ply.write('3 '+str(7+n*8)+' '+str(3+n*8)+' '+str(2+n*8)+'\n')
else:
ply.write('3 '+str(0+n*8+7+8)+' '+str(4+n*8)+' '+str(5+n*8)+'\n') # Connect top right with top left next column (7) #!!ReducedTriangles!!
ply.write('3 '+str(5+n*8)+' '+str(1+n*8+5+8)+' '+str(0+n*8+7+8)+'\n')
#ply.write('3 '+str(2+n*8)+' '+str(6+n*8)+' '+str(7+n*8)+'\n') # No need to build this wall, already written in past column
#ply.write('3 '+str(7+n*8)+' '+str(3+n*8)+' '+str(2+n*8)+'\n')
ply.write(str(n)+'\n')
ply.close()
print '\n\n Done!\n'