11#!/usr/bin/python
2- """
3- Copyright (c) 2011, Nakeisha Schimke. All rights reserved.
4- Redistribution and use in source and binary forms, with or without
5- modification, are permitted provided that the following conditions are met:
6- * Redistributions of source code must retain the above copyright notice, this
7- list of conditions and the following disclaimer.
8- * Redistributions in binary form must reproduce the above copyright notice,
9- this list of conditions and the following disclaimer in the documentation
10- and/or other materials provided with the distribution.
11- * Neither the name of the The University of Tulsa nor the names of its
12- contributors may be used to endorse or promote products derived from this
13- software without specific prior written permission.
14- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
15- ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
16- WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
17- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
18- ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
19- (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
20- LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
21- ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
22- (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
23- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24- """
252import numpy
263import nibabel as nb
274import sys
285import logging
296
7+
308def edge_mask (mask ):
319 """Create an edge of brain mask from a binary brain mask.
32-
10+
3311 Return a two-dimensional edge of brain mask.
3412 """
3513 brain = numpy .zeros (mask .shape [1 :])
3614 # iterate over axial
37- for i in range (0 ,mask .shape [1 ]- 1 ):
15+ for i in range (0 , mask .shape [1 ] - 1 ):
3816 # iterate over coronal
39- for k in range (mask .shape [2 ]- 1 , 0 , - 1 ):
40- brain [i ,k ] = mask [:,i , k ].any ()
41-
42- edgemask = numpy .zeros (brain .shape ,dtype = 'uint8' )
43- for u in range (1 ,brain .shape [0 ]- 2 ):
44- for v in range (1 ,brain .shape [1 ]- 2 ):
45- if brain [u ,v ] + brain [u - 1 , v ] == 1 :
46- edgemask [u ,v ] = 1
47- elif brain [u ,v ] + brain [u ,v - 1 ] == 1 :
48- edgemask [u ,v ] = 1
49- elif brain [u ,v ] + brain [u + 1 , v ] == 1 :
50- edgemask [u ,v ] = 1
51- elif brain [u ,v ] + brain [u ,v + 1 ] == 1 :
52- edgemask [u ,v ] = 1
17+ for k in range (mask .shape [2 ] - 1 , 0 , - 1 ):
18+ brain [i , k ] = mask [:, i , k ].any ()
19+
20+ edgemask = numpy .zeros (brain .shape , dtype = 'uint8' )
21+ for u in range (1 , brain .shape [0 ] - 2 ):
22+ for v in range (1 , brain .shape [1 ] - 2 ):
23+ if brain [u , v ] + brain [u - 1 , v ] == 1 :
24+ edgemask [u , v ] = 1
25+ elif brain [u , v ] + brain [u , v - 1 ] == 1 :
26+ edgemask [u , v ] = 1
27+ elif brain [u , v ] + brain [u + 1 , v ] == 1 :
28+ edgemask [u , v ] = 1
29+ elif brain [u , v ] + brain [u , v + 1 ] == 1 :
30+ edgemask [u , v ] = 1
5331 return edgemask
5432
33+
5534def convex_hull (brain ):
56- """Use Andrew's monotone chain algorithm to find the lower half of the convex hull.
57-
35+ """Use Andrew's monotone chain algorithm to find the lower half of the
36+ convex hull.
37+
5838 Return a two-dimensional convex hull.
5939 """
6040 # convert brain to a list of points
6141 nz = numpy .nonzero (brain )
6242 # transpose so we get an n x 2 matrix where n_i = (x,y)
63- pts = numpy .array ([nz [0 ],nz [1 ]]).transpose ()
64- pts_t = pts .transpose ()
65-
66- def cross (o ,a ,b ):
67- return (a [0 ] - o [0 ])* (b [1 ]- o [1 ]) - (a [1 ] - o [1 ]) * (b [0 ] - o [0 ])
68-
43+ pts = numpy .array ([nz [0 ], nz [1 ]]).transpose ()
44+
45+ def cross (o , a , b ):
46+ return (a [0 ] - o [0 ]) * (b [1 ] - o [1 ]) - (a [1 ] - o [1 ]) * (b [0 ] - o [0 ])
47+
6948 lower = []
70- for i in range (0 ,pts .shape [0 ]):
71- p = (pts [i ,0 ],pts [i ,1 ])
72- while len (lower ) >= 2 and cross (lower [- 2 ],lower [- 1 ],p ) <= 0 :
49+ for i in range (0 , pts .shape [0 ]):
50+ p = (pts [i , 0 ], pts [i , 1 ])
51+ while len (lower ) >= 2 and cross (lower [- 2 ], lower [- 1 ], p ) <= 0 :
7352 lower .pop ()
7453 lower .append (p )
7554
7655 return numpy .array (lower ).transpose ()
77-
56+
7857
7958def deface (anat_filename , mask_filename , defaced_filename , buff = 10 ):
8059 """Deface neuroimage using a binary brain mask.
81-
60+
8261 Keyword arguments:
8362 anat_filename -- the filename of the neuroimage to deface
8463 mask_filename -- the filename of the binary brain mask
8564 defaced_filename -- the filename of the defaced output image
86- buff -- the buffer size between the shearing line and the brain (default value is 10.0)
65+ buff -- the buffer size between the shearing line and the brain
66+ (default value is 10.0)
8767 """
8868 nii_anat = nb .load (anat_filename )
8969 nii_mask = nb .load (mask_filename )
90-
70+
9171 if numpy .equal (nii_anat .shape , nii_mask .shape ).all ():
9272 pass
9373 else :
94- logger .WARNING ("Anatomical and mask images do not have the same dimensions." )
74+ logger .warning (
75+ "Anatomical and mask images do not have the same dimensions." )
9576 sys .exit (- 1 )
9677
9778 anat_ax = nb .orientations .aff2axcodes (nii_anat .get_affine ())
9879 mask_ax = nb .orientations .aff2axcodes (nii_mask .get_affine ())
9980
100-
10181 logger .debug ("Anat image axes: {0}" .format (anat_ax ))
10282 logger .debug ("Mask image axes: {0}" .format (mask_ax ))
10383 logger .debug ("Mask shape!: {0}" .format (nii_mask .shape ))
@@ -109,7 +89,8 @@ def deface(anat_filename, mask_filename, defaced_filename, buff=10):
10989
11090 if anat_ax [0 ] != mask_ax [0 ]:
11191 # align mask to anat space
112- logger .debug ("Aligning mask to anatomical space... {0} -> {1}" .format (mask_ax [0 ],anat_ax [0 ]))
92+ logger .debug ("Aligning mask to anatomical space... {0} -> {1}" .format (
93+ mask_ax [0 ], anat_ax [0 ]))
11394 mask = nb .orientations .flip_axis (mask , 0 )
11495 if anat_ax [1 ] != 'P' :
11596 # flip anatspace
@@ -130,55 +111,54 @@ def deface(anat_filename, mask_filename, defaced_filename, buff=10):
130111 logger .debug ("Aligning mask to +y -> S" )
131112 mask = nb .orientations .flip_axis (mask , 2 )
132113
133-
134114 edgemask = edge_mask (mask )
135115 low = convex_hull (edgemask )
136- slope = (low [1 ][0 ]- low [1 ][1 ]) / (low [0 ][0 ]- low [0 ][1 ])
116+ slope = (low [1 ][0 ] - low [1 ][1 ]) / (low [0 ][0 ] - low [0 ][1 ])
137117
138- yint = low [1 ][0 ]- (low [0 ][0 ]* slope ) - buff
139- ys = numpy .arange (0 ,mask .shape [2 ])* slope + yint
140- defaced_mask = numpy .ones (mask .shape ,dtype = 'uint8' )
118+ yint = low [1 ][0 ] - (low [0 ][0 ] * slope ) - buff
119+ ys = numpy .arange (0 , mask .shape [2 ]) * slope + yint
120+ defaced_mask = numpy .ones (mask .shape , dtype = 'uint8' )
141121
142- for x in range (0 , ys .size - 1 ):
122+ for x in range (0 , ys .size - 1 ):
143123 if ys [x ] < 0 :
144124 break
145125 else :
146- ymax = min (ys [x ],mask .shape [2 ])
147- defaced_mask [:,x ,:ymax ] = 0
148-
149- defaced_img = defaced_mask * anat
150-
151-
152- if (anat_flip [1 ]):
153- newimg = nb .orientations .flip_axis (defaced_img , 1 )
154- if (anat_flip [2 ]):
155- newimg = nb .orientations .flip_axis (defaced_img , 2 )
156- else :
157- newimg = defaced_img
158- new_anat = nb .Nifti1Image (newimg ,nii_anat .get_affine ())
159- nb .save (new_anat ,defaced_filename )
126+ ymax = min (ys [x ], mask .shape [2 ])
127+ defaced_mask [:, x , :ymax ] = 0
128+
129+ defaced_img = defaced_mask * anat
130+
131+ newimg = defaced_img
132+ if anat_flip [1 ]:
133+ newimg = nb .orientations .flip_axis (newimg , 1 )
134+ if anat_flip [2 ]:
135+ newimg = nb .orientations .flip_axis (newimg , 2 )
136+ new_anat = nb .Nifti1Image (newimg , nii_anat .affine , nii_anat .header .copy ())
137+ nb .save (new_anat , defaced_filename )
160138 logger .info ("Defaced file: {0}" .format (defaced_filename ))
161139
140+
162141if __name__ == '__main__' :
163142
164143 logger = logging .getLogger (__name__ )
165- #logging.basicConfig(filename="hull.log",level=logging.DEBUG)
144+ # logging.basicConfig(filename="hull.log",level=logging.DEBUG)
166145 logger .setLevel (logging .DEBUG )
167146 ch = logging .StreamHandler ()
168147 ch .setLevel (logging .DEBUG )
169148 logger .addHandler (ch )
170149
171150 if len (sys .argv ) < 4 :
172- logger .debug ("Usage: quickshear.py anat_file strip_file defaced_file [buffer]" )
151+ logger .debug (
152+ "Usage: quickshear.py anat_file strip_file defaced_file [buffer]" )
173153 sys .exit (- 1 )
174154 else :
175155 anatfile = sys .argv [1 ]
176156 stripfile = sys .argv [2 ]
177157 newfile = sys .argv [3 ]
178158 if len (sys .argv ) >= 5 :
179159 try :
180- buff = sys .argv [4 ]
160+ buff = int ( sys .argv [4 ])
181161 except :
182- raise NumberError
162+ raise ValueError
183163 deface (anatfile , stripfile , newfile , buff )
184164 deface (anatfile , stripfile , newfile )
0 commit comments