OiO.lk Community platform!

Oio.lk is an excellent forum for developers, providing a wide range of resources, discussions, and support for those in the developer community. Join oio.lk today to connect with like-minded professionals, share insights, and stay updated on the latest trends and technologies in the development field.
  You need to log in or register to access the solved answers to this problem.
  • You have reached the maximum number of guest views allowed
  • Please register below to remove this limitation

Issue with Axial and Coronal Planes Display in MPR App using VTK and SimpleITK - DICOM CT

  • Thread starter Thread starter Miguel Nobre Menezes
  • Start date Start date
M

Miguel Nobre Menezes

Guest
Question: I'm developing a Multi-Planar Reconstruction (MPR) application using VTK and SimpleITK to display DICOM images in axial, sagittal, and coronal planes. However, I'm encountering issues with the display of axial and coronal planes. The images do not appear correctly in these orientations.

Problem Description:

  • The axial plane is sort of tilted/rotated but the actual slicing seems correct. The same happens in the sagital plane.
  • The coronal plane display correctly

Here are the relevant parts of my code:

Code for VTK Image Viewer Setup:

Code:
class VTKImageViewer:
    def __init__(self, parent, orientation='axial'):
        self.orientation = orientation
        self.vtkWidget = QVTKRenderWindowInteractor(parent)
        self.renderer = vtk.vtkRenderer()
        self.vtkWidget.GetRenderWindow().AddRenderer(self.renderer)
        self.interactor = self.vtkWidget.GetRenderWindow().GetInteractor()
        
        self.imageSlice = vtk.vtkImageSlice()
        self.sliceMapper = vtk.vtkImageSliceMapper()
        self.imageSlice.SetMapper(self.sliceMapper)
        
        if orientation == 'axial':
            self.sliceMapper.SetOrientationToZ()
        elif orientation == 'sagittal':
            self.sliceMapper.SetOrientationToX()
        elif orientation == 'coronal':
            self.sliceMapper.SetOrientationToY()
        
        self.renderer.AddViewProp(self.imageSlice)
        
        self.setup_camera()
        self.add_orientation_marker()

    def setup_camera(self):
        camera = self.renderer.GetActiveCamera()
        if self.orientation == 'axial':
            camera.SetViewUp(0, -1, 0)  # Changed from (1, 0, 0)
            camera.SetPosition(0, 0, -1)
        elif self.orientation == 'sagittal':
            camera.SetViewUp(0, 0, 1)
            camera.SetPosition(1, 0, 0)
        elif self.orientation == 'coronal':
            camera.SetViewUp(0, 0, 1)
            camera.SetPosition(0, 1, 0)
        self.renderer.ResetCamera()

    def add_orientation_marker(self):
        axes = vtk.vtkAxesActor()
        widget = vtk.vtkOrientationMarkerWidget()
        widget.SetOutlineColor(0.9300, 0.5700, 0.1300)
        widget.SetOrientationMarker(axes)
        widget.SetInteractor(self.interactor)
        widget.SetViewport(0.0, 0.0, 0.4, 0.4)
        widget.SetEnabled(1)
        widget.InteractiveOff()

    def setInputData(self, image_data):
        self.sliceMapper.SetInputData(image_data)
        self.setup_camera()
        self.renderer.ResetCamera()
        self.vtkWidget.GetRenderWindow().Render()

        print(f"--- {self.orientation.capitalize()} Viewer Input Data ---")
        print(f"Image dimensions: {image_data.GetDimensions()}")
        print(f"Image spacing: {image_data.GetSpacing()}")
        print(f"Image origin: {image_data.GetOrigin()}")

    def setSlice(self, slice_number):
        self.sliceMapper.SetSliceNumber(slice_number)
        self.vtkWidget.GetRenderWindow().Render()

    def get_slice_number(self):
        return self.sliceMapper.GetSliceNumber()

    def reset_camera(self):
        self.renderer.ResetCamera()
        self.vtkWidget.GetRenderWindow().Render()

    def set_window_level(self, window, level):
        self.imageSlice.GetProperty().SetColorWindow(window)
        self.imageSlice.GetProperty().SetColorLevel(level)
        self.vtkWidget.GetRenderWindow().Render()

Code for Converting SimpleITK Image to VTK:

Code:
def sitk_to_vtk(self, sitk_image):
    sitk_array = sitk.GetArrayFromImage(sitk_image)

    vtk_image = vtk.vtkImageData()
    vtk_image.SetDimensions(sitk_image.GetSize())
    vtk_image.SetSpacing(sitk_image.GetSpacing())
    vtk_image.SetOrigin(sitk_image.GetOrigin())

    if sitk_array.dtype == np.float64:
        dataType = vtk.VTK_DOUBLE
    elif sitk_array.dtype == np.float32:
        dataType = vtk.VTK_FLOAT
    elif sitk_array.dtype == np.int16:
        dataType = vtk.VTK_SHORT
    else:
        raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

    vtk_image.AllocateScalars(dataType, 1)
    
    vtk_array = vtk_image.GetPointData().GetScalars()
    vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

    transform = vtk.vtkTransform()
    direction = sitk_image.GetDirection()
    matrix = vtk.vtkMatrix4x4()
    for i in range(3):
        for j in range(3):
            matrix.SetElement(i, j, direction[i*3 + j])
    transform.SetMatrix(matrix)

    imageReslice = vtk.vtkImageReslice()
    imageReslice.SetInputData(vtk_image)
    imageReslice.SetResliceTransform(transform)
    imageReslice.SetInterpolationModeToLinear()
    imageReslice.SetOutputDimensionality(3)
    imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
    imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
    imageReslice.SetOutputExtent(vtk_image.GetExtent())
    imageReslice.Update()

    return imageReslice.GetOutput()

