Romeving file not relevant to 4Dflow visualisation

# #! /usr/bin/env python3

# Some references:

# - pydicom anonymization tutorial:

#   https://pydicom.github.io/pydicom/stable/auto_examples/metadata_processing/plot_anonymize.html

# - other script:

#   https://github.com/janpipek/anonymize_dicom

import os

import sys

import traceback

from PyQt5 import QtCore, QtWidgets

import pydicom

TAGS_REMOVE = (

    "PatientBirthDate",

    "InstitutionAddress",

    "OtherPatientIDs",

    "PerformingPhysicianName",

    "OperatorsName",

    "ReferringPhysicianName",

    "StudyDate",

    "AcquisitionDate",

    "ContentDate",

    "AcquisitionDateTime",

    "StudyTime",

    "ContentTime",

    "SeriesTime",

    "AcquisitionTime",

    "AccessionNumber",

    "StationName",

    "StudyID",

    "RequestingService",

    "RequestedProcedureDescription",

    # doplneno

    "PatientSex",

    #"InstanceCreatorUID",

    "StudyDescription",

    #"SeriesDescription",

    "InstitutionalDepartmentName",

    "PhysiciansOfRecord",

    "NameOfPhysicianReadingStudy",

    "AdmittingDiagnosesDescription",

    "DerivationDescription",

    "PatientBirthTime",

    "OtherPatientNames",

    "PatientAge",

    "PatientSize",

    "PatientWeight",

    "EthnicGroup",

    "Occupation",

    "AdditionalPatientHistory",

    "PatientComments",

    "DeviceSerialNumber",

    "ProtocolName",

    "ImageComments",

    #"RequestAttributesSequence"

    # doplneno po posouzeni vysledku

    "SeriesDate",

    "InstanceCreationDate",

    "InstanceCreationTime"

)

TAGS_SUBSTITUTE = {

    "PatientName": "Patient_14",

    "PatientID": "G318_patient14_bodyPart1_visit1",

    "InstitutionName": "deleted",

    # doplneno

    #"StudyInstanceUID": "deleted",

    #"SeriesInstanceUID": "deleted",

    #"FrameOfReferenceUID": "deleted",

    #"SynchronizationFrameOfReferenceUID": "deleted",

    #"SOPInstanceUID": "deleted",

    #"ReferencedSOPInstanceUID": "deleted",

    #"UID": "deleted",

    #"ContentSequence": "deleted", # nefunguje na 99-PHOENIX, vyresit

    #"StorageMediaFileSetUID": "deleted",

    # doplneno, neznamy Type (1, 1C, 2, 2C, 3?)

    "ReferringPhysicianAddress": "deleted",

    "ReferringPhysicianTelephoneNumber": "deleted",

    "MedicalRecordLocator": "deleted",

    #"ReferencedFrameOfReferenceUID": "deleted",

    #"RelatedFrameOfReferenceUID": "deleted"

    # doplneno po posouzeni vysledku

    "RequestingPhysician": "deleted"

}

def is_dicom_file(path):

    """Fast way to check whether file is DICOM.

    :rtype: bool

    """

    if not os.path.isfile(path):

        return False

    try:

        with open(path, "rb") as f:

            return f.read(132).decode("ASCII")[-4:] == "DICM"

    except:

        return False

def dicom_files_in_dir(directory="."):

    """Full paths of all DICOM files in the directory.

    :rtype: list (str)

    """

    directory = os.path.expanduser(directory)

    candidates = [os.path.join(directory, f) for f in sorted(os.listdir(directory))]

    return [f for f in candidates if is_dicom_file(f)]

def anonymize(input_filename, output_filename, *, keep_private_tags=False):

    assert not os.path.exists(output_filename)

    ds = pydicom.dcmread(input_filename)

    def _delete_field(field):

        if field not in ds:

            return

        delattr(ds, field)

    def _update_field(field, value):

        if field not in ds:

            return

        setattr(ds, field, value)

    for field in TAGS_REMOVE:

        _delete_field(field)

    for field, value in TAGS_SUBSTITUTE.items():

        _update_field(field, value)

    # Private DICOM tags are not private metadata, see:

    # http://dicom.nema.org/dicom/2013/output/chtml/part05/sect_7.8.html

    if not keep_private_tags:

        ds.remove_private_tags()

    # create parent directory

    basedir = os.path.dirname(output_filename)

    os.makedirs(basedir, exist_ok=True)

    ds.save_as(output_filename)

