Clinical Helical CT

This example reconstructs measured projections of patients from a clinical helical CT scanner. The same script also reconstructs case C002 (chest, Siemens, flying focal spot, 11072 views). The data are the full-dose projections of LDCT-and-Projection-data (see Data sources). Each view is one detector readout in the DICOM-CT-PD format.

Case L213 (abdomen, GE): 7140 views, 888 × 64 cylindrical detector, 1 mm voxels, 30 CGLS iterations without regularization. Soft-tissue window.

Geometry

The scanner geometry goes to Projector as per-view tensors, without rebinning to a circular orbit or a flat detector:

• Source: the focal spot of each view, including the flying focal spot offsets in angle, radius and z.

• Detector: a cylindrical surface centred on the nominal focal spot, given as a detector_surface callback. Projection and backprojection use these pixel positions directly.

• Views: the views are split over the GPUs with distributed=True; cgls() runs on the distributed operator.

The detector column of the central ray is 368.56 for C002 and 442.2 for L213, as 0-based stored columns. We estimated both values beforehand with gradient-based calibration of the column offset; that calibration is not part of this script.

docs/assets/gallery/compute.py

def helical(args):
    """CGLS of measured full-dose helical projections with diffct's distributed Projector."""
    from torch import distributed as dist
    device = torch.device("cuda", int(os.environ["LOCAL_RANK"]))
    torch.cuda.set_device(device)
    dist.init_process_group("nccl", device_id=device)
    rank = dist.get_rank()
    case = Path(args.case)
    m = dict(np.load(case / "meta.npz"))
    proj = np.load(case / "proj.npy", mmap_mode="r")
    (n_cols, n_rows), z = proj.shape[1:], m["z"]
    sdd, du, dv, rho = (float(m[k][0]) for k in ("sdd", "du", "dv", "rho"))
    margin = math.ceil(n_rows / 2 * dv * (rho + 250) / sdd + np.abs(m["dz"]).max())
    mid = 0.5 * (z.min() + z.max())
    depth = int(math.ceil(z.max() - z.min() + 2 * margin))
    off_u = args.center - (n_cols / 2 - 0.5)
    off_v = (float(m["ce_row"][0]) - 1) - (n_rows / 2 - 0.5)

    def cyl(r, phi, height):
        return torch.stack((-r * phi.sin(), r * phi.cos(), height), 1)

    t = lambda k: torch.as_tensor(m[k], dtype=torch.float64)  # noqa: E731
    phi, zz, r = t("phi"), t("z") - mid, t("rho")
    source = cyl(r + t("drho"), phi + t("dphi"), zz + t("dz"))  # flying focal spot
    s, c, zero = phi.sin(), phi.cos(), torch.zeros_like(phi)
    detector = cyl(r, phi, zz) + sdd * torch.stack((s, -c, zero), 1)
    trajectory = tuple(x.float() for x in (source, detector, torch.stack((c, s, zero), 1),
                                            torch.stack((zero, zero, torch.ones_like(phi)), 1)))

    def surface(u, v):  # cylinder centred on the nominal focal spot; stored rows run against v
        gamma = (u - off_u * du) / sdd
        return torch.stack((sdd * gamma.sin(), -(v - off_v * dv), sdd * (gamma.cos() - 1)), -1)

    operator = Projector(trajectory, (depth, 512, 512), (n_cols, n_rows), detector_spacing=(du, dv),
                         voxel_spacing=1.0, detector_surface=surface, distributed=True, partition="views")
    b = torch.from_numpy(np.ascontiguousarray(proj[operator.view_slice])).to(device).div_(MU_WATER)
    keep, snapshots = set(args.keep), {}

    def callback(iteration, residual, volume):
        if rank == 0 and iteration in keep:
            snapshots[iteration] = volume.detach().float().cpu()
            print(f"it {iteration} residual {residual:.6g}", flush=True)

    start = time.time()
    result = cgls(operator, b, max(args.keep), callback=callback)
    torch.cuda.synchronize()
    seconds = time.time() - start
    if rank == 0:
        axial_z = args.axial_z or [mid - 60.0, mid, mid + 60.0]
        axial_index = [int(round(zv - (mid - depth / 2))) for zv in axial_z]
        arrays = {}
        for it, volume in snapshots.items():
            arrays[f"axial_{it}"] = np.stack([volume[k].numpy().astype(np.float16) for k in axial_index])
            arrays[f"coronal_{it}"] = volume[:, 256].numpy().astype(np.float16)
            arrays[f"sagittal_{it}"] = volume[:, :, 256].numpy().astype(np.float16)
        history = np.array([float(x) for x in result.residual_history])
        save(Path(args.out), f"helical_{case.name}", **arrays, iterations=np.array(sorted(snapshots)),
             residual=history, axial_z=np.array(axial_z), seconds=seconds,
             views=len(z), detector=np.array([n_cols, n_rows]), gpus=dist.get_world_size(),
             manufacturer=str(m["manufacturer"]), depth=depth, z_mid=mid)
        print("HELICAL_DONE", case.name, f"{seconds:.1f} s", flush=True)
    dist.barrier()
    dist.destroy_process_group()

Run it

1. Download the full-dose projection and image series of a Chest or Liver case from The Cancer Imaging Archive.

2. Convert one case to arrays:

   pip install pydicom
   python docs/assets/gallery/ldct_prep.py <projection dir> <image dir> ldct/L213
   

3. Reconstruct on 8 GPUs:

   torchrun --nproc-per-node=8 docs/assets/gallery/compute.py helical \
       --case ldct/L213 --center 442.2 --axial-z -418 -364 -307 --out arrays