-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnifti_python_gadgetron_svrtk.py
445 lines (369 loc) · 19.2 KB
/
nifti_python_gadgetron_svrtk.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
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
import numpy as np
import gadgetron
import ismrmrd
import logging
import time
import io
import os
from datetime import datetime
from scipy.ndimage import zoom
import skimage
from scipy import ndimage
from skimage.measure import label, regionprops
import torch
import monai
from monai.inferers import sliding_window_inference
from monai.networks.nets import DenseNet121, UNet, AttentionUnet
import warnings
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning, module="monai.utils.module")
warnings.filterwarnings("ignore")
warnings.simplefilter("ignore")
from ismrmrd.meta import Meta
import itertools
import ctypes
import copy
import io
import warnings
import scipy.ndimage as ndi
import time
import shutil
import SimpleITK as sitk
warnings.simplefilter('default')
from ismrmrd.acquisition import Acquisition
from ismrmrd.flags import FlagsMixin
from ismrmrd.equality import EqualityMixin
from ismrmrd.constants import *
import matplotlib.image
import matplotlib.pyplot as plt
from matplotlib.widgets import RectangleSelector
import sys
# from python_version import extract_ismrmrd_parameters_from_headers as param, flip_image as fi, set_nii_hdr as tools
import nibabel as nib
import src.utils as utils
from src.utils import ArgumentsTrainTestLocalisation, plot_losses_train
from src import networks as md
from numpy.fft import fftshift, ifftshift, fftn, ifftn
import subprocess
def get_first_index_of_non_empty_header(header):
# if the data is under-sampled, the corresponding acquisition Header will be filled with 0
# in order to catch valuable information, we need to catch a non-empty header
# using the following lines
print(np.shape(header))
dims = np.shape(header)
for ii in range(0, dims[0]):
# print(header[ii].scan_counter)
if header[ii].scan_counter > 0:
break
print(ii)
return ii
def append_new_line(file_name, text_to_append):
"""Append given text as a new line at the end of file"""
with open(file_name, "a+") as file_object:
# Move read cursor to the start of file.
file_object.seek(0)
# If file is not empty then append '\n'
data = file_object.read(100)
if len(data) > 0:
file_object.write("\n")
# Append text at the end of file
file_object.write(text_to_append)
def send_reconstructed_images(connection, data_array, acq_header):
# the function creates a new ImageHeader for each 4D dataset [RO,E1,E2,CHA]
# copy information from the acquisitionHeader
# fill additional fields
# and send the reconstructed image and ImageHeader to the next gadget
# get header info
hdr = connection.header
enc = hdr.encoding[0]
if enc.encodingLimits.slice is not None:
nslices = enc.encodingLimits.slice.maximum + 1
else:
nslices = 1
if enc.encodingLimits.repetition is not None:
nreps = enc.encodingLimits.repetition.maximum + 1
else:
nreps = 1
ncoils = 1
dims = data_array.shape
base_header = ismrmrd.ImageHeader()
base_header.version = acq_header.version
ndims_image = (dims[0], dims[1], dims[2], dims[3])
base_header.channels = ncoils # The coils have already been combined
base_header.matrix_size = (data_array.shape[0], data_array.shape[1], data_array.shape[2])
base_header.position = acq_header.position
base_header.read_dir = acq_header.read_dir
base_header.phase_dir = acq_header.phase_dir
base_header.slice_dir = acq_header.slice_dir
base_header.patient_table_position = acq_header.patient_table_position
base_header.acquisition_time_stamp = acq_header.acquisition_time_stamp
base_header.image_index = 0
base_header.image_series_index = 0
base_header.data_type = ismrmrd.DATATYPE_CXFLOAT
base_header.image_type = ismrmrd.IMTYPE_COMPLEX
base_header.repetition = acq_header.repetition
ninstances = nslices * nreps
r = np.zeros((dims[0], dims[1], ninstances))
r = data_array
base_header.image_type = ismrmrd.IMTYPE_COMPLEX
image_array = ismrmrd.image.Image.from_array(r, headers=base_header)
print("..................................................................................")
logging.info("Sending reconstructed slice to the scanner...")
connection.send(image_array)
logging.info("Sent!")
print("..................................................................................")
def figstring(name):
date_path = datetime.today().strftime("%Y-%m-%d")
timestamp = f"{datetime.today().strftime('%H-%M-%S')}"
if not os.path.isdir("/home/sn21/data/t2-stacks/" + date_path):
os.mkdir("/home/sn21/data/t2-stacks/" + date_path)
final_str_nii = f"{date_path}{os.sep}{timestamp}-{name}"
final_str_no_ext = f"{date_path}{os.sep}{timestamp}-{name}"
return final_str_nii
def IsmrmrdToNiftiGadget(connection):
logging.info("Initializing data processing in Python...")
# start = time.time()
# get header info
hdr = connection.header
enc = hdr.encoding[0]
if enc.encodingLimits.slice is not None:
nslices = enc.encodingLimits.slice.maximum + 1
else:
nslices = 1
if enc.encodingLimits.repetition is not None:
nreps = enc.encodingLimits.repetition.maximum + 1
else:
nreps = 1
if enc.encodingLimits.contrast is not None:
ncontrasts = enc.encodingLimits.contrast.maximum + 1
else:
ncontrasts = 1
print("Number of echoes =", ncontrasts)
ncoils = 1
# Matrix size
eNx = int(enc.encodedSpace.matrixSize.x * 2)
eNy = enc.encodedSpace.matrixSize.y
eNz = enc.encodedSpace.matrixSize.z
print("eNx = ", eNx, "eNy = ", eNy, "eNz = ", eNz)
# Initialise a storage array
# eNy = enc.encodingLimits.kspace_encoding_step_1.maximum + 1
ninstances = nslices * nreps
# print("Number of instances ", ninstances)
im = np.zeros((eNx, eNy, nslices), dtype=np.complex64)
print("Image Shape ", im.shape)
# Get today's date
date_path = datetime.today().strftime("%Y-%m-%d")
# Path to the main DICOMs folder
dicoms_path = f"/home/sn21/data/t2-stacks/{date_path}/dicoms"
# Check if the main DICOMs folder exists
if not os.path.exists(dicoms_path):
logging.error("DICOMs folder not found.")
return
# This can all be uncommented out once the entire automated pipeline is working!
#
# # List all subdirectories in the main DICOMs folder
# subfolders = [folder for folder in os.listdir(dicoms_path) if os.path.isdir(os.path.join(dicoms_path, folder))]
#
# # Iterate over each subdirectory
# for subfolder in subfolders:
# # Path to the current DICOMs subfolder
# subfolder_path = os.path.join(dicoms_path, subfolder)
#
# # Check if the current subfolder is named correctly (dicoms_$number)
# if not subfolder.startswith("dicoms_"):
# logging.warning(f"Ignoring folder '{subfolder}' as it does not follow the naming convention.")
# continue
#
# # Check if NIfTI files already exist in the subfolder
# nifti_files_exist = any(file.lower().endswith('.nii') for file in os.listdir(subfolder_path))
# if nifti_files_exist:
# logging.warning(f"NIfTI files already exist in '{subfolder}'. Skipping DICOM to NIfTI conversion.")
# continue
#
# print("subfolder", subfolder)
#
# # Construct the command to convert DICOMs to NIfTI
# dcm2nii_command = f"dcm2nii {subfolder_path}"
# subprocess.call(dcm2nii_command, shell=True)
#
# # List all .nii.gz files in the directory
# nifti_files = [file for file in os.listdir(subfolder_path) if file.endswith('.nii.gz')]
#
# for nifti_file in nifti_files:
# nifti_file_path = os.path.join(subfolder_path, nifti_file)
#
# # Read the NIfTI file
# image = sitk.ReadImage(nifti_file_path)
# data = sitk.GetArrayFromImage(image)
#
# # Determine the midpoint slice
# nslices = data.shape[0]
# if nslices % 2 != 0:
# mid = int(nslices / 2) + 1
# else:
# mid = int(nslices / 2)
#
# # Split the data into two halves along the slice dimension
# data_1 = data[:mid, :, :]
# data_2 = data[mid:, :, :]
#
# # Interleave the slices from the two halves
# data_interleaved = np.zeros_like(data)
# data_interleaved[::2, :, :] = data_1
# data_interleaved[1::2, :, :] = data_2
#
# # Save the interleaved data as a new NIfTI file with the same name
# output_file_path = nifti_file_path.replace(".nii.gz", ".nii.gz")
# interleaved_image = sitk.GetImageFromArray(data_interleaved)
# interleaved_image.CopyInformation(image)
#
# # Overwrite the slice thickness to be 4.5 mm
# interleaved_image.SetSpacing((image.GetSpacing()[0], image.GetSpacing()[1], 4.5))
# sitk.WriteImage(interleaved_image, output_file_path)
# print("Interleaved NIfTI file saved as:", output_file_path)
# sitk.WriteImage(interleaved_image, output_file_path)
#
# print("Interleaved NIfTI file saved as:", output_file_path)
# #
# # for subfolder in subfolders:
# # subfolder_path = os.path.join(dicoms_path, subfolder)
# # nifti_files = [file for file in os.listdir(subfolder_path) if file.endswith('_interleaved.nii.gz')
# # and not file.startswith('o')]
# # for index, nifti_file in enumerate(nifti_files):
# # source_path = os.path.join(subfolder_path, nifti_file)
# # filename = os.path.splitext(nifti_file)[0]
# # filename = os.path.splitext(filename)[0]
# # new_filename = f"{filename}_{subfolder}_{index + 1}.nii.gz"
# # destination_path = os.path.join(dicoms_path, new_filename)
# # shutil.copy(source_path, destination_path)
# #
# # logging.info("Copied all .nii.gz files to the main DICOMs folder.")
#
# for subfolder in subfolders:
# subfolder_path = os.path.join(dicoms_path, subfolder)
# nifti_files = [file for file in os.listdir(subfolder_path) if
# file.endswith('0.nii.gz') and not file.startswith('o')]
# for index, nifti_file in enumerate(nifti_files):
# source_path = os.path.join(subfolder_path, nifti_file)
# filename = os.path.splitext(nifti_file)[0]
# filename = os.path.splitext(filename)[0]
# new_filename = f"{filename}_{subfolder}_{index + 1}.nii.gz"
# destination_path = os.path.join(dicoms_path, new_filename)
# shutil.move(source_path, destination_path)
#
# # Remove nii.gz files starting with 'o'
# unwanted_files = [file for file in os.listdir(subfolder_path) if file.endswith('.nii.gz')]
# for file in unwanted_files:
# file_path = os.path.join(subfolder_path, file)
# os.remove(file_path)
#
# logging.info("Moved appropriate .nii.gz files to the main DICOMs folder and removed files starting with 'o'.")
for acquisition in connection:
# print(acquisition)
# imag = np.abs(acquisition.data.transpose(3, 2, 1, 0))
# print("Slice Dimensions ", imag.shape)
# ndim = imag.shape
# print("ndim ", ndim)
# Get crop image, flip and rotate to match with true Nifti image
# img = imag[:, :, :, 0]
# Stuff into the buffer
slice = acquisition.slice
repetition = acquisition.repetition
contrast = acquisition.contrast
# print("Repetition ", repetition, "Slice ", slice, "Contrast ", contrast)
logging.info("Initializing data processing in Python...")
# # Copy the generated .nii.gz files from each subfolder to the main DICOMs folder
# for subfolder in subfolders:
# subfolder_path = os.path.join(dicoms_path, subfolder)
# nifti_files = [file for file in os.listdir(subfolder_path) if file.endswith('.nii.gz')]
# for index, nifti_file in enumerate(nifti_files):
# source_path = os.path.join(subfolder_path, nifti_file)
# # Get the filename without the extension
# filename = os.path.splitext(nifti_file)[0]
# filename = os.path.splitext(filename)[0]
# # Create the new filename with the index
# # Construct the new filename with the subfolder name and index
# new_filename = f"{filename}_{subfolder}_{index + 1}.nii.gz"
# destination_path = os.path.join(dicoms_path, new_filename)
# shutil.copy(source_path, destination_path)
#
# logging.info("Copied all .nii.gz files to the main DICOMs folder.")
# if the whole stack of slices has been acquired >> apply network to the entire 3D volume
if slice == nslices - 1:
print("..................................................................................")
# List all files in the directory
date_path = datetime.today().strftime("%Y-%m-%d")
directory_path = '/home/sn21/data/t2-stacks/' + date_path
files = [file for file in os.listdir(directory_path) if file.endswith('.nii.gz')]
# Count the number of files
num_files = len(files)
print(f"There are {num_files} files in the directory.")
print("Launching docker now...")
# Set the DISPLAY and XAUTHORITY environment variables
os.environ['DISPLAY'] = ':1' # Replace with your X11 display, e.g., ':1.0'
os.environ['XAUTHORITY'] = '/home/sn21/.Xauthority'
# Record the start time
start_time = time.time()
# Construct the command for docker
date_path = datetime.today().strftime("%Y-%m-%d")
# command = f'''docker run --rm --mount type=bind,source=/home/sn21/data/t2-stacks,target=/home/data \
# fetalsvrtk/svrtk:general_auto_amd sh -c 'bash
# /home/auto-proc-svrtk/scripts/auto-brain-055t-reconstruction.sh \
# /home/data/{date_path}/dicoms /home/data/{date_path}-result 1 4.5 1.0 1 ;
# chmod 1777 -R /home/data/{date_path}-result;
# bash /home/auto-proc-svrtk/scripts/auto-body-055t-reconstruction.sh \
# /home/data/{date_path}/dicoms /home/data/{date_path}-result-body 1 4.5 1.0 1 ;
# chmod 1777 -R /home/data/{date_path}-result-body; ' '''
# ' '''
command = f'''docker run --rm \
--mount type=bind,source=/home/sn21/data/t2-stacks,target=/home/data \
fetalsvrtk/svrtk:general_auto_amd sh -c 'bash /home/auto-proc-svrtk/scripts/auto-brain-055t-reconstruction.sh \
/home/data/{date_path}/dicoms /home/data/{date_path}-result 1 4.5 1.0 1 ; \
chmod 1777 -R /home/data/{date_path}-result; \
/bin/MIRTK/build/lib/tools/pad-3d /home/data/{date_path}-result/reo-SVR-output-brain.nii.gz /home/ref.nii.gz 160 1 ; \
/bin/MIRTK/build/lib/tools/edit-image /home/ref.nii.gz /home/ref.nii.gz -dx 1 -dy 1 -dz 1 ; \
/bin/MIRTK/build/lib/tools/transform-image /home/data/{date_path}-result/reo-SVR-output-brain.nii.gz \
/home/data/{date_path}-result/grid-reo-SVR-output-brain.nii.gz -target /home/ref.nii.gz -interp BSpline ; \
/bin/MIRTK/build/lib/tools/nan /home/data/{date_path}-result/grid-reo-SVR-output-brain.nii.gz 1000000 ; \
/bin/MIRTK/build/lib/tools/convert-image /home/data/{date_path}-result/grid-reo-SVR-output-brain.nii.gz \
/home/data/{date_path}-result/grid-reo-SVR-output-brain.nii.gz -short ; \
chmod 1777 /home/data/{date_path}-result/grid-reo-SVR-output-brain.nii.gz ; \
bash /home/auto-proc-svrtk/scripts/auto-body-055t-reconstruction.sh \
/home/data/{date_path}/dicoms /home/data/{date_path}-result-body 1 4.5 1.0 1 ; \
chmod 1777 -R /home/data/{date_path}-result-body; \
/bin/MIRTK/build/lib/tools/pad-3d /home/data/{date_path}-result-body/reo-DSVR-output-body.nii.gz /home/ref.nii.gz 256 1 ; \
/bin/MIRTK/build/lib/tools/edit-image /home/ref.nii.gz /home/ref.nii.gz -dx 1 -dy 1 -dz 1 ; \
/bin/MIRTK/build/lib/tools/transform-image /home/data/{date_path}-result-body/reo-DSVR-output-body.nii.gz \
/home/data/{date_path}-result-body/grid-reo-DSVR-output-body.nii.gz -target /home/ref.nii.gz -interp BSpline ; \
/bin/MIRTK/build/lib/tools/nan /home/data/{date_path}-result-body/grid-reo-DSVR-output-body.nii.gz 1000000 ; \
/bin/MIRTK/build/lib/tools/convert-image /home/data/{date_path}-result-body/grid-reo-DSVR-output-body.nii.gz \
/home/data/{date_path}-result-body/grid-reo-DSVR-output-body.nii.gz -short ; \
chmod 1777 /home/data/{date_path}-result-body/grid-reo-DSVR-output-body.nii.gz ; \
' '''
subprocess.Popen(['gnome-terminal', '--', 'bash', '-c', command])
# os.system(command)
# Record the end time
end_time = time.time()
# Calculate the elapsed time
elapsed_time = end_time - start_time
print(f"Execution time: {elapsed_time} seconds")
# Get today's date
date = datetime.now().strftime("%Y-%m-%d")
# Define the file path
file_path = f"/home/sn21/data/t2-stacks/{date}/execution_time.txt"
# Save the execution time to a text file
with open(file_path, "w") as file:
file.write(f"Execution time: {elapsed_time} seconds")
print(f"Execution time: {elapsed_time} seconds")
print(f"Execution time saved to: {file_path}")
# print("Docker command is running in detached mode.")
# command = f'''docker run --rm --mount type=bind,source=/home/sn21/data/t2-stacks,target=/home/data
# fetalsvrtk/svrtk:auto-2.20 sh -c 'bash /home/auto-proc-svrtk/auto-brain-reconstruction.sh /home/data/{
# date_path} /home/data/{date_path}-result 1 4.5 1.0 1 ;' '''
#
# # Use subprocess to run the Docker command in the current terminal
# subprocess.call(command, shell=True)
# # Use the gnome-terminal to run the command in a detached mode
# # subprocess.Popen(['gnome-terminal', '--', 'bash', '-c', docker_command])
# subprocess.Popen(['gnome-terminal', '--', 'bash', '-c', command])
# time.sleep(5)