Hello
Dear Andras and Li
Thank you very much for your guidance and sorry for the late reply.
I want to use the particles system resulting from Chest Imaging Platform (CIP) in classification based on particles properties. These properties can indicate the features of the blood vessels of the lungs.
The output of CIP is in vtk format, which contains the following properties for each particle.
1- scale
2- val
3- eigenvalue
4- eigenvector
5- Hessian elements
6- Spatial coordinates
Comment: At first, I convert vtk formt to vtp format.
I want to transfer this information (particle properties) from the vtp (xml structure) to 3D matrices (the same size as the CT image matrix) or tensors. The end result will be, for example, files with the following names:
scaleMatrix.nrrd , valMatrix.nrrd, eignevalueMatrix.nrrd, eignevectorTensor.nrrd and …
These files must contain matrices or tensors in which the voxels, in the particle position, have values such as scale, val, eigenvalue, and … .
I think that the coordinate system of CT images in nrrd format is the same as the coordinate system in ITK that is LPS. Is it true?
Accordingly, I used the space origin item in the header section of the CT.nrrd file to create a transform matrix.
I enter the following commands:
sn@MP:~$ ipython
Python 3.6.9 (default, Oct 8 2020, 12:12:24)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.14.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: # Load modules
...: import numpy as np
...: import nrrd
...: import vtk
...: import os
In [2]: # Load CT (nrrd file)
...: filenameCT = '/home/sn/testCIP/result1/WholeLung01/wholeLung/CT_wholeLung.nrrd'
...: CT, headerCT = nrrd.read(filenameCT)
...: itemsHeaderCT = list(headerCT.items())
...: print("items in CT header")
...: for i in range(len(itemsHeaderCT)):
...: print(i, itemsHeaderCT[i][0], itemsHeaderCT[3][1])
...:
items in CT header
0 type [509 347 629]
1 dimension [509 347 629]
2 space [509 347 629]
3 sizes [509 347 629]
4 space directions [509 347 629]
5 kinds [509 347 629]
6 endian [509 347 629]
7 encoding [509 347 629]
8 space origin [509 347 629]
In [3]: print("items header of nnd file CT", itemsHeaderCT[3])
items header of nnd file CT ('sizes', array([509, 347, 629]))
In [4]: scaleMatrix = np.zeros((itemsHeaderCT[3][1]))
In [5]: print("matrix size scaleMatrix", scaleMatrix.shape)
matrix size scaleMatrix (509, 347, 629)
In [6]: print("matrix size CT:", CT.shape)
matrix size CT: (509, 347, 629)
In [7]: # Load particles system (vtp file)
...: filenameParticles = "/home/sn/testCIP/result1/WholeLung01/wholeLungVesselParticles.vtp";
...: readerDataVTP = vtk.vtkXMLPolyDataReader();
...: readerDataVTP.SetFileName(filenameParticles);
...: readerDataVTP.Update();
...: polydata = readerDataVTP.GetOutput();
In [8]: # Get particles coordinate
...: particlesCoordinates = polydata.GetPoints();
...: numberOfParticles = readerDataVTP.GetNumberOfPoints();
...: print("number of particles:", numberOfParticles)
...: print("Type of particles coordinates in vtp file: ", type(particlesCoordinates.GetPoint(0)))
number of particles: 48425
Type of particles coordinates in vtp file: <class 'tuple'>
In [9]: # Convert and get first particle coordinate from tuple to numpy array
...: particleCoordinates = np.asarray(particlesCoordinates.GetPoint(0))
...: print("The coordinates of first particle: ", particleCoordinates)
The coordinates of first particle: [-146.41410828 66.66109467 -923.09509277]
In [10]: # Convert and get end particle coordinate from tuple to numpy array
...: particleCoordinates = np.asarray(particlesCoordinates.GetPoint(numberOfParticles-1))
...: print("The coordinates of ", numberOfParticles, " particle: ", particleCoordinates)
The coordinates of 48425 particle: [ -74.92162323 13.57623863 -547.72540283]
In [11]: # Get scales particles
...: pointData = polydata.GetPointData()
...: scaleData = pointData.GetArray(0)
In [12]: # Get number of particles
...: print("The number of particles based on scale: ", scaleData.GetNumberOfTuples())
The number of particles based on scale: 48425
In [13]: # Get the scale of first particles
...: scaleParticle = scaleData.GetComponent(0,0)
In [14]: print("The scale of first particles: ", scaleParticle)
The scale of first particles: 1.5849934816360474
In [15]: # Get the scale of 273 particle
...: scaleParticle = scaleData.GetComponent(273,0)
...: print("The scale of 273 particle: ", scaleParticle)
The scale of 273 particle: 1.9540350437164307
In [16]: # Create scale matrix
...: scaleMatrix = np.zeros((itemsHeaderCT[3][1]))
...: print("The size of scale matrix: ", scaleMatrix.shape)
The size of scale matrix: (509, 347, 629)
In [17]: # Get space origin item in header nrrd CT file for creating transform matrix
...: spaceOrigin = np.array([itemsHeaderCT[8][1][0],itemsHeaderCT[8][1][1] ,itemsHeaderCT[8][1][2]])
...: rowIndex = np.int(np.round(np.absolute(spaceOrigin[0])))
...: columnIndex = np.int(np.round(np.absolute(spaceOrigin[1])))
...: sliceIndex = np.int(np.round(np.absolute(spaceOrigin[2])))
In [18]: # Create Transform matrix for converting the particle coordinates to voxel index in scale matrix
...: transferLat, transferVert, transferLong = ((1, 0, 0, columnIndex), (0, 1, 0, rowIndex), (0, 0, 1, sliceIndex))
...: transformMatrix = np.array([transferLat, transferVert, transferLong])
...: print("Transfer Matrix: ", '\n', transformMatrix)
Transfer Matrix:
[[ 1 0 0 118]
[ 0 1 0 190]
[ 0 0 1 925]]
In [19]: # Set voxel index
...: indexParticle = 1
...: for indexParticle in range(numberOfParticles):
...: particleCoordsTmp = np.asarray(particlesCoordinates.GetPoint(indexParticle))
...: particleCoordinates = np.array([[particleCoordsTmp[0], particleCoordsTmp[1], particleCoordsTmp[2], 1]]).T
...: voxelIndex = (np.absolute(np.round(transformMatrix.dot(particleCoordinates)))).astype(int)
...: scaleMatrix[np.int(voxelIndex[2,0])][np.int(voxelIndex[0,0])][np.int(voxelIndex[1,0])] = scaleData.GetComponent(indexParticle, 0)
...:
In [20]: # Save scale matrix
...: nrrd.write('scaleMatrix.nrrd', scaleMatrix, headerCT)
In [21]: exit
sn@MP:~$
After executing these commands, I noticed that the scale property of all particles has not been transferred into the scaleMatrix. Also, its directions are not the same as CT.
If possible, please guide me.
Best regards.
Shahrokh.