def anonymize_tree(input_directory, output_directory, *, keep_private_tags=False, update_signal=None):

    input_directory = os.path.realpath(input_directory)

    output_directory = os.path.realpath(output_directory)

    def anonymize_dir(directory):

        for file in dicom_files_in_dir(directory):

            output_file = os.path.join(output_directory, os.path.relpath(file, input_directory))

            if update_signal is not None:

                update_signal.emit(os.path.relpath(output_file, start=output_directory))

            anonymize(file, output_file, keep_private_tags=keep_private_tags)

    def walk_dir(directory):

        for entry in os.scandir(directory):

            if entry.is_dir():

                path = os.path.join(directory, entry.name)

                anonymize_dir(path)

                walk_dir(path)

    walk_dir(input_directory)

class RunThread(QtCore.QThread):

    """

    Runs a function in a thread and alerts the parent when done.

    Uses a pyqtSignal to alert the main thread of completion.

    """

    finished = QtCore.pyqtSignal(str, name="finished")

    def __init__(self, func, on_finish, *args, **kwargs):

        super(RunThread, self).__init__()

        self.args = args

        self.kwargs = kwargs

        self.func = func

        self.finished.connect(on_finish)

        self.start()

    def run(self):

        try:

            result = self.func(*self.args, **self.kwargs)

        except Exception as e:

            traceback.print_exc()

            result = e

        finally:

            if result is None:

                self.finished.emit("")

            else:

                self.finished.emit(str(result))

class Anonymizer(QtWidgets.QMainWindow):

    progress_update = QtCore.pyqtSignal(str, name="progress_update")

    def __init__(self, path=None):

        super(Anonymizer, self).__init__()

        self.setWindowTitle("DICOM anonymizer")

        central = QtWidgets.QWidget()

        self.setCentralWidget(central)

        main_layout = QtWidgets.QVBoxLayout()

        central.setLayout(main_layout)

        input_widget = QtWidgets.QWidget()

        main_layout.addWidget(input_widget)

        input_layout = QtWidgets.QHBoxLayout()

        input_widget.setLayout(input_layout)

        input_label = QtWidgets.QLabel("Input:")

        input_layout.addWidget(input_label)

        self.input_line = QtWidgets.QLineEdit(path)

        input_layout.addWidget(self.input_line)

        input_button = QtWidgets.QPushButton("Open")

        input_button.clicked.connect(lambda: self.open_directory(which="input"))

        input_layout.addWidget(input_button)

        output_widget = QtWidgets.QWidget()

        main_layout.addWidget(output_widget)

        output_layout = QtWidgets.QHBoxLayout()

        output_widget.setLayout(output_layout)

        output_label = QtWidgets.QLabel("Output:")

        output_layout.addWidget(output_label)

        self.output_line = QtWidgets.QLineEdit()

        output_layout.addWidget(self.output_line)

        output_button = QtWidgets.QPushButton("Open")

        output_button.clicked.connect(lambda: self.open_directory(which="output"))

        output_layout.addWidget(output_button)

        anonymize_button = QtWidgets.QPushButton("Anonymize")

        anonymize_button.clicked.connect(self.anonymize)

        main_layout.addWidget(anonymize_button)

        self.status_label = QtWidgets.QLabel("")

        self.statusBar().addPermanentWidget(self.status_label)

        self._task = None

        self.progress_update.connect(self.on_progress_update)

    def open_directory(self, *, which):

        dialog = QtWidgets.QFileDialog(self)

        dialog.setFileMode(QtWidgets.QFileDialog.DirectoryOnly)

        dialog.setViewMode(QtWidgets.QFileDialog.List)

        dialog.setOption(QtWidgets.QFileDialog.ShowDirsOnly, True)

        if dialog.exec_():

            directory = str(dialog.selectedFiles()[0])

            if which == "input":

                self.input_line.setText(directory)

            elif which == "output":

                self.output_line.setText(directory)

            else:

                raise Exception("invalid which='{}'".format(which))

            self.status_label.setText("")

    def anonymize(self):

        input_directory = self.input_line.text()

        if not os.path.isdir(input_directory):

            self.status_label.setText("Input directory is invalid.")

            return

        output_directory = self.output_line.text()

        # create if it does not exist

        if not os.path.exists(output_directory):

            os.makedirs(output_directory)

        # check that it is a directory

        if not os.path.isdir(output_directory):

            self.status_label.setText("Output directory is invalid.")

            return

        # check that the directory is empty

        if os.listdir(output_directory):

            self.status_label.setText("Output directory is not empty.")

            return

        if os.path.realpath(output_directory) == os.path.realpath(input_directory):

            self.status_label.setText("The output directory cannot be the same as the input directory.")

            return

        if os.path.dirname(os.path.realpath(output_directory)).startswith(os.path.realpath(input_directory)):

            self.status_label.setText("The output directory cannot be a subdirectory of the input directory.")

            return

        # run asynchronous task

        self.task = RunThread(anonymize_tree, self.on_task_finished, input_directory, output_directory, update_signal=self.progress_update)

    def on_task_finished(self, result):

        self.status_label.setAlignment(QtCore.Qt.AlignLeft)

        if not result:

            self.status_label.setText("Anonymized!")

        else:

            self.status_label.setText("Error: {}".format(result))

        self.task = None

    def on_progress_update(self, item):

        self.status_label.setAlignment(QtCore.Qt.AlignRight)

        self.status_label.setText(item)

