temptools.py 18.2 KB
Newer Older
Christoph Arthofer's avatar
Christoph Arthofer committed
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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
import os
import shutil
import pandas as pd
import nibabel as nib
import numpy as np
import shlex
import subprocess
import sys
from fsl.wrappers import fslmaths,flirt,applyxfm,concatxfm,bet,fast,fslstats
from fsl.wrappers.fnirt import invwarp, applywarp, convertwarp
from file_tree import FileTree
# from fsl.utils.filetree import FileTree
from fsl.utils.fslsub import func_to_cmd
from operator import itemgetter
import tempfile

def correctBiasMidtransWrapper(aff_matrix_paths, temp_dir, ref_path, unbiasing_invmtx_path, unbiased_matrix_paths):
    """! Writes the nonlinear registration parameters for a given resolution level and modalities to a file readable by MMORF.

    @param aff_matrix_paths:              List of filepaths to affine transformations
    @param temp_dir:                      Output directory
    @param ref_path:                      Path to reference template
    @param unbiasing_invmtx_path:         Path to unbiasing matrix
    @param unbiased_matrix_paths:         List of filepaths to unbiased transformations

    """
    separate_path = os.path.join(temp_dir, 'T1_to_unbiased')
    command = 'midtrans -v --separate=' + separate_path + ' --template=' + ref_path + ' -o ' + unbiasing_invmtx_path + ' '
    count = 0
    for omat_path in aff_matrix_paths:
        if os.path.exists(omat_path):
            count += 1
            print(count, ' ', omat_path)
            command += omat_path + ' '

    stream = os.popen(command)
    output = stream.read()
    print(output)

    # Renaming matrices
    for i, sub_unbiasing_mat in enumerate(unbiased_matrix_paths):
        sub_unbiasing_mat_temp = os.path.join(temp_dir, 'T1_to_unbiased%04d.mat' % (i + 1))

        os.rename(sub_unbiasing_mat_temp, sub_unbiasing_mat)
        print(i, ' ', sub_unbiasing_mat_temp, ' renamed to ', sub_unbiasing_mat)

    print('T1 unbiasing matrices constructed!')


def soft_clamp(x, k):
    """! Piecewise function for soft intensity clamping of T1 images. Takes a single parameter k which defines the transition to the clamping part of the function.

    f(x) = 0                                  | x <= 0
    f(x) = x                                  | 0 < x <= k
    f(x) = 3k/4 + k/(2(1 + exp(-8(x - k)/k))) | x > k

    @param x:              Image as numpy array
    @param k:              Defines the transition to the clamping part of the function

    Date: 08/02/2021
    Author: Frederik J Lange
    Copyright: FMRIB 2021

    """

    return np.piecewise(x,
                        [x <= 0, (0 < x) & (x <= k), x > k],
                        [lambda x: 0, lambda x: x, lambda x: k / (2 * (1 + np.exp(-8 * (x - k) / k))) + 0.75 * k])


def clampImage(img_path, out_path):
    """! Performs preprocessing steps and clamping on an image.

    @param img_path:              Path to input image
    @param out_path:              Path to clamped output image

    """

    out_dir = os.path.split(out_path)[0]
    mask_path = os.path.splitext(os.path.splitext(os.path.basename(out_path))[0])[0] + '_brain.nii.gz'
    mask_path = os.path.join(out_dir, mask_path)
    bet(img_path, mask_path, robust=True)
    with tempfile.TemporaryDirectory(dir=out_dir) as tmpdirname:
        fast(mask_path, tmpdirname + '/fast', iter=0, N=True, g=True, v=False)
        wm_intensity_mean = fslstats(mask_path).k(tmpdirname + '/fast_seg_2').M.run()
        print('White matter mean intensity is: ', wm_intensity_mean)

    img_nib = nib.load(img_path)
    img_clamped_np = soft_clamp(img_nib.get_fdata(), wm_intensity_mean)
    img_clamped_nib = nib.Nifti1Image(img_clamped_np, affine=img_nib.affine, header=img_nib.header)
    img_clamped_nib.to_filename(out_path)


