Skip to content

Commit 5a60ea1

Browse files
kylekmakylekma
committed
add CT2X information
1 parent c376704 commit 5a60ea1

File tree

5 files changed

+372
-2
lines changed

5 files changed

+372
-2
lines changed

‎CT2XRAY/ReadMe.md

+30-1
Original file line numberDiff line numberDiff line change
@@ -1 +1,30 @@
1-
# Will be released soon.
1+
# Introduction
2+
This is the method used to genereate synthetic Xray images from CT volumes, as introduced in the CVPR 2019 paper: [X2CT-GAN: Reconstructing CT from Biplanar X-Rays with Generative Adversarial Networks.](https://arxiv.org/abs/1905.06902)
3+
4+
As many researchers are interested in the method, we finally took sometime to organzie the source code and release them. Sorry for the delay and let us know if you have any questions.
5+
6+
It includes three files:
7+
- ctpro.py CT pre-processing functions, mainly normalization and resampling methods.
8+
- xraypro.py DRR software script to generate synthetic Xray images from normalized CT volumes.
9+
- ctxray_utils.py Utility functions.
10+
11+
# Dependence
12+
13+
glob, scipy, numpy, cv2, pfm, matplotlib, pydicom, SiimpleITK
14+
15+
16+
# Usage Guidance
17+
18+
1. Install the DRR software from [here](https://www.plastimatch.org/) on a Windows computer, and we use 1.7.3.12-win64 version in our original work.
19+
2. In ctpro.py, first set data path;
20+
* root_path ： original CT volume root
21+
* mask_root_path ： CT mask without table
22+
* save_root_path ： save dir
23+
* then execute ctpro.py to generate the normalized CT volumes.
24+
3. In xraypro.py, again set the data path;
25+
* root_path ： the normalized data path, same as the save_root_path in ctpro.py file
26+
* save_root_path ：the Xray output path
27+
* plasti_path ： DDR software executable file location
28+
* then execute xraypro.py to generate the synthetic Xrays.
29+
4. Simple, right? If you feel our work is helpful to your project, please refer our CVPR work in your literature.
30+
5. And yes, it still not done. :) To have realistic looking Xrays, we did another CycleGAN trick for the style transfer. You can directly use the official release or waiting for our version, coming soon.

‎CT2XRAY/ctpro.py

+93
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
# ------------------------------------------------------------------------------
2+
# Copyright (c) Tencent
3+
# Licensed under the GPLv3 License.
4+
# Created by Kai Ma ([email protected])
5+
# ------------------------------------------------------------------------------
6+
import os
7+
import time
8+
import glob
9+
import scipy
10+
import numpy as np
11+
from ctxray_utils import load_scan_mhda, save_scan_mhda, bedstriping
12+
13+
# resample to stantard spacing (keep physical scale stable, change pixel numbers)
14+
def resample(image, spacing, new_spacing=[1,1,1]):
15+
# .mhd image order : z, y, x
16+
if not isinstance(spacing, np.ndarray):
17+
spacing = np.array(spacing)
18+
if not isinstance(new_spacing, np.ndarray):
19+
new_spacing = np.array(new_spacing)
20+
spacing = spacing[::-1]
21+
new_spacing = new_spacing[::-1]
22+
23+
spacing = np.array(list(spacing))
24+
resize_factor = spacing / new_spacing
25+
new_real_shape = image.shape * resize_factor
26+
new_shape = np.round(new_real_shape)
27+
real_resize_factor = new_shape / image.shape
28+
new_spacing = spacing / real_resize_factor
29+
image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')
30+
return image, new_spacing
31+
32+
# crop to stantard shape (scale x scale x scale) (keep pixel scale consistency)
33+
def crop_to_standard(scan, scale):
34+
z, y, x = scan.shape
35+
if z >= scale:
36+
ret_scan = scan[z-scale:z, :, :]
37+
else:
38+
temp1 = np.zeros(((scale-z)//2, y, x))
39+
temp2 = np.zeros(((scale-z)-(scale-z)//2, y, x))
40+
ret_scan = np.concatenate((temp1, scan, temp2), axis=0)
41+
z, y, x = ret_scan.shape
42+
if y >= scale:
43+
ret_scan = ret_scan[:, (y-scale)//2:(y+scale)//2, :]
44+
else:
45+
temp1 = np.zeros((z, (scale-y)//2, x))
46+
temp2 = np.zeros((z, (scale-y)-(scale-y)//2, x))
47+
ret_scan = np.concatenate((temp1, ret_scan, temp2), axis=1)
48+
z, y, x = ret_scan.shape
49+
if x >= scale:
50+
ret_scan = ret_scan[:, :, (x-scale)//2:(x+scale)//2]
51+
else:
52+
temp1 = np.zeros((z, y, (scale-x)//2))
53+
temp2 = np.zeros((z, y, (scale-x)-(scale-x)//2))
54+
ret_scan = np.concatenate((temp1, ret_scan, temp2), axis=2)
55+
return ret_scan
56+
57+
58+
if __name__ == '__main__':
59+
root_path = 'F:/Data/LIDC-IDRI_Volumes/'
60+
mask_root_path = 'F:/Data/LIDC-IDRI_BodyMask/'
61+
save_root_path = 'D:/LIDC_TEST'
62+
files_list = glob.glob(root_path + '*.mhd')
63+
files_list = sorted(files_list)
64+
65+
start = time.time()
66+
for fileIndex, filePath in enumerate(files_list):
67+
t0 = time.time()
68+
_, file = os.path.split(filePath)
69+
fileName = os.path.splitext(file)[0]
70+
print('Begin {}/{}: {}'.format(fileIndex+1, len(files_list), fileName))
71+
maskFile = fileName + '_outMultiLabelMask' + '.mhd'
72+
maskPath = os.path.join(mask_root_path, maskFile)
73+
saveDir = os.path.join(save_root_path, fileName)
74+
if not os.path.exists(saveDir):
75+
os.makedirs(saveDir)
76+
savePath = os.path.join(saveDir, '{}.mha'.format(fileName))
77+
ct_itk, ct_scan, ori_origin, ori_size, ori_spacing= load_scan_mhda(filePath)
78+
print("Old : ", ori_size)
79+
_, ct_mask, _, _, _ = load_scan_mhda(maskPath)
80+
81+
bedstrip_scan = bedstriping(ct_scan, ct_mask)
82+
new_scan, new_spacing = resample(bedstrip_scan, ori_spacing, [1, 1, 1])
83+
84+
print("Std : ", new_scan.shape[::-1])
85+
std_scan = crop_to_standard(new_scan, scale=320)
86+
87+
save_scan_mhda(std_scan, (0, 0, 0), new_spacing, savePath)
88+
_, _, _, new_size, new_spacing = load_scan_mhda(savePath)
89+
print("New : ", new_size)
90+
t1 = time.time()
91+
print('End! Case time: {}'.format(t1-t0))
92+
end = time.time()
93+
print('Finally! Total time: {}'.format(end-start))

‎CT2XRAY/ctxray_utils.py

+146
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
# ------------------------------------------------------------------------------
2+
# Copyright (c) Tencent
3+
# Licensed under the GPLv3 License.
4+
# Created by Kai Ma ([email protected])
5+
# ------------------------------------------------------------------------------
6+
import os
7+
import cv2
8+
import time
9+
import pydicom
10+
import numpy as np
11+
import SimpleITK as sitk
12+
13+
# Load the scans in given folder path
14+
def load_scan_dcm(path):
15+
slices = [pydicom.dcmread(os.path.join(path, s)) for s in os.listdir(path) if 'dcm' in s]
16+
slices.sort(key=lambda x: int(x.ImagePositionPatient[2]))
17+
18+
try:
19+
slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
20+
except:
21+
slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
22+
23+
for s in slices:
24+
if hasattr(s, 'SliceThickness'):
25+
pass
26+
else:
27+
s.SliceThickness = slice_thickness
28+
return slices
29+
30+
# input could be .mhd/.mha format
31+
def load_scan_mhda(path):
32+
img_itk = sitk.ReadImage(path)
33+
img = sitk.GetArrayFromImage(img_itk)
34+
return img_itk, img, img_itk.GetOrigin(), img_itk.GetSize(), img_itk.GetSpacing()
35+
36+
# output coule be .mhd/.mha format
37+
def save_scan_mhda(volume, origin, spacing, path):
38+
itkimage = sitk.GetImageFromArray(volume, isVector=False)
39+
itkimage.SetSpacing(spacing)
40+
itkimage.SetOrigin(origin)
41+
sitk.WriteImage(itkimage, path, True)
42+
43+
# do betstriping
44+
def bedstriping(ct_scan, ct_mask):
45+
bedstrip_scan = ct_scan * ct_mask
46+
return bedstrip_scan
47+
48+
# build CycleGAN trainSet A
49+
def buildTrainA(raw_root_path, save_root_path):
50+
files_list = os.listdir(raw_root_path)
51+
if not os.path.exists(save_root_path):
52+
os.makedirs(save_root_path)
53+
start = time.time()
54+
for fileIndex, fileName in enumerate(files_list):
55+
t0 = time.time()
56+
print('Begin {}/{}: {}'.format(fileIndex+1, len(files_list), fileName))
57+
imageDir = os.path.join(raw_root_path, fileName)
58+
imagePath = os.path.join(imageDir, 'xray1.png')
59+
image = cv2.imread(imagePath, 0)
60+
savePath = os.path.join(save_root_path, '{}.png'.format(fileName))
61+
cv2.imwrite(savePath, image)
62+
t1 = time.time()
63+
print('End! Case time: {}'.format(t1-t0))
64+
end = time.time()
65+
print('Finally! Total time of build TrainA: {}'.format(end-start))
66+
67+
# build CycleGAN trainSet B
68+
def buildTrainB(raw_root_path, save_root_path):
69+
BSize = 1000
70+
files_list = os.listdir(raw_root_path)
71+
if not os.path.exists(save_root_path):
72+
os.makedirs(save_root_path)
73+
tb = np.arange(0, len(files_list))
74+
np.random.shuffle(tb)
75+
start = time.time()
76+
for i in range(0, BSize):
77+
t0 = time.time()
78+
imageName = files_list[tb[i]]
79+
print('Begin {}/{}: {}'.format(i+1, BSize, imageName))
80+
imagePath = os.path.join(raw_root_path, imageName)
81+
image = cv2.imread(imagePath, 0)
82+
image = cv2.resize(image, (320, 320))
83+
savePath = os.path.join(save_root_path, imageName)
84+
cv2.imwrite(savePath, image)
85+
t1 = time.time()
86+
print('End! Case time: {}'.format(t1-t0))
87+
end = time.time()
88+
print('Finlly! Total time of build TrainB: {}'.format(end-start))
89+
90+
# resize all images in a directory
91+
def resize_to_standard(path):
92+
images = os.listdir(path)
93+
for i in images:
94+
imageName = os.path.join(path, i)
95+
image = cv2.imread(imageName, 0)
96+
image = cv2.resize(image, (320, 320))
97+
cv2.imwrite(imageName, image)
98+
print('All images have been resized!')
99+
100+
def convert2hu(path):
101+
img, scan, origin, size, spacing = load_scan_mhda(path)
102+
#print(scan.dtype)
103+
scan = scan.astype(np.int16)
104+
scanhu = scan - 1024
105+
root, file = os.path.split(path)
106+
fileName = file.replace('.mha', 'hu')
107+
#print(root, file)
108+
savePath = os.path.join(root, fileName+'.mha')
109+
save_scan_mhda(scanhu, origin, spacing, savePath)
110+
111+
def psnr(img_1, img_2, PIXEL_MAX = 255.0):
112+
mse = np.mean((img_1 - img_2) ** 2)
113+
if mse == 0:
114+
return 100
115+
return 20 * np.log10(PIXEL_MAX / np.sqrt(mse))
116+
117+
118+
if __name__ == '__main__':
119+
raw_root_pathA = 'D:/LIDC_TEST'
120+
save_root_pathA = 'D:/xrayv2r/trainA'
121+
raw_root_pathB = 'F:/Data/NIHCXR/images_002/images'
122+
save_root_pathB = 'D:/xrayv2r/trainB'
123+
#buildTrainA(raw_root_pathA, save_root_pathA)
124+
#buildTrainB(raw_root_pathB, save_root_pathB)
125+
126+
'''
127+
raw_root_pathA = 'D:/xrayv2r/trainA500'
128+
ASize = 260
129+
files_list = os.listdir(raw_root_pathA)
130+
if not os.path.exists(save_root_pathA):
131+
os.makedirs(save_root_pathA)
132+
start = time.time()
133+
for i in range(0, ASize):
134+
t0 = time.time()
135+
imageName = files_list[i]
136+
print('Begin {}/{}: {}'.format(i+1, ASize, imageName))
137+
imagePath = os.path.join(raw_root_pathA, imageName)
138+
image = cv2.imread(imagePath, 0)
139+
#image = cv2.resize(image, (320, 320))
140+
savePath = os.path.join(save_root_pathA, imageName)
141+
cv2.imwrite(savePath, image)
142+
t1 = time.time()
143+
print('End! Case time: {}'.format(t1-t0))
144+
end = time.time()
145+
print('Finlly! Total time of build TrainB: {}'.format(end-start))
146+
'''

‎CT2XRAY/xraypro.py

+102
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
# ------------------------------------------------------------------------------
2+
# Copyright (c) Tencent
3+
# Licensed under the GPLv3 License.
4+
# Created by Kai Ma ([email protected])
5+
# ------------------------------------------------------------------------------
6+
import os
7+
import cv2
8+
import pfm
9+
import time
10+
import numpy as np
11+
import matplotlib.pyplot as plt
12+
from subprocess import check_output as qx
13+
from ctxray_utils import load_scan_mhda, save_scan_mhda
14+
15+
16+
# compute xray source center in world coordinate
17+
def get_center(origin, size, spacing):
18+
origin = np.array(origin)
19+
size = np.array(size)
20+
spacing = np.array(spacing)
21+
center = origin + (size - 1) / 2 * spacing
22+
return center
23+
24+
# convert a ndarray to string
25+
def array2string(ndarray):
26+
ret = ""
27+
for i in ndarray:
28+
ret = ret + str(i) + " "
29+
return ret[:-2]
30+
31+
# save a .pfm file as a .png file
32+
def savepng(filename, direction):
33+
raw_data , scale = pfm.read(filename)
34+
max_value = raw_data.max()
35+
im = (raw_data / max_value * 255).astype(np.uint8)
36+
# PA view should do additional left-right flip
37+
if direction == 1:
38+
im = np.fliplr(im)
39+
40+
savedir, _ = os.path.split(filename)
41+
outfile = os.path.join(savedir, "xray{}.png".format(direction))
42+
# plt.imshow(im, cmap=plt.cm.gray)
43+
plt.imsave(outfile, im, cmap=plt.cm.gray)
44+
# plt.imsave saves an image with 32bit per pixel, but we only need one channel
45+
image = cv2.imread(outfile)
46+
gray = cv2.split(image)[0]
47+
cv2.imwrite(outfile, gray)
48+
49+
50+
if __name__ == '__main__':
51+
root_path = 'D:/LIDC_TEST/'
52+
save_root_path = 'D:/LIDC_TEST'
53+
plasti_path = 'D:/Program Files/Plastimatch/bin'
54+
55+
files_list = os.listdir(root_path)
56+
start = time.time()
57+
for fileIndex, fileName in enumerate(files_list):
58+
t0 = time.time()
59+
print('Begin {}/{}: {}'.format(fileIndex+1, len(files_list), fileName))
60+
saveDir = os.path.join(root_path, fileName)
61+
if not os.path.exists(saveDir):
62+
os.makedirs(saveDir)
63+
# savePath is the .mha file store location
64+
savePath = os.path.join(saveDir, '{}.mha'.format(fileName))
65+
ct_itk, ct_scan, ori_origin, ori_size, ori_spacing = load_scan_mhda(savePath)
66+
# compute isocenter
67+
center = get_center(ori_origin, ori_size, ori_spacing)
68+
# map the .mha file value to (-1000, 1000)
69+
cmd_str = plasti_path + '/plastimatch adjust --input {} --output {}\
70+
--pw-linear "0, -1000"'.format(savePath, saveDir+'/out.mha')
71+
output = qx(cmd_str)
72+
# get virtual xray file
73+
directions = [1, 2]
74+
for i in directions:
75+
if i == 1:
76+
nrm = "0 1 0"
77+
else:
78+
nrm = "1 0 0"
79+
'''
80+
plastimatch usage
81+
-t : save format
82+
-g : sid sad [DistanceSourceToPatient]:541
83+
[DistanceSourceToDetector]:949.075012
84+
-r : output image resolution
85+
-o : isocenter position
86+
-z : physical size of imager
87+
-I : input file in .mha format
88+
-O : output prefix
89+
'''
90+
cmd_str = plasti_path + '/drr -t pfm -nrm "{}" -g "541 949" \
91+
-r "320 320" -o "{}" -z "500 500" -I {} -O {}'.format\
92+
(nrm, array2string(center), saveDir+'/out.mha', saveDir+'/{}'.format(i))
93+
output = qx(cmd_str)
94+
# plastimatch would add a 0000 suffix
95+
pfmFile = saveDir+'/{}'.format(i) + '0000.pfm'
96+
savepng(pfmFile, i)
97+
# remove the temp .mha file couse it is so large
98+
os.remove(saveDir+'/out.mha')
99+
t1 = time.time()
100+
print('End! Case time: {}'.format(t1-t0))
101+
end = time.time()
102+
print('Finally! Total time: {}'.format(end-start))

‎ReadMe.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,7 @@ Qualitative results from our original paper.
135135
- [x] Testing code example and script of the algorithm
136136
- [x] Visualization code example and script of the algorithm
137137
- [x] Pre-processed LIDC data upload to cloud (70G, training and test data used in the CVPR work)
138-
- [ ] Source code to generate synthesized X-Rays from CT volumes
138+
- [x] Source code to generate synthesized X-Rays from CT volumes
139139
- [ ] Source code to generate realistic X-Rays by using CycleGAN
140140

141141
## Acknowledgement

0 commit comments

Comments
 (0)