DICOM Houdini to VDB (you will need to remap values, result is not normalized 0-1 density)

uses file size to determine if it is an image with a trasform that needs to be included.

import pydicom
import numpy as np
import os
import hou

def create_global_attrib(geo, name, value):
    """Safely create and set a global attribute"""
    if geo.findGlobalAttrib(name) is None:
        attrib = geo.addAttrib(hou.attribType.Global, name, value)
    geo.setGlobalAttribValue(name, value)

def load_dicom_volume():
    dicom_dir = r"C:\Users\Reference\Desktop\dcom\DICOM\PA000001\ST000001\SE000001"
    slices = []

    # Load DICOM files with size filtering
    for fname in os.listdir(dicom_dir):
        file_path = os.path.join(dicom_dir, fname)
        if os.path.getsize(file_path) < 10 * 1024:
            continue  # Skip small headers

        try:
            ds = pydicom.dcmread(file_path)
            if not hasattr(ds, 'ImagePositionPatient') or not hasattr(ds, 'ImageOrientationPatient'):
                continue
            slices.append(ds)
        except Exception as e:
            print(f"Skipped {fname}: {str(e)}")

    if not slices:
        raise ValueError("No valid DICOM slices found")

    # Validate consistent orientation
    ref_orientation = np.array(slices[0].ImageOrientationPatient, dtype=np.float64)
    for ds in slices[1:]:
        if not np.allclose(np.array(ds.ImageOrientationPatient, dtype=np.float64), ref_orientation):
            raise ValueError("Inconsistent slice orientations detected")

    # Calculate slice normal for sorting
    row_dir = ref_orientation[:3]
    col_dir = ref_orientation[3:]
    slice_normal = np.cross(row_dir, col_dir)

    # Sort slices
    slices.sort(key=lambda x: np.dot(
        np.array(x.ImagePositionPatient, dtype=np.float64),
        slice_normal
    ))

    # Create geometry and volume
    geo = hou.pwd().geometry()

    # Store metadata
    create_global_attrib(geo, "PatientID", str(slices[0].get("PatientID", "")))
    create_global_attrib(geo, "StudyDescription", str(slices[0].get("StudyDescription", "")))
    create_global_attrib(geo, "SeriesDescription", str(slices[0].get("SeriesDescription", "")))

    # Store position and orientation data
    pos_array = np.array([ds.ImagePositionPatient for ds in slices]).flatten()
    ori_array = np.array([ds.ImageOrientationPatient for ds in slices]).flatten()
    create_global_attrib(geo, "ImagePositionPatient", pos_array.tolist())
    create_global_attrib(geo, "ImageOrientationPatient", ori_array.tolist())

    # Get pixel spacing (validate consistency)
    if not hasattr(slices[0], 'PixelSpacing'):
        raise ValueError("DICOM slices missing PixelSpacing attribute")
    
    pixel_spacing = slices[0].PixelSpacing
    for ds in slices[1:]:
        if not np.allclose(ds.PixelSpacing, pixel_spacing):
            raise ValueError("Inconsistent pixel spacing between slices")

    dx = float(pixel_spacing[0])  # Column spacing
    dy = float(pixel_spacing[1])  # Row spacing

    # Calculate slice spacing from positions
    slice_positions = [np.array(ds.ImagePositionPatient) for ds in slices]
    if len(slice_positions) > 1:
        distances = [np.linalg.norm(slice_positions[i+1] - slice_positions[i]) 
                    for i in range(len(slice_positions)-1)]
        dz = float(np.mean(distances))
    else:
        dz = 1.0  # Default value when spacing can't be determined

    # Create 3D volume
    pixel_data = np.stack([ds.pixel_array for ds in slices], axis=-1)
    volume_data = pixel_data.transpose(1, 0, 2)  # Houdini's XYZ ordering

    # Convert to proper data type and flatten in Fortran order
    if volume_data.dtype != np.float32:
        volume_data = volume_data.astype(np.float32)
    voxel_values = volume_data.flatten(order='F').tolist()

    # Create volume
    volume = geo.createVolume(volume_data.shape[0], 
                            volume_data.shape[1], 
                            volume_data.shape[2])

    # Set voxels
    try:
        volume.setAllVoxels(voxel_values)
    except TypeError:
        # Fallback method if the above fails
        for i, val in enumerate(voxel_values):
            volume.setVoxel(i, float(val))

    # Create proper transform matrix
    transform = hou.Matrix4()
    
    # Set the scale components
    transform.setAt(0, 0, dx)  # X scale (column spacing)
    transform.setAt(1, 1, dy)  # Y scale (row spacing)
    transform.setAt(2, 2, dz)  # Z scale (slice spacing)
    
    # Set the translation from first slice position
    first_pos = slice_positions[0]
    transform.setAt(3, 0, float(first_pos[0]))  # X position
    transform.setAt(3, 1, float(first_pos[1]))  # Y position
    transform.setAt(3, 2, float(first_pos[2]))  # Z position
    
    volume.setTransform(transform)

    # Set volume visualization properties (modern Houdini method)
    # Create visualization attributes directly on the volume primitive
    volume.setDisplayMode(hou.volumeDisplayMode.Raytraced)
    volume.setVisualization(hou.volumeVisualization.Smoke)
    volume.setShadingMode(hou.volumeShadingMode.PhysicallyBased)
    
    # Set reasonable default rendering parameters
    volume.setUniformSampling(True)
    volume.setStepSize(0.5)
    volume.setStepSizeShadow(1.0)

    print(f"Successfully loaded {len(slices)} slices")
    print(f"Volume dimensions: {volume_data.shape}")
    print(f"Voxel spacing: {dx:.4f} × {dy:.4f} × {dz:.4f} mm")
    print(f"First slice position: {first_pos}")
    print("Metadata stored in global attributes")

# Execute the function
load_dicom_volume()