Main Application Code:

Code:
class MainWindow(QMainWindow):
    def __init__(self):
        super().__init__()
        uic.loadUi('main.ui', self)
        self.image = None
        self.dicom_series = {}
        self.pixel_spacing = None
        self.z_spacing = None
        self.size = None

        # Create VTK viewers for displaying images
        self.axialViewer = VTKImageViewer(self.axialContainer, 'axial')
        self.sagittalViewer = VTKImageViewer(self.sagittalContainer, 'sagittal')
        self.coronalViewer = VTKImageViewer(self.coronalContainer, 'coronal')

        # Create layouts for the containers
        self.axialLayout = QVBoxLayout(self.axialContainer)
        self.sagittalLayout = QVBoxLayout(self.sagittalContainer)
        self.coronalLayout = QVBoxLayout(self.coronalContainer)

        # Add VTK widgets to the layouts
        self.axialLayout.addWidget(self.axialViewer.vtkWidget)
        self.sagittalLayout.addWidget(self.sagittalViewer.vtkWidget)
        self.coronalLayout.addWidget(self.coronalViewer.vtkWidget)

        self.setup_connections()

    def setup_connections(self):
        self.loadButton.clicked.connect(self.load_dicom_folder)
        self.axialScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.sagittalScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.coronalScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.seriesList.currentRowChanged.connect(self.series_selected)

    def load_dicom_folder(self):
        folder = QFileDialog.getExistingDirectory(self, "Select DICOM folder")
        if folder:
            try:
                dicom_files = self.find_dicom_files(folder)
                
                if not dicom_files:
                    raise RuntimeError("No DICOM files found in the selected folder or its subfolders")
                
                series_dict = self.group_dicom_by_series(dicom_files)
                
                if not series_dict:
                    raise RuntimeError("No valid DICOM series found in the selected folder or its subfolders")
                
                self.dicom_series = series_dict
                self.seriesList.clear()
                
                total_series = len(series_dict)
                for idx, (series_uid, files) in enumerate(series_dict.items(), start=1):
                    ds = pydicom.dcmread(files[0])
                    description = ds.SeriesDescription if 'SeriesDescription' in ds else 'No Description'
                    item = QListWidgetItem(f"Series {idx}/{total_series} - <p>Question:
I'm developing a Multi-Planar Reconstruction (MPR) application using VTK and SimpleITK to display DICOM images in axial, sagittal, and coronal planes. However, I'm encountering issues with the display of axial and coronal planes. The images do not appear correctly in these orientations.</p>
<p>Problem Description:</p>
<ul>
<li>The axial plane is sort of tilted/rotated but the actual slicing seems correct. The same happens in the sagital plane.</li>
<li>The coronal plane display correctly</li>
</ul>
<p>Here are the relevant parts of my code:</p>
<p>Code for VTK Image Viewer Setup:</p>
<pre><code>class VTKImageViewer:
    def __init__(self, parent, orientation='axial'):
        self.orientation = orientation
        self.vtkWidget = QVTKRenderWindowInteractor(parent)
        self.renderer = vtk.vtkRenderer()
        self.vtkWidget.GetRenderWindow().AddRenderer(self.renderer)
        self.interactor = self.vtkWidget.GetRenderWindow().GetInteractor()
        
        self.imageSlice = vtk.vtkImageSlice()
        self.sliceMapper = vtk.vtkImageSliceMapper()
        self.imageSlice.SetMapper(self.sliceMapper)
        
        if orientation == 'axial':
            self.sliceMapper.SetOrientationToZ()
        elif orientation == 'sagittal':
            self.sliceMapper.SetOrientationToX()
        elif orientation == 'coronal':
            self.sliceMapper.SetOrientationToY()
        
        self.renderer.AddViewProp(self.imageSlice)
        
        self.setup_camera()
        self.add_orientation_marker()

    def setup_camera(self):
        camera = self.renderer.GetActiveCamera()
        if self.orientation == 'axial':
            camera.SetViewUp(0, -1, 0)  # Changed from (1, 0, 0)
            camera.SetPosition(0, 0, -1)
        elif self.orientation == 'sagittal':
            camera.SetViewUp(0, 0, 1)
            camera.SetPosition(1, 0, 0)
        elif self.orientation == 'coronal':
            camera.SetViewUp(0, 0, 1)
            camera.SetPosition(0, 1, 0)
        self.renderer.ResetCamera()

    def add_orientation_marker(self):
        axes = vtk.vtkAxesActor()
        widget = vtk.vtkOrientationMarkerWidget()
        widget.SetOutlineColor(0.9300, 0.5700, 0.1300)
        widget.SetOrientationMarker(axes)
        widget.SetInteractor(self.interactor)
        widget.SetViewport(0.0, 0.0, 0.4, 0.4)
        widget.SetEnabled(1)
        widget.InteractiveOff()

    def setInputData(self, image_data):
        self.sliceMapper.SetInputData(image_data)
        self.setup_camera()
        self.renderer.ResetCamera()
        self.vtkWidget.GetRenderWindow().Render()

        print(f"--- {self.orientation.capitalize()} Viewer Input Data ---")
        print(f"Image dimensions: {image_data.GetDimensions()}")
        print(f"Image spacing: {image_data.GetSpacing()}")
        print(f"Image origin: {image_data.GetOrigin()}")

    def setSlice(self, slice_number):
        self.sliceMapper.SetSliceNumber(slice_number)
        self.vtkWidget.GetRenderWindow().Render()

    def get_slice_number(self):
        return self.sliceMapper.GetSliceNumber()

    def reset_camera(self):
        self.renderer.ResetCamera()
        self.vtkWidget.GetRenderWindow().Render()

    def set_window_level(self, window, level):
        self.imageSlice.GetProperty().SetColorWindow(window)
        self.imageSlice.GetProperty().SetColorLevel(level)
        self.vtkWidget.GetRenderWindow().Render()
</code></pre>
<p>Code for Converting SimpleITK Image to VTK:</p>
<pre><code>def sitk_to_vtk(self, sitk_image):
    sitk_array = sitk.GetArrayFromImage(sitk_image)

    vtk_image = vtk.vtkImageData()
    vtk_image.SetDimensions(sitk_image.GetSize())
    vtk_image.SetSpacing(sitk_image.GetSpacing())
    vtk_image.SetOrigin(sitk_image.GetOrigin())

    if sitk_array.dtype == np.float64:
        dataType = vtk.VTK_DOUBLE
    elif sitk_array.dtype == np.float32:
        dataType = vtk.VTK_FLOAT
    elif sitk_array.dtype == np.int16:
        dataType = vtk.VTK_SHORT
    else:
        raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

    vtk_image.AllocateScalars(dataType, 1)
    
    vtk_array = vtk_image.GetPointData().GetScalars()
    vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

    transform = vtk.vtkTransform()
    direction = sitk_image.GetDirection()
    matrix = vtk.vtkMatrix4x4()
    for i in range(3):
        for j in range(3):
            matrix.SetElement(i, j, direction[i*3 + j])
    transform.SetMatrix(matrix)

    imageReslice = vtk.vtkImageReslice()
    imageReslice.SetInputData(vtk_image)
    imageReslice.SetResliceTransform(transform)
    imageReslice.SetInterpolationModeToLinear()
    imageReslice.SetOutputDimensionality(3)
    imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
    imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
    imageReslice.SetOutputExtent(vtk_image.GetExtent())
    imageReslice.Update()

    return imageReslice.GetOutput()
</code></pre>
<p>Main Application Code:</p>
<pre><code>class MainWindow(QMainWindow):
    def __init__(self):
        super().__init__()
        uic.loadUi('main.ui', self)
        self.image = None
        self.dicom_series = {}
        self.pixel_spacing = None
        self.z_spacing = None
        self.size = None

        # Create VTK viewers for displaying images
        self.axialViewer = VTKImageViewer(self.axialContainer, 'axial')
        self.sagittalViewer = VTKImageViewer(self.sagittalContainer, 'sagittal')
        self.coronalViewer = VTKImageViewer(self.coronalContainer, 'coronal')

        # Create layouts for the containers
        self.axialLayout = QVBoxLayout(self.axialContainer)
        self.sagittalLayout = QVBoxLayout(self.sagittalContainer)
        self.coronalLayout = QVBoxLayout(self.coronalContainer)

        # Add VTK widgets to the layouts
        self.axialLayout.addWidget(self.axialViewer.vtkWidget)
        self.sagittalLayout.addWidget(self.sagittalViewer.vtkWidget)
        self.coronalLayout.addWidget(self.coronalViewer.vtkWidget)

        self.setup_connections()

    def setup_connections(self):
        self.loadButton.clicked.connect(self.load_dicom_folder)
        self.axialScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.sagittalScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.coronalScrollBar.valueChanged.connect(lambda value: self.update_views())
        self.seriesList.currentRowChanged.connect(self.series_selected)

    def load_dicom_folder(self):
        folder = QFileDialog.getExistingDirectory(self, "Select DICOM folder")
        if folder:
            try:
                dicom_files = self.find_dicom_files(folder)
                
                if not dicom_files:
                    raise RuntimeError("No DICOM files found in the selected folder or its subfolders")
                
                series_dict = self.group_dicom_by_series(dicom_files)
                
                if not series_dict:
                    raise RuntimeError("No valid DICOM series found in the selected folder or its subfolders")
                
                self.dicom_series = series_dict
                self.seriesList.clear()
                
                total_series = len(series_dict)
                for idx, (series_uid, files) in enumerate(series_dict.items(), start=1):
                    ds = pydicom.dcmread(files[0])
                    description = ds.SeriesDescription if 'SeriesDescription' in ds else 'No Description'
                    item = QListWidgetItem(f"Series {idx}/{total_series} - {description} ({len(files)} images)")
                    item.setData(Qt.UserRole, series_uid)
                    self.seriesList.addItem(item)

                if self.seriesList.count() > 0:
                    self.seriesList.setCurrentRow(0)
                    self.series_selected(0)
            except Exception as e:
                QMessageBox.critical(self, "Error", f"Error loading DICOM: {str(e)}")
                print(f"Error loading DICOM: {str(e)}")

    def find_dicom_files(self, root_dir):
        dicom_files = []
        for root, _, files in os.walk(root_dir):
            for file in files:
                if file.lower().endswith('.dcm'):
                    full_path = os.path.join(root, file)
                    if self.is_dicom_image(full_path):
                        dicom_files.append(full_path)
        return dicom_files

    def is_dicom_image(self, file_path):
        try:
            ds = pydicom.dcmread(file_path)
            return 'PixelData' in ds
        except:
            return False

    def group_dicom_by_series(self, dicom_files):
        series_dict = {}
        for file in dicom_files:
            try:
                ds = pydicom.dcmread(file)
                series_id = ds.SeriesInstanceUID
                if series_id not in series_dict:
                    series_dict[series_id] = []
                series_dict[series_id].append(file)
            except Exception as e:
                print(f"Warning: Could not read file {file} with pydicom. Error: {str(e)}")
        return series_dict

    def series_selected(self, index):
        series_uid = self.seriesList.item(index).data(Qt.UserRole)
        self.load_series(self.dicom_series[series_uid])

    def load_series(self, dicom_files):
        sorted_files = self.sort_dicom_files(dicom_files)
        print(f"Number of DICOM files: {len(sorted_files)}")

        if not sorted_files:
            raise RuntimeError("No valid DICOM files found in the series")

        reader = sitk.ImageSeriesReader()
        reader.SetFileNames(sorted_files)
        self.image = reader.Execute()

        print(f"Image size: {self.image.GetSize()}")
        print(f"Image spacing: {self.image.GetSpacing()}")
        print(f"Image origin: {self.image.GetOrigin()}")
        print(f"Image direction: {self.image.GetDirection()}")

        # Ensure the image is in LPS orientation
        self.image = sitk.DICOMOrient(self.image, 'LPS')

        print("--- After DICOM Orientation ---")
        print(f"Image size: {self.image.GetSize()}")
        print(f"Image spacing: {self.image.GetSpacing()}")
        print(f"Image origin: {self.image.GetOrigin()}")
        print(f"Image direction: {self.image.GetDirection()}")

        spacing = self.image.GetSpacing()
        self.pixel_spacing = spacing[:2]
        self.z_spacing = spacing[2]

        self.size = self.image.GetSize()

        self.axialScrollBar.setMaximum(self.size[2] - 1)
        self.sagittalScrollBar.setMaximum(self.size[0] - 1)
        self.coronalScrollBar.setMaximum(self.size[1] - 1)

        # Convert SimpleITK image to VTK image
        vtk_image = self.sitk_to_vtk(self.image)

        # Set the VTK image as input to the viewers
        self.axialViewer.setInputData(vtk_image)
        self.sagittalViewer.setInputData(vtk_image)
        self.coronalViewer.setInputData(vtk_image)

        self.update_views()

    def sort_dicom_files(self, dicom_files):
        valid_files = []
        for file in dicom_files:
            try:
                ds = pydicom.dcmread(file)
                position = tuple(map(float, ds.ImagePositionPatient))
                valid_files.append((file, position))
            except Exception as e:
                print(f"Warning: Could not read file {file}. Error: {str(e)}")
                continue

        return [file for file, _ in sorted(valid_files, key=lambda x: x[1][2])]

    def update_views(self):
        if self.image is None:
            return

        axial_slice = self.axialScrollBar.value()
        sagittal_slice = self.sagittalScrollBar.value()
        coronal_slice = self.coronalScrollBar.value()

        self.axialViewer.setSlice(axial_slice)
        self.sagittalViewer.setSlice(sagittal_slice)
        self.coronalViewer.setSlice(coronal_slice)

    def sitk_to_vtk(self, sitk_image):
        sitk_array = sitk.GetArrayFromImage(sitk_image)

        vtk_image = vtk.vtkImageData()
        vtk_image.SetDimensions(sitk_image.GetSize())
        vtk_image.SetSpacing(sitk_image.GetSpacing())
        vtk_image.SetOrigin(sitk_image.GetOrigin())

        if sitk_array.dtype == np.float64:
            dataType = vtk.VTK_DOUBLE
        elif sitk_array.dtype == np.float32:
            dataType = vtk.VTK_FLOAT
        elif sitk_array.dtype == np.int16:
            dataType = vtk.VTK_SHORT
        else:
            raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

        vtk_image.AllocateScalars(dataType, 1)
        
        vtk_array = vtk_image.GetPointData().GetScalars()
        vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

        transform = vtk.vtkTransform()
        direction = sitk_image.GetDirection()
        matrix = vtk.vtkMatrix4x4()
        for i in range(3):
            for j in range(3):
                matrix.SetElement(i, j, direction[i*3 + j])
        transform.SetMatrix(matrix)

        imageReslice = vtk.vtkImageReslice()
        imageReslice.SetInputData(vtk_image)
        imageReslice.SetResliceTransform(transform)
        imageReslice.SetInterpolationModeToLinear()
        imageReslice.SetOutputDimensionality(3)
        imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
        imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
        imageReslice.SetOutputExtent(vtk_image.GetExtent())
        imageReslice.Update()

        return imageReslice.GetOutput()
</code></pre>
<p>What I've Tried:
Adjusting the camera orientation for the axial and coronal views.
Ensuring the image is in LPS orientation before converting to VTK.
Verifying the DICOM files and their order.</p>
<p>Screenshots attached<a href="https://i.sstatic.net/Aw34IG8J.png" rel="nofollow noreferrer">enter image description here</a></p>
<p>Despite these efforts, the axial and coronal planes are still not displayed correctly.</p>
<p>Question:
What could be causing the incorrect display in the axial and coronal planes?
Are there any specific steps or configurations in VTK or SimpleITK that I might be missing?
How can I ensure the correct orientation and display of these planes?
Any help or suggestions would be greatly appreciated!</p> ({len(files)} images)")
                    item.setData(Qt.UserRole, series_uid)
                    self.seriesList.addItem(item)

                if self.seriesList.count() > 0:
                    self.seriesList.setCurrentRow(0)
                    self.series_selected(0)
            except Exception as e:
                QMessageBox.critical(self, "Error", f"Error loading DICOM: {str(e)}")
                print(f"Error loading DICOM: {str(e)}")

    def find_dicom_files(self, root_dir):
        dicom_files = []
        for root, _, files in os.walk(root_dir):
            for file in files:
                if file.lower().endswith('.dcm'):
                    full_path = os.path.join(root, file)
                    if self.is_dicom_image(full_path):
                        dicom_files.append(full_path)
        return dicom_files

    def is_dicom_image(self, file_path):
        try:
            ds = pydicom.dcmread(file_path)
            return 'PixelData' in ds
        except:
            return False

    def group_dicom_by_series(self, dicom_files):
        series_dict = {}
        for file in dicom_files:
            try:
                ds = pydicom.dcmread(file)
                series_id = ds.SeriesInstanceUID
                if series_id not in series_dict:
                    series_dict[series_id] = []
                series_dict[series_id].append(file)
            except Exception as e:
                print(f"Warning: Could not read file {file} with pydicom. Error: {str(e)}")
        return series_dict

    def series_selected(self, index):
        series_uid = self.seriesList.item(index).data(Qt.UserRole)
        self.load_series(self.dicom_series[series_uid])

    def load_series(self, dicom_files):
        sorted_files = self.sort_dicom_files(dicom_files)
        print(f"Number of DICOM files: {len(sorted_files)}")

        if not sorted_files:
            raise RuntimeError("No valid DICOM files found in the series")

        reader = sitk.ImageSeriesReader()
        reader.SetFileNames(sorted_files)
        self.image = reader.Execute()

        print(f"Image size: {self.image.GetSize()}")
        print(f"Image spacing: {self.image.GetSpacing()}")
        print(f"Image origin: {self.image.GetOrigin()}")
        print(f"Image direction: {self.image.GetDirection()}")

        # Ensure the image is in LPS orientation
        self.image = sitk.DICOMOrient(self.image, 'LPS')

        print("--- After DICOM Orientation ---")
        print(f"Image size: {self.image.GetSize()}")
        print(f"Image spacing: {self.image.GetSpacing()}")
        print(f"Image origin: {self.image.GetOrigin()}")
        print(f"Image direction: {self.image.GetDirection()}")

        spacing = self.image.GetSpacing()
        self.pixel_spacing = spacing[:2]
        self.z_spacing = spacing[2]

        self.size = self.image.GetSize()

        self.axialScrollBar.setMaximum(self.size[2] - 1)
        self.sagittalScrollBar.setMaximum(self.size[0] - 1)
        self.coronalScrollBar.setMaximum(self.size[1] - 1)

        # Convert SimpleITK image to VTK image
        vtk_image = self.sitk_to_vtk(self.image)

        # Set the VTK image as input to the viewers
        self.axialViewer.setInputData(vtk_image)
        self.sagittalViewer.setInputData(vtk_image)
        self.coronalViewer.setInputData(vtk_image)

        self.update_views()

    def sort_dicom_files(self, dicom_files):
        valid_files = []
        for file in dicom_files:
            try:
                ds = pydicom.dcmread(file)
                position = tuple(map(float, ds.ImagePositionPatient))
                valid_files.append((file, position))
            except Exception as e:
                print(f"Warning: Could not read file {file}. Error: {str(e)}")
                continue

        return [file for file, _ in sorted(valid_files, key=lambda x: x[1][2])]

    def update_views(self):
        if self.image is None:
            return

        axial_slice = self.axialScrollBar.value()
        sagittal_slice = self.sagittalScrollBar.value()
        coronal_slice = self.coronalScrollBar.value()

        self.axialViewer.setSlice(axial_slice)
        self.sagittalViewer.setSlice(sagittal_slice)
        self.coronalViewer.setSlice(coronal_slice)

    def sitk_to_vtk(self, sitk_image):
        sitk_array = sitk.GetArrayFromImage(sitk_image)

        vtk_image = vtk.vtkImageData()
        vtk_image.SetDimensions(sitk_image.GetSize())
        vtk_image.SetSpacing(sitk_image.GetSpacing())
        vtk_image.SetOrigin(sitk_image.GetOrigin())

        if sitk_array.dtype == np.float64:
            dataType = vtk.VTK_DOUBLE
        elif sitk_array.dtype == np.float32:
            dataType = vtk.VTK_FLOAT
        elif sitk_array.dtype == np.int16:
            dataType = vtk.VTK_SHORT
        else:
            raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

        vtk_image.AllocateScalars(dataType, 1)
        
        vtk_array = vtk_image.GetPointData().GetScalars()
        vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

        transform = vtk.vtkTransform()
        direction = sitk_image.GetDirection()
        matrix = vtk.vtkMatrix4x4()
        for i in range(3):
            for j in range(3):
                matrix.SetElement(i, j, direction[i*3 + j])
        transform.SetMatrix(matrix)

        imageReslice = vtk.vtkImageReslice()
        imageReslice.SetInputData(vtk_image)
        imageReslice.SetResliceTransform(transform)
        imageReslice.SetInterpolationModeToLinear()
        imageReslice.SetOutputDimensionality(3)
        imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
        imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
        imageReslice.SetOutputExtent(vtk_image.GetExtent())
        imageReslice.Update()

        return imageReslice.GetOutput()

What I've Tried: Adjusting the camera orientation for the axial and coronal views. Ensuring the image is in LPS orientation before converting to VTK. Verifying the DICOM files and their order.

Screenshots attachedenter image description here

Despite these efforts, the axial and coronal planes are still not displayed correctly.

Question: What could be causing the incorrect display in the axial and coronal planes? Are there any specific steps or configurations in VTK or SimpleITK that I might be missing? How can I ensure the correct orientation and display of these planes? Any help or suggestions would be greatly appreciated!
<p>Question:
I'm developing a Multi-Planar Reconstruction (MPR) application using VTK and SimpleITK to display DICOM images in axial, sagittal, and coronal planes. However, I'm encountering issues with the display of axial and coronal planes. The images do not appear correctly in these orientations.</p>
<p>Problem Description:</p>
<ul>
<li>The axial plane is sort of tilted/rotated but the actual slicing seems correct. The same happens in the sagital plane.</li>
<li>The coronal plane display correctly</li>
</ul>
<p>Here are the relevant parts of my code:</p>
<p>Code for VTK Image Viewer Setup:</p>
<pre><code>class VTKImageViewer:
def __init__(self, parent, orientation='axial'):
self.orientation = orientation
self.vtkWidget = QVTKRenderWindowInteractor(parent)
self.renderer = vtk.vtkRenderer()
self.vtkWidget.GetRenderWindow().AddRenderer(self.renderer)
self.interactor = self.vtkWidget.GetRenderWindow().GetInteractor()

self.imageSlice = vtk.vtkImageSlice()
self.sliceMapper = vtk.vtkImageSliceMapper()
self.imageSlice.SetMapper(self.sliceMapper)

if orientation == 'axial':
self.sliceMapper.SetOrientationToZ()
elif orientation == 'sagittal':
self.sliceMapper.SetOrientationToX()
elif orientation == 'coronal':
self.sliceMapper.SetOrientationToY()

self.renderer.AddViewProp(self.imageSlice)

self.setup_camera()
self.add_orientation_marker()

def setup_camera(self):
camera = self.renderer.GetActiveCamera()
if self.orientation == 'axial':
camera.SetViewUp(0, -1, 0) # Changed from (1, 0, 0)
camera.SetPosition(0, 0, -1)
elif self.orientation == 'sagittal':
camera.SetViewUp(0, 0, 1)
camera.SetPosition(1, 0, 0)
elif self.orientation == 'coronal':
camera.SetViewUp(0, 0, 1)
camera.SetPosition(0, 1, 0)
self.renderer.ResetCamera()

def add_orientation_marker(self):
axes = vtk.vtkAxesActor()
widget = vtk.vtkOrientationMarkerWidget()
widget.SetOutlineColor(0.9300, 0.5700, 0.1300)
widget.SetOrientationMarker(axes)
widget.SetInteractor(self.interactor)
widget.SetViewport(0.0, 0.0, 0.4, 0.4)
widget.SetEnabled(1)
widget.InteractiveOff()

def setInputData(self, image_data):
self.sliceMapper.SetInputData(image_data)
self.setup_camera()
self.renderer.ResetCamera()
self.vtkWidget.GetRenderWindow().Render()

print(f"--- {self.orientation.capitalize()} Viewer Input Data ---")
print(f"Image dimensions: {image_data.GetDimensions()}")
print(f"Image spacing: {image_data.GetSpacing()}")
print(f"Image origin: {image_data.GetOrigin()}")

def setSlice(self, slice_number):
self.sliceMapper.SetSliceNumber(slice_number)
self.vtkWidget.GetRenderWindow().Render()

def get_slice_number(self):
return self.sliceMapper.GetSliceNumber()

def reset_camera(self):
self.renderer.ResetCamera()
self.vtkWidget.GetRenderWindow().Render()

def set_window_level(self, window, level):
self.imageSlice.GetProperty().SetColorWindow(window)
self.imageSlice.GetProperty().SetColorLevel(level)
self.vtkWidget.GetRenderWindow().Render()
</code></pre>
<p>Code for Converting SimpleITK Image to VTK:</p>
<pre><code>def sitk_to_vtk(self, sitk_image):
sitk_array = sitk.GetArrayFromImage(sitk_image)

vtk_image = vtk.vtkImageData()
vtk_image.SetDimensions(sitk_image.GetSize())
vtk_image.SetSpacing(sitk_image.GetSpacing())
vtk_image.SetOrigin(sitk_image.GetOrigin())

if sitk_array.dtype == np.float64:
dataType = vtk.VTK_DOUBLE
elif sitk_array.dtype == np.float32:
dataType = vtk.VTK_FLOAT
elif sitk_array.dtype == np.int16:
dataType = vtk.VTK_SHORT
else:
raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

vtk_image.AllocateScalars(dataType, 1)

vtk_array = vtk_image.GetPointData().GetScalars()
vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

transform = vtk.vtkTransform()
direction = sitk_image.GetDirection()
matrix = vtk.vtkMatrix4x4()
for i in range(3):
for j in range(3):
matrix.SetElement(i, j, direction[i*3 + j])
transform.SetMatrix(matrix)

imageReslice = vtk.vtkImageReslice()
imageReslice.SetInputData(vtk_image)
imageReslice.SetResliceTransform(transform)
imageReslice.SetInterpolationModeToLinear()
imageReslice.SetOutputDimensionality(3)
imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
imageReslice.SetOutputExtent(vtk_image.GetExtent())
imageReslice.Update()

return imageReslice.GetOutput()
</code></pre>
<p>Main Application Code:</p>
<pre><code>class MainWindow(QMainWindow):
def __init__(self):
super().__init__()
uic.loadUi('main.ui', self)
self.image = None
self.dicom_series = {}
self.pixel_spacing = None
self.z_spacing = None
self.size = None

# Create VTK viewers for displaying images
self.axialViewer = VTKImageViewer(self.axialContainer, 'axial')
self.sagittalViewer = VTKImageViewer(self.sagittalContainer, 'sagittal')
self.coronalViewer = VTKImageViewer(self.coronalContainer, 'coronal')

# Create layouts for the containers
self.axialLayout = QVBoxLayout(self.axialContainer)
self.sagittalLayout = QVBoxLayout(self.sagittalContainer)
self.coronalLayout = QVBoxLayout(self.coronalContainer)

# Add VTK widgets to the layouts
self.axialLayout.addWidget(self.axialViewer.vtkWidget)
self.sagittalLayout.addWidget(self.sagittalViewer.vtkWidget)
self.coronalLayout.addWidget(self.coronalViewer.vtkWidget)

self.setup_connections()

def setup_connections(self):
self.loadButton.clicked.connect(self.load_dicom_folder)
self.axialScrollBar.valueChanged.connect(lambda value: self.update_views())
self.sagittalScrollBar.valueChanged.connect(lambda value: self.update_views())
self.coronalScrollBar.valueChanged.connect(lambda value: self.update_views())
self.seriesList.currentRowChanged.connect(self.series_selected)

def load_dicom_folder(self):
folder = QFileDialog.getExistingDirectory(self, "Select DICOM folder")
if folder:
try:
dicom_files = self.find_dicom_files(folder)

if not dicom_files:
raise RuntimeError("No DICOM files found in the selected folder or its subfolders")

series_dict = self.group_dicom_by_series(dicom_files)

if not series_dict:
raise RuntimeError("No valid DICOM series found in the selected folder or its subfolders")

self.dicom_series = series_dict
self.seriesList.clear()

total_series = len(series_dict)
for idx, (series_uid, files) in enumerate(series_dict.items(), start=1):
ds = pydicom.dcmread(files[0])
description = ds.SeriesDescription if 'SeriesDescription' in ds else 'No Description'
item = QListWidgetItem(f"Series {idx}/{total_series} - {description} ({len(files)} images)")
item.setData(Qt.UserRole, series_uid)
self.seriesList.addItem(item)

if self.seriesList.count() > 0:
self.seriesList.setCurrentRow(0)
self.series_selected(0)
except Exception as e:
QMessageBox.critical(self, "Error", f"Error loading DICOM: {str(e)}")
print(f"Error loading DICOM: {str(e)}")

def find_dicom_files(self, root_dir):
dicom_files = []
for root, _, files in os.walk(root_dir):
for file in files:
if file.lower().endswith('.dcm'):
full_path = os.path.join(root, file)
if self.is_dicom_image(full_path):
dicom_files.append(full_path)
return dicom_files

def is_dicom_image(self, file_path):
try:
ds = pydicom.dcmread(file_path)
return 'PixelData' in ds
except:
return False

def group_dicom_by_series(self, dicom_files):
series_dict = {}
for file in dicom_files:
try:
ds = pydicom.dcmread(file)
series_id = ds.SeriesInstanceUID
if series_id not in series_dict:
series_dict[series_id] = []
series_dict[series_id].append(file)
except Exception as e:
print(f"Warning: Could not read file {file} with pydicom. Error: {str(e)}")
return series_dict

def series_selected(self, index):
series_uid = self.seriesList.item(index).data(Qt.UserRole)
self.load_series(self.dicom_series[series_uid])

def load_series(self, dicom_files):
sorted_files = self.sort_dicom_files(dicom_files)
print(f"Number of DICOM files: {len(sorted_files)}")

if not sorted_files:
raise RuntimeError("No valid DICOM files found in the series")

reader = sitk.ImageSeriesReader()
reader.SetFileNames(sorted_files)
self.image = reader.Execute()

print(f"Image size: {self.image.GetSize()}")
print(f"Image spacing: {self.image.GetSpacing()}")
print(f"Image origin: {self.image.GetOrigin()}")
print(f"Image direction: {self.image.GetDirection()}")

# Ensure the image is in LPS orientation
self.image = sitk.DICOMOrient(self.image, 'LPS')

print("--- After DICOM Orientation ---")
print(f"Image size: {self.image.GetSize()}")
print(f"Image spacing: {self.image.GetSpacing()}")
print(f"Image origin: {self.image.GetOrigin()}")
print(f"Image direction: {self.image.GetDirection()}")

spacing = self.image.GetSpacing()
self.pixel_spacing = spacing[:2]
self.z_spacing = spacing[2]

self.size = self.image.GetSize()

self.axialScrollBar.setMaximum(self.size[2] - 1)
self.sagittalScrollBar.setMaximum(self.size[0] - 1)
self.coronalScrollBar.setMaximum(self.size[1] - 1)

# Convert SimpleITK image to VTK image
vtk_image = self.sitk_to_vtk(self.image)

# Set the VTK image as input to the viewers
self.axialViewer.setInputData(vtk_image)
self.sagittalViewer.setInputData(vtk_image)
self.coronalViewer.setInputData(vtk_image)

self.update_views()

def sort_dicom_files(self, dicom_files):
valid_files = []
for file in dicom_files:
try:
ds = pydicom.dcmread(file)
position = tuple(map(float, ds.ImagePositionPatient))
valid_files.append((file, position))
except Exception as e:
print(f"Warning: Could not read file {file}. Error: {str(e)}")
continue

return [file for file, _ in sorted(valid_files, key=lambda x: x[1][2])]

def update_views(self):
if self.image is None:
return

axial_slice = self.axialScrollBar.value()
sagittal_slice = self.sagittalScrollBar.value()
coronal_slice = self.coronalScrollBar.value()

self.axialViewer.setSlice(axial_slice)
self.sagittalViewer.setSlice(sagittal_slice)
self.coronalViewer.setSlice(coronal_slice)

def sitk_to_vtk(self, sitk_image):
sitk_array = sitk.GetArrayFromImage(sitk_image)

vtk_image = vtk.vtkImageData()
vtk_image.SetDimensions(sitk_image.GetSize())
vtk_image.SetSpacing(sitk_image.GetSpacing())
vtk_image.SetOrigin(sitk_image.GetOrigin())

if sitk_array.dtype == np.float64:
dataType = vtk.VTK_DOUBLE
elif sitk_array.dtype == np.float32:
dataType = vtk.VTK_FLOAT
elif sitk_array.dtype == np.int16:
dataType = vtk.VTK_SHORT
else:
raise ValueError(f"Unsupported data type: {sitk_array.dtype}")

vtk_image.AllocateScalars(dataType, 1)

vtk_array = vtk_image.GetPointData().GetScalars()
vtk_array.SetVoidArray(sitk_array.ravel(), sitk_array.size, 1)

transform = vtk.vtkTransform()
direction = sitk_image.GetDirection()
matrix = vtk.vtkMatrix4x4()
for i in range(3):
for j in range(3):
matrix.SetElement(i, j, direction[i*3 + j])
transform.SetMatrix(matrix)

imageReslice = vtk.vtkImageReslice()
imageReslice.SetInputData(vtk_image)
imageReslice.SetResliceTransform(transform)
imageReslice.SetInterpolationModeToLinear()
imageReslice.SetOutputDimensionality(3)
imageReslice.SetOutputOrigin(vtk_image.GetOrigin())
imageReslice.SetOutputSpacing(vtk_image.GetSpacing())
imageReslice.SetOutputExtent(vtk_image.GetExtent())
imageReslice.Update()

return imageReslice.GetOutput()
</code></pre>
<p>What I've Tried:
Adjusting the camera orientation for the axial and coronal views.
Ensuring the image is in LPS orientation before converting to VTK.
Verifying the DICOM files and their order.</p>
<p>Screenshots attached<a href="https://i.sstatic.net/Aw34IG8J.png" rel="nofollow noreferrer">enter image description here</a></p>
<p>Despite these efforts, the axial and coronal planes are still not displayed correctly.</p>
<p>Question:
What could be causing the incorrect display in the axial and coronal planes?
Are there any specific steps or configurations in VTK or SimpleITK that I might be missing?
How can I ensure the correct orientation and display of these planes?
Any help or suggestions would be greatly appreciated!</p>
 

Latest posts

Top