Probably the simplest would be to get the volume as a numpy array and then write out the voxel indices and values you are interested in. Something like this (this outputs all voxels that has value >260):
import numpy as np voxelArray = slicer.util.arrayFromVolume(getNode('MRHead')) indices = np.where(voxelArray>260) numberOfVoxels = len(indices[0]) for pointIndex in range(numberOfVoxels): i = indices[0][pointIndex] j = indices[1][pointIndex] k = indices[2][pointIndex] print("%d %d %d %d" % (i, j, k, voxelArray[i,j,k]))
See some more examples in the script repository: https://www.slicer.org/wiki/Documentation/Nightly/ScriptRepository#Modify_voxels_in_a_volume
If you need physical (RAS) coordinates instead of voxel coordinates (IJK) then you have to multiply the (x,y,z,1) vector from the left with the volume’s IJKToRAS matrix.