Skip to content

Commit 90590c5

Browse files
tashrifbillahtashrifbillah
authored
Merge pull request #17 from pnlbwh/fit-shm-ref
Fit shm ref
2 parents 7d3cca7 + 46d4576 commit 90590c5

4 files changed

Lines changed: 132 additions & 81 deletions

File tree

‎lib/harmonization.py‎

Lines changed: 45 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -163,13 +163,13 @@ def createTemplate(self):
163163
# createTemplate steps -----------------------------------------------------------------------------------------
164164

165165
# read image lists
166-
refImgs, refMasks= common_processing(self.ref_csv)
166+
refImgs, refMasks= common_processing(self.ref_unproc_csv)
167167
if not self.ref_csv.endswith('.modified'):
168168
self.ref_csv += '.modified'
169169
# debug: use the following line to omit processing again
170170
# refImgs, refMasks = read_caselist(self.ref_csv)
171171

172-
targetImgs, targetMasks= common_processing(self.target_csv)
172+
targetImgs, targetMasks= common_processing(self.tar_unproc_csv)
173173
if not self.target_csv.endswith('.modified'):
174174
self.target_csv += '.modified'
175175
# debug: use the following line to omit processing again
@@ -248,8 +248,8 @@ def createTemplate(self):
248248

249249
def harmonizeData(self):
250250

251-
from reconstSignal import reconst
252-
from preprocess import dti_harm
251+
from reconstSignal import reconst, approx
252+
from preprocess import dti_harm, preprocessing, common_processing
253253

254254
# check the templatePath
255255
if not exists(self.templatePath):
@@ -258,49 +258,59 @@ def harmonizeData(self):
258258
if not listdir(self.templatePath):
259259
raise ValueError(f'{self.templatePath} is empty')
260260

261+
262+
# fit spherical harmonics on reference site
263+
if self.debug and self.ref_csv:
264+
check_csv(self.ref_unproc_csv, self.force)
265+
refImgs, refMasks= read_caselist(self.ref_unproc_csv)
266+
267+
# reference data is not manipulated in multi-shell-dMRIharmonization i.e. bvalMapped, resampled, nor denoised
268+
# this block may be uncommented in a future design
269+
# res= []
270+
# pool = multiprocessing.Pool(self.N_proc)
271+
# for imgPath, maskPath in zip(refImgs, refMasks):
272+
# res.append(pool.apply_async(func=preprocessing, args=(imgPath, maskPath)))
273+
#
274+
# attributes = [r.get() for r in res]
275+
#
276+
# pool.close()
277+
# pool.join()
278+
#
279+
# for i in range(len(refImgs)):
280+
# refImgs[i] = attributes[i][0]
281+
# refMasks[i] = attributes[i][1]
282+
283+
pool = multiprocessing.Pool(self.N_proc)
284+
for imgPath, maskPath in zip(refImgs, refMasks):
285+
pool.apply_async(func= approx, args=(imgPath,maskPath,))
286+
287+
pool.close()
288+
pool.join()
289+
290+
291+
261292
# go through each file listed in csv, check their existence, create dti and harm directories
262293
check_csv(self.target_csv, self.force)
263-
264-
# target data is not manipulated in multi-shell-dMRIharmonization i.e. bvalMapped, resampled, nor denoised
265-
# this block may be uncommented in a future design
266-
# from preprocess import dti_harm
267-
# if self.debug:
268-
# # calcuate diffusion measures of target site before any processing so we are able to compare
269-
# # with the ones after harmonization
270-
# imgs, masks= read_caselist(self.tar_unproc_csv)
271-
# pool = multiprocessing.Pool(self.N_proc)
272-
# for imgPath, maskPath in zip(imgs, masks):
273-
# imgPath= convertedPath(imgPath)
274-
# maskPath= convertedPath(maskPath)
275-
# pool.apply_async(func= dti_harm, args= ((imgPath, maskPath, )))
276-
#
277-
# pool.close()
278-
# pool.join()
294+
targetImgs, targetMasks= common_processing(self.tar_unproc_csv)
279295