def averageImages(img_paths, out_path, norm_bool=False):
    """! Creates an average image from individual (non)normalised images.

    @param img_paths:             List of filepaths
    @param out_path:              Path to average output image
    @param norm_bool:             Normalise intensities of each image before averaging true or false

    """

    n_exist = 0
    n_imgs = len(img_paths)
    for i, img_path in enumerate(img_paths):
        if os.path.exists(img_path):
            n_exist += 1
            print(i, ' ', img_path)
            img_nib = nib.load(img_path)
            if norm_bool:
                img_nib = fslmaths(img_nib).inm(1000).run()
            if i == 0:
                sum_img = img_nib
            else:
                sum_img = fslmaths(sum_img).add(img_nib).run()
        else:
            print(i, ' ', img_path, ' does not exist!')

    if n_exist > 0:
        mean_img = fslmaths(sum_img).div(n_exist).run()
        mean_img.to_filename(out_path)

    assert n_exist == n_imgs, "Not all images available!"


def applyWarpWrapper(img_path, ref_path, warped_path, warp_path, interp='spline', norm_bool=False):
    """! Wrapper for FSL applywarp - applies a warp (deformation field) to an image.

    @param img_path:              Path to input image
    @param ref_path:              Path to reference image
    @param warped_path:           Path to warped output image
    @param warp_path:             Path to warp (deformation field)
    @param interp:                Interpolation method (same options as FSL applywarp)
    @param norm_bool:             Normalise intensities of each image before averaging true or false

    """

    print(img_path, warp_path)
    if os.path.exists(img_path):
        img_nib = nib.load(img_path)
        if norm_bool:
            img_nib = fslmaths(img_nib).inm(1000).run()
        print('applywarp(src=img_nib,ref=ref_path,out=warped_path,warp=warp_path,interp=interp)')
        applywarp(src=img_nib, ref=ref_path, out=warped_path, warp=warp_path, interp=interp)


# def submitJob_fsl_sub(name, log_dir, queue, wait_for=[], script=None, command=None, coprocessor_class=None, export_var=None,
#                   debug=False):
#     """! Wrapper for fslsub - submits a job to the cluster. This function can be easily extended to work with other workload managers.
#
#     @param name:                  Job name
#     @param log_dir:               Directory where output log-files will be saved
#     @param queue:                 Name of queue to submit the job to
#     @param wait_for:              List of IDs of jobs required to finish before running this job.
#     @param script:                Path to a shell script, which contains one command per line - commands will be submitted as an array job
#     @param command:               Alternatively a single command can be provided as a string - command will be submitted as single job
#     @param coprocessor_class:     Coprocessor class, if not None cuda will be selected
#     @param export_var:            Environment variables to be exported to the submission node
#     @param debug:                 If True, information about job will be written to output
#
#     @return  The job ID.
#     """
#
#     fsl_sub.submit()
#
#     cmd = 'fsl_sub'
#     if wait_for and any(job != '' for job in wait_for):
#         cmd += ' -j '
#         for j, job in enumerate(wait_for):
#             if job != '':
#                 cmd += job.replace("\n", "")
#                 if j < len(wait_for) - 1:
#                     cmd += ','
#
#     cmd += ' -N ' + name + \
#            ' -l ' + log_dir + \
#            ' -q ' + queue
#
#     if coprocessor_class is not None :
#         cmd += ' --coprocessor cuda'
#
#     if export_var is not None :
#         cmd += ' --export ' + export_var
#
#     if debug:
#         cmd += ' --debug'
#
#     if script is not None and os.path.exists(script):
#         cmd += ' -t ' + script
#     elif command is not None :
#         cmd += shlex.split(command)
#
#     # stream = os.popen(cmd)
#     # job_id = stream.read()
#
#     try:
#         result = subprocess.run(command, capture_output, text=True, check=True)
#     except subprocess.CalledProcessError as e:
#         print(str(e), file=sys.stderr)
#         return None
#
#     job_id = result.stdout.strip()
#
#     return job_id


