Unable to reconstruct - industrial cone-beam scan

[Resinato]

Dear all, I am trying without success to reconstruct a cone beam scan made with an industrial scanner (RX solutions Easytom XL).
I am trying with github copilot: succesfully installed and running trigrepy with anaconda, binned the projections ot 4x to save time, and trying this script:

import numpy as np
import glob
from PIL import Image
import tigre
import tigre.algorithms as algs
import matplotlib.pyplot as plt
import sys
import os

# -------- USER PARAMETERS --------
folder = r"G:\Proj_bin4"
pattern = "*.tiff"
num_projs = 1200
nVoxel = np.array([256,256,256])

bin_factor = 4
flat_path = None    # e.g., r"G:flat_bin4.tiff"
dark_path = None    # e.g., r"G:dark_bin4.tiff"

# --- LOAD FLAT/DARK (OPTIONAL) ---
flat = np.array(Image.open(flat_path)).astype(np.float32) if flat_path else None
dark = np.array(Image.open(dark_path)).astype(np.float32) if dark_path else None

# --- LOAD & LOG-CORRECT PROJECTIONS ---
filelist = sorted(glob.glob(f"{folder}\\{pattern}"))[:num_projs]
print(f"Found {len(filelist)} projection files.")
if len(filelist) < num_projs:
    print("Warning: Not enough projection files found!")

logcorr_list = []
for idx, filename in enumerate(filelist):
    img = np.array(Image.open(filename)).astype(np.float32)
    if flat is not None:
        _flat = flat
        _dark = dark if dark is not None else 0.0
        corrected = (img - _dark) / (_flat - _dark + 1e-8)
    else:
        # No flat: use mean of the top rows as open beam estimate (air region)
        open_beam_est = np.max(img)  # il massimo di tutte le immagini
        corrected = img / (open_beam_est + 1e-8)
    corrected = np.clip(corrected, 1e-6, 1.0)
    log_proj = -np.log(corrected)
    logcorr_list.append(log_proj)
    if idx == 0 or (idx+1) % 50 == 0 or idx+1 == len(filelist):
        print(f"Processed {idx+1} / {len(filelist)}")

projections = np.stack(logcorr_list)
print("Projections shape (native):", projections.shape)
projections = np.transpose(projections, (0,2,1))
print("Projections shape (after transpose):", projections.shape)
print("dtype:", projections.dtype)
print("min/max:", projections.min(), projections.max())
print("Any NaN?", np.isnan(projections).any())
print("Any Inf?", np.isinf(projections).any())
sys.stdout.flush()

# Quick check: show result
plt.imshow(projections[0], cmap="gray")
plt.title("Log-corrected projection 0")
plt.show()
plt.hist(projections[0].ravel(), bins=200)
plt.title("Histogram of log-corrected projection 0")
plt.show()

angles = np.linspace(0, 2*np.pi, projections.shape[0], endpoint=False)
print("Angles shape:", angles.shape)

# --- GEOMETRY & ISOTROPIC VOXEL CALC ---
original_nDetector = np.array([2456, 1944])
original_dDetector = np.array([0.124, 0.124])
geo = tigre.geometry()
geo.nDetector = original_nDetector // bin_factor
geo.dDetector = original_dDetector * bin_factor
geo.sDetector = geo.nDetector * geo.dDetector
geo.DSO = 936.444
geo.DSD = 1131.2

# --- Isotropic voxel calculation ---
det_width = geo.nDetector[0] * geo.dDetector[0]
magnif = geo.DSD / geo.DSO
fov_object = det_width / magnif
isotropic_voxel = fov_object / nVoxel[0]
geo.nVoxel = nVoxel
geo.dVoxel = np.array([isotropic_voxel] * 3)
geo.sVoxel = geo.nVoxel * geo.dVoxel
geo.offOrigin = np.array([0,0,0])
geo.offDetector = np.array([0,0])

print("geo.nDetector:", geo.nDetector)
print("geo.dDetector:", geo.dDetector)
print("geo.DSO:", geo.DSO, "geo.DSD:", geo.DSD)
print(f"Magnification: {magnif:.3f}")
print(f"Field of view (object plane): {fov_object:.2f} mm")
print(f"Isotropic voxel size: {isotropic_voxel:.4f} mm")
print("geo.nVoxel:", geo.nVoxel)
print("geo.dVoxel:", geo.dVoxel)
sys.stdout.flush()

# --- RECONSTRUCTION ---
print("\nRunning FDK reconstruction ...")
try:
    recon = algs.fdk(projections, geo, angles)
    print("Recon shape:", recon.shape)
except Exception as e:
    print(f"Recon error: {e}")