280296
# reconstSignal steps ------------------------------------------------------------------------------------------
281297

282298
# read target image list
283299
moving= pjoin(self.templatePath, f'Mean_{self.target}_FA_b{self.bshell_b}.nii.gz')
284-
imgs, masks= read_caselist(self.tar_unproc_csv)
285300

286-
fm= None
287301
if not self.target_csv.endswith('.modified'):
288302
self.target_csv += '.modified'
289-
fm = open(self.target_csv, 'w')
290303

291304

292305
self.harm_csv= self.target_csv+'.harmonized'
293306
fh= open(self.harm_csv, 'w')
294307
pool = multiprocessing.Pool(self.N_proc)
295308
res= []
296-
for imgPath, maskPath in zip(imgs, masks):
309+
for imgPath, maskPath in zip(targetImgs, targetMasks):
297310
res.append(pool.apply_async(func= reconst, args= (imgPath, maskPath, moving, self.templatePath,)))
298311

299312
for r in res:
300-
imgPath, maskPath, harmImg, harmMask= r.get()
301-
302-
if isinstance(fm, io.IOBase):
303-
fm.write(imgPath + ',' + maskPath + '\n')
313+
harmImg, harmMask= r.get()
304314
fh.write(harmImg + ',' + harmMask + '\n')
305315

306316

@@ -314,23 +324,18 @@ def harmonizeData(self):
314324
# res.append(reconst(imgPath, maskPath, moving, self.templatePath))
315325
#
316326
# for r in res:
317-
# imgPath, maskPath, harmImg, harmMask= r
318-
#
319-
# if isinstance(fm, io.IOBase):
320-
# fm.write(imgPath + ',' + maskPath + '\n')
327+
# harmImg, harmMask= r
321328
# fh.write(harmImg + ',' + harmMask + '\n')
322329

323330

324-
if isinstance(fm, io.IOBase):
325-
fm.close()
326331
fh.close()
327332

328333

329334
if self.debug:
330335
harmImgs, harmMasks= read_caselist(self.harm_csv)
331336
pool = multiprocessing.Pool(self.N_proc)
332337
for imgPath,maskPath in zip(harmImgs,harmMasks):
333-
pool.apply_async(func= dti_harm, args= (imgPath,maskPath))
338+
pool.apply_async(func= dti_harm, args= (imgPath,maskPath,))
334339
pool.close()
335340
pool.join()
336341

@@ -440,10 +445,10 @@ def main(self):
440445
if self.N_proc==-1:
441446
self.N_proc= N_CPU
442447

443-
if self.target_csv.endswith('.modified'):
444-
self.tar_unproc_csv= str(self.target_csv).split('.modified')[0]
445-
else:
446-
self.tar_unproc_csv= str(self.target_csv)
448+
449+
if self.ref_csv:
450+
self.ref_unproc_csv= self.ref_csv.strip('.modified')
451+
self.tar_unproc_csv= self.target_csv.strip('.modified')
447452

448453

449454
# check appropriateness of N_shm
@@ -525,3 +530,4 @@ def main(self):
525530
if __name__ == '__main__':
526531
pipeline.run()
527532

533+

‎lib/multi-shell-harmonization.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -196,6 +196,7 @@ def main(self):
196196

197197
check_call((' ').join([pjoin(SCRIPTDIR, 'harmonization.py'),
198198
'--tar_list', tarListOutPrefix + f'_b{int(bval)}.csv',
199+
f'--ref_list {refListOutPrefix}_b{int(bval)}.csv' if self.ref_csv else '',
199200
'--bshell_b', str(int(bval)),
200201
'--process'] + pipeline_vars), shell=True)
201202

