I find myself wanting to run through a list of (x,y,z) coordinates of some data volume (here called “d”) to do some sort of processing on each voxel. What I have come up with is the following…

First, find the set of coordinates that match some criterion. For example, find all coordinates in “d” that are greater than the 70^th percentile:
[cc lang="python"]coords = array( nonzero( d > prctile( d, 70 ) ) ).transpose()[/cc]
Now that we have the list of coordinates, we can run through each coordinate and do some sort of processing on it:
[cc lang="python"]
for ii,coord in enumerate( coords ):
r = coord[0]
c = coord[1]

# more stuff here
[/cc]
Obviously if you are using a 3 dimensional volume “d” then you would use:
[cc lang="python"]
s = coord[0]
r = coord[2]
c = coord[2]
[/cc]

Below is a Python class that will read in a Varian FDF file, or a Varian “.img” directory (which contains the FDF files). I have used this in the past, but can’t make any claims about it. I offer it up in hopes it is useful to someone.

[cc lang="python"]
import os
import re
from numpy import *
import struct

class Varian:

def __init__(self):
pass

def read( self, filename ):
if filename.endswith(‘.fdf’):
data = self.readFDF( filename )
elif filename.endswith(‘.img’):
data = self.readIMG( filename )
else:
print “Unknown filename %s ” % (filename)

return data

def readFDF(self, filename ):

fp = open( filename, ‘rb’ )

xsize = -1
ysize = -1
zsize = 1
bigendian = -1
done = False

while not done :

line = fp.readline()

if( len( line ) >= 1 and line[0] == chr(12) ):
break

if( len( line ) >= 1 and line[0] != chr(12) ):

if( line.find(‘bigendian’) > 0 ):
endian = line.split(‘=’)[-1].rstrip(‘\n; ‘).strip(‘ ‘)

if( line.find(‘echos’) > 0 ):
nechoes = line.split(‘=’)[-1].rstrip(‘\n; ‘).strip(‘ ‘)

if( line.find(‘echo_no’) > 0 ):
echo_no = line.split(‘=’)[-1].rstrip(‘\n; ‘).strip(‘ ‘)

if( line.find(‘nslices’) > 0 ):
nslices = line.split(‘=’)[-1].rstrip(‘\n; ‘).strip(‘ ‘)

if( line.find(‘slice_no’) > 0 ):
sl = line.split(‘=’)[-1].rstrip(‘\n; ‘).strip(‘ ‘)

if( line.find(‘matrix’) > 0 ):
m = re.findall(‘(\d+)’, line.rstrip())

if len(m) == 2:
xsize, ysize = int(m[0]), int(m[1])
elif len(m) == 3:
xsize, ysize, zsize = int(m[0]), int(m[1]), int(m[2])

fp.seek(-xsize*ysize*zsize*4,2)

if bigendian == 1:
fmt = “>%df” % (xsize*ysize*zsize)
else:
fmt = “<%df” % (xsize*ysize*zsize)

data = struct.unpack(fmt, fp.read(xsize*ysize*zsize*4))
data = array( data ).reshape( [xsize, ysize, zsize ] ).squeeze()

fp.close()

return data

def readIMG(self, directory):

# Get a list of all the FDF files in the directory
try:
files = os.listdir(directory)
except:
print “Could not find the directory %s” % directory
return

files = [ file for file in files if file.endswith('.fdf') ]

data = []
for file in files:
data.append( self.readFDF( directory+’/'+file ) )

data = transpose( array( data ), (1,2,0) )

return data

[/cc]

There is a great Python package pydicom that implements a nice interface in order to be able to access data within Dicom files.

One application which I wrote up was a dicom directory summarizer which goes through a list of dicom files and summarizes the types of MRI data in the directory.  I found myself getting frustrated trying to figure out which series of data was which given the huge number of dicom files (with really long names too!) in a directory.

The code below may be run within a Dicom directory and should run on Siemens Dicom data (IMA) files. It has been a while that I have run it so I can’t guarantee that it will work, but it should be a good place to start.

[cc lang="python"]
#! /usr/bin/python

import dicom
import os
import re

def blah(val):
return re.compile(‘[\-\w]+\.MR\.[\-\w]+\.\d+\.1\..*’).match(val, 1)

# Get a list of all the files
files = []
for entry in os.listdir(‘.’):
if ~os.path.isdir(entry) & entry.endswith(‘IMA’):
files.append(entry)

# Filter to find the first of each series
firsts = filter( blah, files )

firsts.sort(key=lambda s: int( re.compile(‘[\-\w]+\.MR\.[\-\w]+\.(\d+)\.1\..*’).search(s).group(1)) )

# Read the first and output some interesting stuff
d = dicom.ReadFile(firsts[1])
print ” Patient: ” + d.PatientsName
print “Acquired: ” + d.StudyDate[0:4]+”-”+d.StudyDate[4:6]+”-”+d.StudyDate[6:8] \
+ ” ” + d.StudyTime[0:2] + “:” + d.StudyTime[2:4] + “:” + d.StudyTime[4:6]
print “Comments: ” + d.ImageComments

# Run through the first file of each of the series
for entry in firsts:
d = dicom.ReadFile(entry)

num = re.compile(‘[\-\w]+\.MR\.[\-\w]+\.(\d+)\.1\..*’).search(entry).group(1)

out = “\t” + str(num) + “) ” + d.SeriesDescription

tt = ‘[_\-\w]+\.MR\.[_\-\w]+\.’+str(num)+’\..*’
count = 0
r = re.compile(tt)
for f in files:
if( r.match(f, 1) ):
#print “%s matches %d” % (f, ii)
count = count + 1

if( not re.compile(“.*(FA|TRACEW|TENSOR|ADC|MoCoSeries)$”).match(d.SeriesDescription, 1 ) ):
out += ” (vols=” + str(count)

if( ‘RepetitionTime’ in d ):
out += “, TR=” + str(d.RepetitionTime)

if( ‘EchoTime’ in d ):
out += “, TE=” + str(d.EchoTime)

out += “)”

print out

[/cc]