plt.imshow(recon[recon.shape[0]//2], cmap='gray')
plt.title("Reconstruction center slice")
plt.show()

# --- SAVE RAW FOR FIJI ---
out_dir = os.path.join(folder, "recon_single")
os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, "recon_fdk_isotropic.raw")
recon_norm = recon - recon.min()
if recon_norm.max() > 0:
    recon_norm /= recon_norm.max()
else:
    print("Warning: recon is flat!")
recon_uint16 = (recon_norm * 65535).astype(np.uint16)
recon_uint16.tofile(out_path)
print(f"Saved: {out_path}")
print("FIJI import shape: (Z,Y,X) =", recon_uint16.shape)
print("--- FIJI RAW Import ---")
print(f"Width = {nVoxel[2]}")
print(f"Height = {nVoxel[1]}")
print(f"Number = {nVoxel[0]}")
print("16-bit unsigned, little-endian")
print(f"Output folder: {out_dir}")`

I have no significant results:

[AnderBiguri]

hi @Resinato . It looks like a geometry definition problem, to be fair. But its quite impossible for me to know what exactly...

[Resinato]

thanks, I am new to this field.
It was a scan made of 4 tiles, 2400 proj each.

I tried sweeping from -30 to +30 mm the X offset with no good result.

Can this xml help?

<conebeam> <reconstruct methode="Filtered Backprojection (Tomography)" iterations="0" filter="Tukey" filterValue="28" horizontalBorderFilterEnabled="false" verticalBorderFilterEnabled="true" appodizationFilter="Sinus" appodization="75" vgexport="1" rawfileexport="1" saveroi="0" rangeMin="-0.18904000000000001" rangeMax="3.1308600000000002" /> <volume_acquisition> <voxelSize X="0.102646001" Y="0.102646001" Z="0.102646001" /> <size X="2434" Y="5574" Z="2434" /> <offset X="0" Y="225.182007" Z="0" /> <rotate X="0" Y="0" Z="0" /> <testSlice X="1217" Y="2788" Z="1217" /> </volume_acquisition> <volume> <voxelSize X="0.100950077" Y="0.100950077" Z="0.100950077" /> <size X="1014" Y="5810" Z="2015" /> <offset X="38.6875458" Y="226.799561" Z="-4.8849225" /> <rotate X="1.59390676" Y="1.41304243" Z="11.5893631" /> <testSlice X="530" Y="700" Z="1003" /> </volume> <profile acquisition="Stacked tomography with n turns" images="9600" turns="4" filterCorrectionPerStack="0"> <turnParam nbImages="2400" index="0"> <avr Vieworks="1" /> </turnParam> <turnParam nbImages="2400" index="1"> <avr Vieworks="1" /> </turnParam> <turnParam nbImages="2400" index="2"> <avr Vieworks="1" /> </turnParam> <turnParam nbImages="2400" index="3"> <avr Vieworks="1" /> </turnParam> </profile> <acquisitionInfo> <xray current="882" voltage="100" type="X-Ray Source" /> <geometry sod="936.463013" sdd="1131.28003" /> </acquisitionInfo> <correction offsetXCorrectionIndex="1" helicoilOffsetXCorrection="0" offsetX="0" offsetY="0" offsetXStack="-0.324310 -0.486465 -0.486465 -0.486465" offsetXShift="0" offsetYShift="0" TrapXAndScaleWithOffsetX="1" tiltZ="0" scale="0" trapezoidX="0" trapezoidY="0" shearing="0" aspectRatio="0" /> <ring_filter ringEnabled="0" kernel="5" falseDetectionLevel="0" defectEnabled="0" replaceStatic="0" replaceDefect="0" replaceNoisy="0" selectionEnabled="0" /> <beam black_level="0" clamp="2000" shift="0" spectrum="1.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000" /> <phase_filter enabled="0" type="MANU" manu_para="0.049999997" vector_auto="" /> <spot horizontal="0.000000 0.187500 0.125000 0.000000 -0.125000 -0.375000 -0.187500 -0.187500 0.062500 0.125000 0.062500 0.062500 -0.125000 -0.250000 0.000000 -0.125000 0.062500 0.125000 0.062500 0.062500 -0.125000 -0.125000 -0.125000 -0.062500 0.250000 0.125000 0.125000 0.125000 0.000000 -0.062500 0.312500 0.187500 0.000000" vertical="0.500000 0.437500 0.500000 0.500000 0.500000 0.437500 0.375000 0.312500 0.312500 0.250000 0.250000 0.250000 0.250000 0.187500 0.250000 0.187500 0.125000 0.125000 0.125000 0.125000 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.062500 0.000000 0.000000" scale="0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000" /> <projections cropright="0" cropleft="0" start="0" end="9599" step="1" applyLimitedAngleCorrection="false" unused="" /> <imagej scriptPath="" parameters="" imagesProcessed="0" baseProjName="" modifiedProjName="" /> <remove_outliers isEnabled="0" radius="1" threshold="100" outlierType="2" /> <scatter enabled="0" threshold="32000" parameters="" bottomLimit="-1" offset="0" automaticThreshold="0" /> <mesh enabled="0" ratio="0" advancedEnabled="0" gradientThreshold="12" filterSize="8" /> </conebeam>

[AnderBiguri]

Sure this helps. I can't run your code because it requires data, but what I can see is that in the xml file you have explicit dimensions of the geoemtry (eg voxel size) while in the code you compute them. This doesn't meant your computing is wrong of course, but do verify that the values are correct. Minor errors in geometry definition always lead to catastrophic failure in reconstruction.

[Resinato]

super happy!
-the width-length of the detector were simply swapped.
-the script was using the first 300 pictures (not one every 8..)

nowthe first of the 4 tiles is now something meaningful:

next step: reconstruct all the 4 tiles in a single volume. let's see how it goes. Thanks for your kind reply!