@@ -216,7 +217,9 @@ def main(self):
216217
## join harmonized data
217218
if self.process:
218219
joinAllBshells(self.target_csv, ref_bvals_file, 'harmonized_', self.N_proc)
219-
220+
221+
if self.debug and self.ref_csv:
222+
joinAllBshells(self.ref_csv, ref_bvals_file, 'reconstructed_', self.N_proc)
220223

221224
if __name__== '__main__':
222225
multi_shell_pipeline.run()

‎lib/preprocess.py‎

Lines changed: 40 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -84,10 +84,6 @@ def preprocessing(imgPath, maskPath):
8484
copyfile(inPrefix + '.bval', outPrefix + '.bval')
8585

8686
maskPath= maskPath
87-
88-
if debug:
89-
dti_harm(outPrefix+'.nii.gz', maskPath)
90-
9187
imgPath= outPrefix+'.nii.gz'
9288

9389

@@ -105,10 +101,6 @@ def preprocessing(imgPath, maskPath):
105101
write_bvals(outPrefix + '.bval', bvals)
106102

107103
maskPath= maskPath
108-
109-
if debug:
110-
dti_harm(outPrefix+'.nii.gz', maskPath)
111-
112104
imgPath= outPrefix+'.nii.gz'
113105

114106
try:
@@ -130,9 +122,6 @@ def preprocessing(imgPath, maskPath):
130122
else:
131123
maskPath= maskPath.split('.nii')[0]+ '_resampled.nii.gz'
132124

133-
if debug:
134-
dti_harm(outPrefix+'.nii.gz', maskPath)
135-
136125
imgPath= outPrefix+'.nii.gz'
137126

138127

@@ -141,39 +130,53 @@ def preprocessing(imgPath, maskPath):
141130

142131

143132
def common_processing(caselist):
133+
144134
imgs, masks = read_caselist(caselist)
145135

146-
res=[]
136+
# compute dti_harm of unprocessed data
147137
pool = multiprocessing.Pool(N_proc)
148-
for imgPath,maskPath in zip(imgs,masks):
149-
res.append(pool.apply_async(func= preprocessing, args= (imgPath,maskPath)))
150-
151-
attributes= [r.get() for r in res]
152-
153-
154-
pool.close()
155-
pool.join()
156-
157-
158-
f = open(caselist + '.modified', 'w')
159-
for i in range(len(imgs)):
160-
imgs[i] = attributes[i][0]
161-
masks[i] = attributes[i][1]
162-
f.write(f'{imgs[i]},{masks[i]}\n')
163-
f.close()
164-
165-
166-
# the following imgs, masks is for diagnosing MemoryError i.e. computing rish w/o preprocessing
167-
# to diagnose, comment all the above and uncomment the following
168-
# imgs, masks = read_caselist(caselist+'.modified')
169-
170-
# experimentally found ncpu=4 to be memroy optimal
171-
pool = multiprocessing.Pool(4)
172138
for imgPath,maskPath in zip(imgs,masks):
173139
pool.apply_async(func= dti_harm, args= (imgPath,maskPath))
174-
175140
pool.close()
176141
pool.join()
142+
143+
144+
try:
145+
copyfile(caselist, caselist + '.modified')
146+
except SameFileError:
147+
pass
148+
149+
# data is not manipulated in multi-shell-dMRIharmonization i.e. bvalMapped, resampled, nor denoised
150+
# this block may be uncommented in a future design
151+
# preprocess data
152+
# res=[]
153+
# pool = multiprocessing.Pool(N_proc)
154+
# for imgPath,maskPath in zip(imgs,masks):
155+
# res.append(pool.apply_async(func= preprocessing, args= (imgPath,maskPath)))
156+
#
157+
# attributes= [r.get() for r in res]
158+
#
159+
# pool.close()
160+
# pool.join()
161+
#
162+
# f = open(caselist + '.modified', 'w')
163+
# for i in range(len(imgs)):
164+
# imgs[i] = attributes[i][0]
165+
# masks[i] = attributes[i][1]
166+
# f.close()
167+
#
168+
#
169+
# # compute dti_harm of preprocessed data
170+
# pool = multiprocessing.Pool(N_proc)
171+
# for imgPath,maskPath in zip(imgs,masks):
172+
# pool.apply_async(func= dti_harm, args= (imgPath,maskPath))
173+
# pool.close()
174+
# pool.join()
175+
#
176+
#
177+
# if debug:
178+
# #TODO compute dti_harm for all intermediate data _denoised, _denoised_bmapped, _bmapped
179+
# pass
177180