def main():

    if len(sys.argv) < 2:

        path = "."

    else:

        path = sys.argv[1]

        if not os.path.isdir(path):

            path = "."

    path = os.path.abspath(path)

    app = QtWidgets.QApplication(sys.argv)

    QtCore.QCoreApplication.setApplicationName("DICOM anonymizer")

    QtCore.QCoreApplication.setOrganizationName("MMG FJFI")

    anonymizer = Anonymizer(path)

    anonymizer.show()

    sys.exit(app.exec_())

if __name__ == "__main__":

    main()

"""Sorts DICOMs from all four series of a 4D flow sequence

according to series number, trigger time and slice position.

Then copies the files to four new folders and renames the files so

that the slice position and trigger time are part of the file name."""

#TODO: all paths rewrite using os.sep.join

import pydicom

#import matplotlib.pyplot as plt

import numpy as numpy

import os

import shutil

import sys

import glob

#TODO: Possibly simplify times_set sorting code using bisect.insort

if __name__ == '__main__':

    import declarations as dcl

    base = dcl.base # or change to e.g.: base = ".."  +  os.sep + "new_base"

    dicom_input_folder = dcl.dicom_input_folder # or change

    dicom_output_folder = dcl.dicom_output_folder # or change, etc.

    series_tuple = dcl.series_tuple

else:

    import declarations as dcl # never change this block

    base = dcl.base

    dicom_input_folder = dcl.dicom_input_folder

    dicom_output_folder = dcl.dicom_output_folder

    series_tuple = dcl.series_tuple

    #from declarations import * # shorter alternative

#########################################################

print("Sorting dicoms to series lists... ")

series_file_lists = [[],[],[],[]]

listcount = 0

skipcount = 0

times_set = set()

slices_set = set()

for filename in glob.glob(dicom_input_folder + os.sep  +  "**"  + \

                          os.sep +  "*.dcm", recursive=True):

    f = pydicom.dcmread(filename)

    if (f.SeriesNumber in series_tuple) and hasattr(f,"SliceLocation") \

            and hasattr(f,"TriggerTime"):

        #sorts dicom filenames to four lists according to series num

        series_file_lists[series_tuple.index(f.SeriesNumber)].\

                            append(filename)

        #create set of trigger times and slice positions

        times_set.add(float(f.TriggerTime))

        slices_set.add(float(f.SliceLocation))

        print(filename)

        listcount = listcount + 1

    else:

        skipcount = skipcount + 1

#create dictionary of trigger time order number and value

times_list = list(times_set)

times_list.sort()

times_dict = { item : i for i, item in enumerate(times_list) }

slices_list = list(slices_set)

slices_list.sort()

slices_dict = { item : i for i, item in enumerate(slices_list) }

n_times = len(times_dict)   # number of time frames per each slice

n_slices = len(slices_dict) # number of slices

print("Number of files sorted to series: " + str(listcount))

print("Number of files skipped: " + str(skipcount))

print("Number of slices: " + str(n_slices))

print("Number of time frames: " + str(n_times))

##################################################################

# Copy and rename DICOM files

for i, num_of_a_series in enumerate(series_tuple):

    print("Sorting with respect to slice location and time... ")

    series_file_lists[i] = sorted(series_file_lists[i], \

                    key = lambda s: pydicom.dcmread(s).InstanceNumber)

    series_file_lists[i] = sorted(series_file_lists[i], \

                    key = lambda s: pydicom.dcmread(s).TriggerTime)

    # create output folders if they don't exists yet:

    dest_dir = dicom_output_folder + "series_" + str(num_of_a_series) + os.sep

    if not os.path.isdir(dest_dir):

        if not os.path.isdir(dicom_output_folder):

            os.mkdir(dicom_output_folder)

        os.mkdir(dest_dir)

    print("Copying and renaming files...")

    for filename in series_file_lists[i]:

        f = pydicom.dcmread(filename)

        name_core = filename[filename.rindex(os.sep)+1:] #original name

        print("Original file name: " + name_core)

        shutil.copy(filename, dest_dir) #copy the file to destin. dir

        dst_file = dest_dir + name_core

        new_dst_file_name = os.path.join(dest_dir, \

            "time_"+str(times_dict[(f.TriggerTime)]).zfill(2) + \

            "_slice_"+str(slices_dict[f.SliceLocation] + 1).zfill(2) + \

            ".dcm")

        print("New file name: " + new_dst_file_name)

        os.rename(dst_file, new_dst_file_name)

    print("Files copied and renamed.")