def submitJob(name, log_dir, queue, wait_for=[], script=None, command=None, coprocessor_class=None, export_var=None,
                  debug=False):
    """! Wrapper for fslsub - submits a job to the cluster. This function can be easily extended to work with other workload managers.

    @param name:                  Job name
    @param log_dir:               Directory where output log-files will be saved
    @param queue:                 Name of queue to submit the job to
    @param wait_for:              List of IDs of jobs required to finish before running this job.
    @param script:                Path to a shell script, which contains one command per line - commands will be submitted as an array job
    @param command:               Alternatively a single command can be provided as a string - command will be submitted as single job
    @param coprocessor_class:     Coprocessor class, if not None cuda will be selected
    @param export_var:            Environment variables to be exported to the submission node
    @param debug:                 If True, information about job will be written to output

    @return  The job ID.
    """
    cmd = 'fsl_sub'
    if wait_for and any(job != '' for job in wait_for):
        cmd += ' -j '
        for j, job in enumerate(wait_for):
            if job != '':
                cmd += job.replace("\n", "")
                if j < len(wait_for) - 1:
                    cmd += ','

    cmd += ' -N ' + name + \
           ' -l ' + log_dir + \
           ' -q ' + queue

    if coprocessor_class is not None :
        cmd += ' --coprocessor cuda'

    if export_var is not None :
        cmd += ' --export ' + export_var

    if debug:
        cmd += ' --debug'

    if script is not None and os.path.exists(script):
        cmd += ' -t ' + script
    elif command is not None :
        cmd += ' ' + command + ' '
        # cmd += ' "' + command + '"'

    # stream = os.popen(cmd)
    # job_id = stream.read()

    print(cmd)

    try:
        result = subprocess.run(shlex.split(cmd), capture_output=True, text=True, check=True)
    except subprocess.CalledProcessError as e:
        print(str(e), file=sys.stderr)
        return None

    job_id = result.stdout.strip()

    return job_id


def RMSdifference(img1_path, img2_path, mask1_path=None, mask2_path=None, rms_path=None):
    """! Calculates the difference between two images or warps as the root mean squared (RMS)

    @param img1_path:                 Path to first image or deformation field
    @param img2_path:                 Path to second image or deformation field
    @param mask1_path:                Path to mask for first image
    @param mask2_path:                Path to mask for second image
    @param rms_path:                  Path to output text file that RMS is written to

    """

    img1_arr = nib.load(img1_path).get_fdata()
    img2_arr = nib.load(img2_path).get_fdata()

    if mask1_path is not None and mask2_path is not None:
        mask1_arr = nib.load(mask1_path).get_fdata()
        mask2_arr = nib.load(mask2_path).get_fdata()

        if len(img1_arr.shape) > 3:
            n_dim = img1_arr.shape[-1]
            img1_mask_stack_arr = np.stack((mask1_arr,) * n_dim, -1)
            n_dim = img2_arr.shape[-1]
            img2_mask_stack_arr = np.stack((mask2_arr,) * n_dim, -1)
        else:
            img1_mask_stack_arr = mask1_arr
            img2_mask_stack_arr = mask2_arr

        img_mask_stack_arr = np.logical_or(img1_mask_stack_arr, img2_mask_stack_arr)
        img1_masked_arr = img1_arr[img_mask_stack_arr > 0]
        img2_masked_arr = img2_arr[img_mask_stack_arr > 0]

        diff_img = img1_masked_arr - img2_masked_arr
    else:
        diff_img = img1_arr - img2_arr

    dim = np.prod(diff_img.shape)
    rms = np.sqrt((diff_img ** 2).sum() / dim)

    print('RMS difference between {} and {}: {}'.format(img1_path, img2_path, rms))
    if rms_path is not None:
        with open(rms_path, 'w+') as f:
            f.write('{}'.format(rms))


def RMSstandardDeviation(img_paths, mean_img_path, mask_path, sd_img_out_path=None, rms_out_path=None):
    """! Calculates the standard deviation of images as the root mean squared (RMS) (== coefficient of variation)

    @param img_paths:                     List of paths to images
    @param mean_img_path:                 Path to average image
    @param mask_path:                     Path to mask
    @param sd_img_out_path:               Path to standard deviation output image
    @param rms_out_path:                  Path to output text file that RMS is written to

    """

    mean_img_nib = nib.load(mean_img_path)

    for i, path in enumerate(img_paths):
        print(i, ' ', path)
        diff_img = fslmaths(path).inm(1000).sub(mean_img_nib).run()
        if i == 0:
            diffsum_img = fslmaths(diff_img).mul(diff_img).run()
        else:
            diffsum_img = fslmaths(diff_img).mul(diff_img).add(diffsum_img).run()

    stdtemp_img = fslmaths(diffsum_img).div(len(img_paths)).run()
    stdtemp_img_np = np.sqrt(stdtemp_img.get_fdata())
    if sd_img_out_path is not None:
        std_img = nib.Nifti1Image(stdtemp_img_np, affine=stdtemp_img.affine, header=stdtemp_img.header)
        std_img.to_filename(sd_img_out_path)

    mask_np = nib.load(mask_path).get_fdata()
    mean_img_np = mean_img_nib.get_fdata()

    stdtemp_img_masked_np = stdtemp_img_np[mask_np > 0]
    mean_img_masked_np = mean_img_np[mask_np > 0]

    cv = stdtemp_img_masked_np / mean_img_masked_np  # coefficient of variation
    dim = np.prod(cv.shape)
    rms = np.sqrt((cv ** 2).sum() / dim)

    print('RMS standard deviation: {}'.format(rms))
    if rms_out_path is not None:
        with open(rms_out_path, 'w+') as f:
            f.write('{}'.format(rms))