178181
return (imgs, masks)
179182

‎lib/reconstSignal.py‎

Lines changed: 43 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,45 @@ def findLargestConnectMask(img, mask):
105105
return mask
106106

107107

108+
def approx(imgPath, maskPath):
109+
110+
print(f'Fitting spherical harmonics on {imgPath} ...')
111+
112+
directory = dirname(imgPath)
113+
inPrefix = imgPath.split('.nii')[0]
114+
prefix = psplit(inPrefix)[-1]
115+
outPrefix = pjoin(directory, 'harm', prefix)
116+
117+
b0, shm_coeff, qb_model = rish(imgPath, maskPath, inPrefix, outPrefix, N_shm)
118+
B = qb_model.B
119+
120+
img= load(imgPath)
121+
hdr= img.header
122+
affine= img.affine
123+
124+
S_hat= np.dot(shm_coeff, B.T)
125+
# keep only upper half of the reconstructed signal
126+
S_hat= S_hat[..., :int(S_hat.shape[3]/2)]
127+
np.nan_to_num(S_hat).clip(min= 0., max= 1., out= S_hat)
128+
129+
# affine= templateAffine for all Scale_L{i}
130+
mappedFile= pjoin(directory, f'{prefix}_mapped_cs.nii.gz')
131+
save_nifti(mappedFile, S_hat, affine, hdr)
132+
133+
# un-normalize approximated data
134+
S_hat_dwi= applymask(S_hat, b0) # overriding applymask function with a nonbinary mask b0
135+
136+
# place b0s in proper indices
137+
S_hat_final= stack_b0(qb_model.gtab.b0s_mask, S_hat_dwi, b0)
138+
139+
# save approximated data
140+
harmImg= pjoin(directory, f'reconstructed_{prefix}.nii.gz')
141+
if force or not isfile(harmImg):
142+
save_nifti(harmImg, S_hat_final, affine, hdr)
143+
copyfile(inPrefix + '.bvec', harmImg.split('.nii')[0] + '.bvec')
144+
copyfile(inPrefix + '.bval', harmImg.split('.nii')[0] + '.bval')
145+
146+
108147
def ring_masking(directory, prefix, maskPath, shm_coeff, b0, qb_model, hdr):
109148

110149
B = qb_model.B
@@ -178,9 +217,7 @@ def ring_masking(directory, prefix, maskPath, shm_coeff, b0, qb_model, hdr):
178217

179218

180219
def reconst(imgPath, maskPath, moving, templatePath):
181-
182-
imgPath, maskPath = preprocessing(imgPath, maskPath)
183-
220+
184221
img = load(imgPath)
185222

186223
directory = dirname(imgPath)
@@ -208,7 +245,7 @@ def reconst(imgPath, maskPath, moving, templatePath):
208245
copyfile(inPrefix + '.bval', harmImg.split('.nii')[0] + '.bval')
209246

210247

211-
return (imgPath, maskPath, harmImg, harmMask)
248+
return (harmImg, harmMask)
212249

213250
def stack_b0(b0s_mask, dwi, b0):
214251

@@ -225,3 +262,5 @@ def stack_b0(b0s_mask, dwi, b0):
225262
j+=1
226263

227264
return np.moveaxis(np.array(S_hat_final), 0, -1)
265+
266+

0 commit comments

Comments
 (0)