"""Loads dicom files sorted with respect to time and slice

position by a preceding script. Than converts them to several

3D .vti files, one for each time frame. This is done for

one subtype of 4D flow sequence (e.g. THR, inplane1, inplane2,

magnitude). Velicity images are scaled to cm/s, magnitude images

ale left unscaled. XYZ-scanner-space to image-space

transformation matrix and offset are stored to .txt files.

"""

import pydicom

import numpy as np

import glob

import numpy

from scipy.ndimage import zoom

import myVTKPythonLibrary as myvtk

import os

from pom_funkce_VTK import numpy2VTK, writeImage

#TODO: shouldn't bounds be num_of_slices*SliceThickness, etc?

#TODO: Should img3d.npy be saved?

    #import declarations as dcl # never change this block

    #base = dcl.base

    #dicom_input_folder = dcl.dicom_input_folder

    #dicom_output_folder = dcl.dicom_output_folder

    #series_tuple = dcl.series_tuple

    #seq_subtype_index = dcl.seq_subtype_index

    #seq_subtype_tuple = dcl.seq_subtype_tuple

    #VENC_tuple = dcl.VENC_tuple

    #vti_folder_prefix = dcl.vti_folder_prefix

    #from declarations import * # shorter alternative

base = "/home/katka/mnt/local/projects/IKEM/AZV-2019/data/G318_patient33_bodyPart1_visit0/G318_patient33_bodyPart1_visit0/dicom_output/series_32/"

out_folder = "/home/katka/mnt/local/projects/IKEM/AZV-2019/data/G318_patient33_bodyPart1_visit0/G318_patient33_bodyPart1_visit0/vti_32/"

arg = base + "*.dcm"

if not os.path.isdir(out_folder):

    os.mkdir(out_folder)

files = []

print('glob: {}'.format(arg))

for fname in glob.glob(arg, recursive=False):

    print("loading: {}".format(fname))

    files.append(pydicom.dcmread(fname))

    dcm = pydicom.dcmread(fname)

print("file count: {}".format(len(files)))

# skip files with no SliceLocation (eg scout views)

slices = []

skipcount = 0

n_slices = 0

for f in files:

    if hasattr(f, 'SliceLocation'):

        slices.append(f)

        n_slices = n_slices +1

    else:

        skipcount = skipcount + 1

print("n slices = ", n_slices)

print("skipped, no SliceLocation: {}".format(skipcount))

# ensure they are in the correct order

slices = sorted(slices, key=lambda s: s.SliceLocation)

# pixel aspects, assuming all slices are the same

ps = slices[0].PixelSpacing

ss = slices[0].SliceThickness

print("PixelSpacing = ", ps)

print("SliceThickness = ", ss)

# create 3D array

img_shape = list(slices[0].pixel_array.shape)

img_shape.append(len(slices))

img3d = np.zeros(img_shape)

# fill 3D array with the images from the files

for i, s in enumerate(slices):

    img2d = s.pixel_array

    img3d[:, :, n_slices - i -1] = img2d

    #img3d[:, :, i] = img2d

########################################################

vtk_image = numpy2VTK(img3d)

vtk_image.SetOrigin( (0.0, 0.0, 0.0))

vtk_image.SetSpacing((ps[0], ps[1], ss))

#TODO: shouldn't bounds be num_of_slices*SliceThickness, etc?

print("Bounds: ", vtk_image.GetBounds())

myvtk.writeImage(vtk_image, out_folder + str("img3d.vti"))

#################################################

#offest must be from last slice, in accordance with img3D filling

offset = slices[-1].ImagePositionPatient # (0020, 0032)

vectors = slices[-1].ImageOrientationPatient # (0020, 0037)

dim = slices[-1].PixelSpacing

z_dim = slices[-1].SliceThickness

vect_row = vectors[0:3]

vect_col = vectors[3:7]

vect_z = numpy.cross(vect_row, vect_col)

vect_z = numpy.cross(vect_row, vect_col)

vect_z = vect_z/ numpy.linalg.norm(vect_z)

print("offset = ", offset)

print("vect row = ", vect_row)

print("vect col = ", vect_col)

print("vect z = ", vect_z)

matrix = numpy.zeros((3,3))

matrix[1,:] = vect_row

matrix[0,:] = vect_col

matrix[2,:] = vect_z

numpy.savetxt(out_folder + "offset.txt",  offset)

numpy.savetxt(out_folder + "matrix_poz.txt",  matrix)

matrix[2,:] = -vect_z

numpy.savetxt(out_folder + "matrix_neg.txt",  matrix)

#TODO: Should img3d.npy be saved?

#numpy.save(base + out_folder + "img3d", img3d)

spacing_all = numpy.zeros(3)

spacing_all[0:2] = ps

spacing_all[2]=ss

numpy.savetxt(out_folder + "spacing.txt",  spacing_all)