microCT_auto_processing/select_proximal_TF_junction_slices_registrated.py at main · Tiankuo528/microCT_auto_processing · GitHub

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
import numpy as np
from skimage import io
import cv2
import os
import matplotlib.pyplot as plt

def midshaft(img2d):

    image = img2d

#smoothing images
    image=np.array(image)
    ret, thresh = cv2.threshold(image, 0, 255, cv2.THRESH_BINARY)
    smoothedimage = cv2.medianBlur(thresh,9)  # 9  21    45
    #The size of the kernel affects which holes get closed
    kernel_size = 45 #5    13 20 49
    kernel = np.ones((kernel_size, kernel_size), np.uint8)
    smoothedimage = cv2.morphologyEx(smoothedimage, cv2.MORPH_CLOSE, kernel)


    # Define the region (mask) at the bottom of the image where particles are located
    height, width = smoothedimage.shape
    mask_height = int(0.2*height)  # Define how tall the mask is from the bottom
    mask = np.zeros_like(thresh)
    mask[height-mask_height:height, :] = 1
    # Remove the small particles by applying the mask
    smoothedimage[mask == 1] = 0

    # plt.imshow(smoothedimage,cmap="gray")
    # plt.show()

 #identify the mid-shaft position

    contours, hierarchy = cv2.findContours(smoothedimage,cv2.RETR_EXTERNAL ,cv2.CHAIN_APPROX_NONE)    #cv2.RETR_EXTERNAL
    count= len(contours)


    if count > 1 :

        # print(i)
        # smoothedimage = np.stack((smoothedimage,)*3, axis=-1)
        # cv2.drawContours(smoothedimage, contours, -1, (0, 255, 0), 2)
        # cv2.imshow('Contours', smoothedimage)
        # cv2.waitKey(0)
        # cv2.destroyAllWindows()

        return False

    elif count == 1 :
        # print(i)
        # smoothedimage = np.stack((smoothedimage,)*3, axis=-1)
        # cv2.drawContours(smoothedimage, contours, -1, (0, 255, 0), 2)
        # cv2.imshow('Contours', smoothedimage)
        # cv2.waitKey(0)
        # cv2.destroyAllWindows()


        return True

config = {
    "foldername" : "selected below proximal t-f junction 0-300" ,       #40-110
    "first_slice_selected"  : "first slice selected",
    "last_slice_selected"  : "last slice selected",
    "first_slice_selected_above_t-f_junction" : 0  #x: first selected slice index above mid-shaft    :50, 450
}
input= r"D:\1c.  Aged Rad+Yoda1+WBV+ Tumor 101223\registration\Tibia w2w5composite"    #folder of input composite w0w_composite
foldername = config["foldername"]
first_slice_selected = config["first_slice_selected"]     # we select 70 slices each time,
last_slice_selected = config["last_slice_selected"]
first_slice_selected_above_mid_shaft = config["first_slice_selected_above_t-f_junction"] #the # of slices starting select above tf- junction

if not  os.path.exists(os.path.join(input,foldername)):
    os.mkdir(os.path.join(input,foldername))
if not  os.path.exists(os.path.join(input,foldername,"slice of t-f junction")):
    os.mkdir(os.path.join(input,foldername,"slice of t-f junction"))
if not  os.path.exists(os.path.join(input,foldername,first_slice_selected)):
    os.mkdir(os.path.join(input,foldername,first_slice_selected))
if not  os.path.exists(os.path.join(input,foldername,last_slice_selected)):
    os.mkdir(os.path.join(input,foldername,last_slice_selected))

for folder in os.listdir(input):
    if not  os.path.exists(os.path.join(input,foldername,folder)):
        os.mkdir(os.path.join(input,foldername,folder))
    folderpath=os.path.join(input,folder)

    try:
        img=[]
        for filename in os.listdir(folderpath):
            if filename.endswith('.tif'):
                rawimage = cv2.imread(os.path.join(folderpath,filename),0)

        #rawimage = cv2.cvtColor(rawimage, cv2.COLOR_BGR2GRAY)
                img.append(rawimage)

        for i in range(len(img)): # shape[0]: read image in grayscale
            #shape = [z,y,x] 3d
            if i > 1250 : #1250: # 1150:
                img2d = img[i]

                if midshaft(img2d):

                    # select slice at mid-shaft and last image selected
                    io.imsave(os.path.join(input,foldername,"slice of t-f junction",os.listdir(folderpath)[i]),img[i])    # output images of t-f junction
                    io.imsave(os.path.join(input,foldername,first_slice_selected,os.listdir(folderpath)[i+first_slice_selected_above_mid_shaft]),img[i + first_slice_selected_above_mid_shaft])
                    io.imsave(os.path.join(input,foldername,last_slice_selected,os.listdir(folderpath)[i+first_slice_selected_above_mid_shaft-300]),img[i + first_slice_selected_above_mid_shaft-300])
                    #io.imsave(os.path.join(input,foldername,slice_above_mid_shaft_folder,os.listdir(folderpath)[i+first_slice_selected_above_mid_shaft+50]),img[i + first_slice_selected_above_mid_shaft +50])    # output last one of selected 50 slices of mid-shaft
                    # select 50 slices from xxx (ex.450) slices above mid-shaft
                    j=i+first_slice_selected_above_mid_shaft  # j=i+x  x: first selected slice index ############################### start select
                    for k in range(300):  #x : # of slices selected
                        per_selected_img= img[j-k]
                        io.imsave(os.path.join(input,foldername,folder,os.listdir(folderpath)[j-k]),per_selected_img)    # output images
                    break
    except:
        print("Error happens:cannot select enough slices",folder)


#folder = r"C:\Users\SplLabMEEG\Desktop\03.07.2022 Aged non tumor Yoda1 Veh loading\selected mice\5. mkcomposite\Tibia w0w4composite\101 week 4 left tibia w0w4composite"   #shape = [z,y,x]  3d
#img=[]
# for filename in os.listdir(folder):
#     rawimage = cv2.imread(os.path.join(folder,filename),0)
#     #rawimage = cv2.cvtColor(rawimage, cv2.COLOR_BGR2GRAY)
#     img.append(rawimage)


#save images