Commit 280a10fb authored by Radek Galabov's avatar Radek Galabov
Browse files

Resolved troubles with matrix_neg/matrix_poz and order of slices when creating VTIs.

parent e88277de

dicom_codes_without_MG/_02_dicom_to_vti_4D.py

+41 −7
Original line number Diff line number Diff line
@@ -7,6 +7,18 @@ ale left unscaled. XYZ-scanner-space to image-space 
transformation matrix and offset are stored to .txt files.

"""



def relative_direction_of_vectors(x,y):

        '''Returns True if  x  points in the same

        semispace  as  y.  Returns False  if  not.

        Returns False if perpendicular.

        '''

        cosine = np.dot(x,y) / (np.linalg.norm(x) * np.linalg.norm(y))

        if  (int(np.sign(cosine))  ==  1):

            return  True

        elif  ((int(np.sign(cosine)) ==  -1)  or

                   (int(np.sign(cosine))  ==  0)):

            return  False



import pydicom

import numpy as np

import glob
@@ -96,14 +108,32 @@ for n in range(n_times): 
    # create 3D array

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

    img_shape.append(len(slices))

    img_shape  = [img_shape[1],  img_shape[0], img_shape[2]]

    img3d = np.zeros(img_shape)



    # check  chirality  of coordinate system

    slice_1  =  slices[0]

    slice_2  =  slices[1]

    slice_1_row_vec = [float(element)  for  element

            in  (slice_1.ImageOrientationPatient)[0:3]]

    slice_1_col_vec = [float(element)  for  element

            in  (slice_1.ImageOrientationPatient)[3:6]]

    acquisition_thr_vec =  (np.array(slice_2.ImagePositionPatient)

                            - np.array(slice_1.ImagePositionPatient))

    proper_right_hand_thr_vec  =  np.cross(slice_1_row_vec,slice_1_col_vec)

                                           

    are_paralell  =  relative_direction_of_vectors(

            proper_right_hand_thr_vec, acquisition_thr_vec)

    

        

    # fill 3D array with the images from the files

    for i, s in enumerate(slices):

        img2d = s.pixel_array

        img3d[:, :, i] = img2d #n_slices -i -1 misto i bylo puvodne

        #img3d[:, :, i] = img2d

        if  are_paralell:

            img3d[:, :, i] = np.matrix.transpose(img2d) 

        elif  not  are_paralell:

            img3d[:, :, n_slices -i -1] = np.matrix.transpose(img2d) 

        





    ########################################################
@@ -199,7 +229,11 @@ for n in range(n_times): 




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

    print('are_paralell ==  ',are_paralell)

    if  (are_paralell):

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

    elif  (not  are_paralell):

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

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

    dim = slices[0].PixelSpacing

    z_dim = slices[0].SliceThickness
@@ -213,7 +247,7 @@ 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 = numpy.cross(vect_row, vect_col)

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



print("vect row = ", vect_row)
@@ -221,8 +255,8 @@ print("vect col = ", vect_col) 
print("vect z = ", vect_z)



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

matrix[1,:] = vect_row

matrix[0,:] = vect_col

matrix[0,:] = vect_row

matrix[1,:] = vect_col

matrix[2,:] = vect_z



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