def mmorfWrapper(mmorf_run_cmd, config_path, img_warp_space,
                 img_ref_scalar, img_mov_scalar, aff_ref_scalar, aff_mov_scalar,
                 mask_ref_scalar, mask_mov_scalar,
                 img_ref_tensor, img_mov_tensor, aff_ref_tensor, aff_mov_tensor,
                 mask_ref_tensor, mask_mov_tensor,
                 warp_out, jac_det_out, bias_out):
    """! Wrapper function for running MMORF.

    @param mmorf_run_cmd:                     Singularity command to run MMORF
    @param config_path:                       Path to config file with fixed parameters
    @param img_warp_space:                    Path to image defining the space in which the warp field will be calculated
    @param img_ref_scalar:                    List of paths to scalar reference images
    @param img_mov_scalar:                    List of paths to scalar moving images
    @param aff_ref_scalar:                    List of paths to affine transformations for scalar reference images
    @param aff_mov_scalar:                    List of paths to affine transformations for scalar moving images
    @param mask_ref_scalar:                   List of paths to masks in reference image spaces
    @param mask_mov_scalar:                   List of paths to masks in moving image spaces
    @param img_ref_tensor:                    List of paths to reference tensors
    @param img_mov_tensor:                    List of paths to moving tensors
    @param aff_ref_tensor:                    List of paths to affine transformations for reference tensors
    @param aff_mov_tensor:                    List of paths to affine transformations for moving tensors
    @param mask_ref_tensor:                   List of paths to masks in reference tensor spaces
    @param mask_mov_tensor:                   List of paths to masks in moving tensor spaces
    @param warp_out:                          Path to output warp field
    @param jac_det_out:                       Path to output Jacobian determinant of final warp field
    @param bias_out:                          Path to output bias field for scalar image pairs

    @return  The command as a string and a dictionary of environment variables.

    """

    export_var = []
    cmd = mmorf_run_cmd
    cmd += ' --config ' + config_path
    split = os.path.split(config_path)
    export_var.append(split[0])
    cmd += ' --img_warp_space ' + img_warp_space
    split = os.path.split(img_warp_space)
    export_var.append(split[0])
    for path in img_ref_scalar:
        cmd += ' --img_ref_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in img_mov_scalar:
        cmd += ' --img_mov_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in aff_ref_scalar:
        cmd += ' --aff_ref_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in aff_mov_scalar:
        cmd += ' --aff_mov_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in mask_ref_scalar:
        cmd += ' --mask_ref_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in mask_mov_scalar:
        cmd += ' --mask_mov_scalar ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in img_ref_tensor:
        cmd += ' --img_ref_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in img_mov_tensor:
        cmd += ' --img_mov_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in aff_ref_tensor:
        cmd += ' --aff_ref_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in aff_mov_tensor:
        cmd += ' --aff_mov_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in mask_ref_tensor:
        cmd += ' --mask_ref_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    for path in mask_mov_tensor:
        cmd += ' --mask_mov_tensor ' + path
        split = os.path.split(path)
        export_var.append(split[0])
    cmd += ' --warp_out ' + warp_out
    split = os.path.split(warp_out)
    export_var.append(split[0])
    cmd += ' --jac_det_out ' + jac_det_out
    split = os.path.split(jac_det_out)
    export_var.append(split[0])
    cmd += ' --bias_out ' + bias_out
    split = os.path.split(bias_out)
    export_var.append(split[0])

    cmd += '\n'

    export_var = list(filter(None, list(set(export_var))))
    export_var = {'SINGULARITY_BIND': export_var}

    return cmd, export_var