Is there a problem with the cone beam CT data import results that are all nan values?

[wwwaaa32498]

I use a Metrotom 1500 industrial CT scanner from ZEISS Measuring machine. The STL model generated from the reconstructed slice map using ZEISS software is shown in the figure below. It is 6 single balls.

The following figure shows the parameters and projected data. The projected data is 2050 images. That is, a projection every 360/2050 degrees.

I only modified the parameters and the import of the image, why the reconstructed 3D is all "nan". Below is a slice of my code and the reconstructed 102410241024 3D volume data

import matplotlib.pyplot as plt import numpy as np from scipy.spatial.transform import Rotation as R import astra import os import imageio.v2 as imageio import time import utils import cv2 as cv

x_min, x_max = -204.8,204.8 # [mm] reconstruction volume range y_min, y_max = -204.8, 204.8 # [mm] reconstruction volume range z_min, z_max = -204.8, 204.8 # [mm] reconstruction volume range x_vol_sz, y_vol_sz, z_vol_sz = 1024, 1024, 1024 # volex number for x, y, z. Please make sure x_vol_sz=y_vol_sz

n_pro = 2050

ard_rot = True

projs_cols = 1024 projs_rows = 1024 det_spacing_x = 0.4 # [mm] det_spacing_y = 0.4 # [mm]

dist_orgin_detector = 1064.629 dist_source_origin = 314.341 # [mm] dist_origin_detector = dist_orgin_detector - dist_source_origin # [mm] magnification = dist_orgin_detector / dist_source_origin # [mm]

src_x_det_crd = -1.25 # [mm] src_y_det_crd = 0.88 # [mm] src_z_det_crd = dist_orgin_detector # [mm]

rot_x_det_crd = -0.586413 # [mm] rot_y_det_crd = -0.74428 # [mm] rot_z_det_crd = dist_origin_detector # [mm]

t = time.time() print('load data', flush=True) angles = np.linspace(0, 2 * np.pi, num=n_pro, endpoint=False) projs = np.zeros((projs_rows, n_pro, projs_cols), dtype=np.float32) # (n_pro, projs_rows, projs_cols)

imageSize = (projs_cols, projs_rows) if ard_rot: imageSize = list(reversed(imageSize))

for i in range(n_pro):

id = 'D:/shiyanshuju/JLDX_SLQ_01 2024-1-8 13-18-38/Corrected' + str(i) + '.uint16'
fid = open(id, 'rb')
data = np.fromfile(fid, dtype=np.uint16, count=2048 * 2048)
# input_file = os.path.join(output_dir, f'projection_{i:03d}.jpg')
# projection = imread(input_file, as_gray=True)
# projections[i] = projection
data = data.reshape((2048, 2048))
data = np.transpose(data)
# data = misc.imresize(data, 0.5)
data = cv.resize(data, (1024, int(1024)))
data = -np.log(data / 65535)
projs[:, i, :] = data  # (projs_rows, n_pro, projs_cols)
print(i)

print(np.round_(time.time() - t, 3), 'sec elapsed')

t = time.time() print('pre-process data', flush=True)

np.log(projs, out=projs) np.negative(projs, out=projs) projs = np.ascontiguousarray(projs) print(np.round_(time.time() - t, 3), 'sec elapsed')

t = time.time() print('compute reconstruction', flush=True)

vol_sz = [x_vol_sz, y_vol_sz, z_vol_sz] # CStalk vol_rec = np.zeros(vol_sz, dtype=np.float32)

vol_geom = astra.create_vol_geom(vol_sz) vol_geom['option']['WindowMinX'] = x_min vol_geom['option']['WindowMaxX'] = x_max vol_geom['option']['WindowMinY'] = y_min vol_geom['option']['WindowMaxY'] = y_max vol_geom['option']['WindowMinZ'] = z_min vol_geom['option']['WindowMaxZ'] = z_max

angles = np.linspace(0, 2*np.pi, n_pro, False) def cal_vecs(src_x_det_crd, src_y_det_crd, src_z_det_crd, rot_x_det_crd, rot_y_det_crd, rot_z_det_crd, angles, det_spacing_x, det_spacing_y):

src_pos = [-rot_x_det_crd, -(src_z_det_crd - rot_z_det_crd),
           -src_y_det_crd * (src_z_det_crd - rot_z_det_crd) / src_z_det_crd]
