Hi,
I’m having trouble updating the direction on an existing 3d itkImage.
Basically I want to reset the origin to [0,0,0] and direction to eye(3).
This first step seems to work but then when I use TransformPhysicalPointToContinuousIndex it returns the wrong index position. Here is a minimal example reproducing the error.
from __future__ import print_function
import itk
import sys
import numpy as np
def saveDirection(image):
''' copy direction to numpy array '''
direction = np.zeros((3,3))
for i in range(9):
direction[i//3, i%3] = image.GetDirection().GetVnlMatrix().get(i//3, i%3)
return direction
def readImageStack(folder):
''' load a single dicom series from a folder and reset origin and direction '''
ImageType = itk.Image[itk.SS, 3]
namesGenerator = itk.GDCMSeriesFileNames.New()
namesGenerator.SetUseSeriesDetails(True)
namesGenerator.SetGlobalWarningDisplay(False)
namesGenerator.SetDirectory(folder)
seriesUID = namesGenerator.GetSeriesUIDs()
if len(seriesUID) != 1:
print("Could not load images in folder %s. No image or multiple series."%folder)
sys.exit(1)
uid = seriesUID[0]
fileNames = namesGenerator.GetFileNames(uid)
reader = itk.ImageSeriesReader[ImageType].New()
dicomIO = itk.GDCMImageIO.New()
reader.SetImageIO(dicomIO)
reader.SetFileNames(fileNames)
reader.ForceOrthogonalDirectionOff()
reader.Update()
image = reader.GetOutput()
oldDirection = saveDirection(image)
print("old direction :\n", oldDirection)
image.SetOrigin((0,0,0)) # reset origin : works
# reset direction
np_dir_vnl = itk.GetVnlMatrixFromArray(np.eye(3))
direction = image.GetDirection()
direction.GetVnlMatrix().copy_in(np_dir_vnl.data_block())
image.SetDirection(direction)
# old code, gives the same bad result
#tmp = image.GetDirection().GetVnlMatrix()
#for i in range(9):
# tmp.set(i//3, i%3, 0)
#for i in range(3):
# tmp.set(i, i, 1)
#direction = image.GetDirection()
#direction.SetVnlMatrix(tmp)
#image.SetDirection(direction)
return image, oldDirection
if __name__ == "__main__":
DATA = "F:/Data/Images/Pat_John Doe/Study_1.3.12.2.1107.5.4.5.164089.30000018061106463446200000037/Series_1.3.12.2.1107.5.4.5.164089.30000018061106500168800005613"
image, oldDirection = readImageStack(DATA)
print(image)
# prob intensity at a given 3d physical point
p = np.array([11.44909474, 23.79051299, 40.75303588])
# since Direction is now identity and origin is [0,0,0], continuous index should be p/spacing
spacing = image.GetSpacing()
print("expected continuous index", [c/s for c,s in zip(p, spacing)])
# using TransformPhysicalPointToContinuousIndex -> off by a few voxels
pixelIndex = image.TransformPhysicalPointToContinuousIndex(p)
print(pixelIndex)
#using saved old direction :
#S-1*D-1*p
print("using old direction", np.linalg.inv(np.diag(spacing)).dot(oldDirection.T).dot(p))
The code above gives the following output :
old direction :
[[ 0.99998595 0.00174611 -0.00500522]
[-0.00175405 0.99999721 -0.00158081]
[ 0.00500245 0.00158957 0.99998622]]
Image (0000020CF83CC080)
RTTI typeinfo: class itk::Image<short,3>
Reference Count: 1
Modified Time: 326554
Debug: Off
Object Name:
Observers:
none
Source: (none)
Source output name: (none)
Release Data: Off
Data Released: False
Global Release Data: Off
PipelineMTime: 310064
UpdateMTime: 326552
RealTimeStamp: 0 seconds
LargestPossibleRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
BufferedRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
RequestedRegion:
Dimension: 3
Index: [0, 0, 0]
Size: [512, 512, 469]
Spacing: [0.108077, 0.108077, 0.108077]
Origin: [0, 0, 0]
Direction:
1 0 0
0 1 0
0 0 1
IndexToPointMatrix:
0.108076 0.000188715 -0.000540951
-0.000189573 0.108077 -0.00017085
0.000540651 0.000171796 0.108076
PointToIndexMatrix:
9.25251 -0.0162296 0.0462858
0.0161562 9.25261 0.0147077
-0.0463115 -0.0146267 9.25251
Inverse Direction:
0.999986 -0.00175405 0.00500245
0.00174611 0.999997 0.00158957
-0.00500522 -0.00158081 0.999986
PixelContainer:
ImportImageContainer (0000020D24629F90)
RTTI typeinfo: class itk::ImportImageContainer<unsigned __int64,short>
Reference Count: 1
Modified Time: 310132
Debug: Off
Object Name:
Observers:
none
Pointer: 0000020D09DF4040
Container manages memory: true
Size: 122945536
Capacity: 122945536
expected continuous index [105.93429698107786, 220.1249378791367, 377.0729742247302]
itkContinuousIndexD3 ([107.433, 220.909, 376.19])
using old direction [107.43298801 220.9086811 376.18957873]
We can see that origin and direction of the image are properly set to the expected values. But somehow IndexToPointMatrix, PointToIndexMatrix and Inverse Direction are not.
This then leads to incorrect behavior of image.TransformPhysicalPointToContinuousIndex() as it uses this out-of-date matrices.
Is this a bug or am I doing something wrong when updating the direction ? I did not try this in C++ but would the result be the same ?
I’m using itk for python version 5.0b3.