det_pos = [-src_x_det_crd-rot_x_det_crd, rot_z_det_crd, src_y_det_crd * rot_z_det_crd / src_z_det_crd]

alpha = 0
theta = 0
gmma = 0

u = [det_spacing_x, 0, 0]
v = [0, 0, det_spacing_y]

r = R.from_rotvec([alpha, 0, 0])
u = r.apply(u)
v = r.apply(v)
r = R.from_rotvec([0, theta, 0])
u = r.apply(u)
v = r.apply(v)
r = R.from_rotvec([0, 0, gmma])
u = r.apply(u)
v = r.apply(v)

vectors = np.zeros((len(angles), 12))

for i in range(len(angles)):
    vectors[i, 0] = np.cos(angles[i]) * src_pos[0] - np.sin(angles[i]) * src_pos[1]
    vectors[i, 1] = np.sin(angles[i]) * src_pos[0] + np.cos(angles[i]) * src_pos[1]
    vectors[i, 2] = src_pos[2]

    vectors[i, 3] = np.cos(angles[i]) * det_pos[0] - np.sin(angles[i]) * det_pos[1]
    vectors[i, 4] = np.sin(angles[i]) * det_pos[0] + np.cos(angles[i]) * det_pos[1]
    vectors[i, 5] = det_pos[2]

    vectors[i, 6] = np.cos(angles[i]) * u[0]
    vectors[i, 7] = np.sin(angles[i]) * u[0]
    vectors[i, 8] = u[2]

    vectors[i, 9] = np.cos(angles[i]) * v[0]
    vectors[i, 10] = np.sin(angles[i]) * v[0]
    vectors[i, 11] = v[2]

return vectors

vecs = cal_vecs(src_x_det_crd, src_y_det_crd, src_z_det_crd, rot_x_det_crd, rot_y_det_crd, rot_z_det_crd, angles, det_spacing_x, det_spacing_y) proj_geom = astra.create_proj_geom('cone_vec', projs_rows, projs_cols, vecs)

vol_id = astra.data3d.link('-vol', vol_geom, vol_rec) proj_id = astra.data3d.link('-sino', proj_geom, projs)

cfg_fdk = astra.astra_dict('FDK_CUDA')

cfg_fdk['ProjectionDataId'] = proj_id cfg_fdk['ReconstructionDataId'] = vol_id cfg_fdk['option'] = {} cfg_fdk['option']['ShortScan'] = False alg_id = astra.algorithm.create(cfg_fdk)

astra.algorithm.run(alg_id, 1)

ccc= vol_rec[512,:,:] astra.algorithm.delete(alg_id) astra.data3d.delete(proj_id) astra.data3d.delete(vol_id)

print(np.round_(time.time() - t, 3), 'sec elapsed') mask = utils.gen_mask(vol_rec, x_min, x_max, magnification, projs_cols, det_spacing_x, src_x_det_crd, rot_x_det_crd) vol_rec = utils.remove_edge_noise(vol_rec, mask)

for i in range(1024): ccc = vol_rec[i,:,:] filename = 'C:/Users/HP/Desktop/try/18/png/proj' + str(i+1) + '.png' # 修改文件名以包含i plt.imsave(file_name, ccc, cmap='gray',format='png') t = time.time()

Brother, can you help me? Much appreciated

[mingdianliu]

First of all, could you please star this repo if it is helpful to your project?

Just want to know if you double-check the imported raw data and ensure all the raw images are correct. You can visualize them one by one to check them.

I can also share our example data if you need it.

[wwwaaa32498]

Thanks for the reply,i had star this repo already. Is this caused by my data import?What should I do with the data to import properly？ data /= 65535.0； or data = -np.log(data / 65535)； I looked at the projection in python and the ball was low, but the rest of the place was close to 65535

[mingdianliu]

It seems your data is reasonable. I can share my example data with you as a reference. Could you please send me an email? lmdvigor@gmail.com, I will send you the data link.

[wwwaaa32498]

Thank you. I was able to rebuild the package in "TIGRE-master", but in "ASTRA Toolbox" I don't know why

[mingdianliu]

You can install astra package in your environment. The environment can be built up by environment.yml file.