from collections import namedtuple from utils.path import path import numpy as np class Dimension(object): def __init__(self): self.scale = 1.0 self.x = self.y = self.z = 0 # dimensions of world (resolution) self.sx = self.sy = self.sz = 0. #voxel thickness in m self.tx = self.ty = self.tz = 0. # world center def __eq__ (self, other): ''' support for comparing Dimension equality''' return (self.scale == other.scale and self.x == other.x and self.y == other.y and self.z == other.z and self.sx == other.sx and self.sy == other.sy and self.sz == other.sz and self.tx == other.tx and self.ty == other.ty and self.tz == other.tz) def __neq__ (self, other): return not self.__eq__(other) Source = namedtuple('Source', 'x y z') class VolumeImage(object): #cdef public fname # name of .dat file #cdef public rawfile_name # name of .raw file #cdef public coeff_name # unknown file #cdef public mode # datatype of scalars - only 32bit floats atm #cdef public np.ndarray data # the 3d numpy data with the scalars #cdef public source # 3-tuple with source coordinates def __init__(self, fname=None): self.di = Dimension() self.data = np.zeros((0, 0, 0)) self.rawfile_name = "" self.coeff_name = "" self.mode = "UCHAR" self.source = Source(0, 0, 0) if fname is not None: self.read(fname) def getVoxel(self, x, y, z): if x < 0 or y < 0 or z < 0 or x >= self.di.x or y >= self.di.y or z >= self.di.z: return None else: return self.data[x][y][z] * self.di.scale def getPoint(self, x, y, z): return self.getVoxel(*self.translate(x, y, z)) def translate(self, x, y, z): '''translate float real world coordinates into int coordinates for space (self.di.x, self.di.y, self.di.z)''' return (int(round((x - self.di.tx) / self.di.sx)), int(round((y - self.di.ty) / self.di.sy)), int(round((z - self.di.tz) / self.di.sz))) def asDecibelData(self, clamp): dbdata = 10 * np.log10(self.data) # filter values that are -inf'ed by np.log(0) dbdata[np.isinf(dbdata)] = clamp return dbdata def read(self, fname, offset=1): ''' parse header file (.dat) into self.di and read the content of the data file (.raw) into self.data (a numpy array) ''' self.fname = path(fname).abspath() assert self.fname.ext == '.dat' with open(self.fname) as f: for line in f: buf = line.split() if buf[0] == "ObjectFileName:": self.rawfile_name = buf[1] elif buf[0] == "CoeffFileName:": self.coeff_name = buf[1] elif buf[0] == "Resolution:": self.di.x = int(buf[1]) self.di.y = int(buf[2]) self.di.z = int(buf[3]) elif buf[0] == "SliceThickness:": self.di.sx = float(buf[1]) self.di.sy = float(buf[2]) self.di.sz = float(buf[3]) elif buf[0] == "Origin:": self.di.tx = float(buf[1]) self.di.ty = float(buf[2]) self.di.tz = float(buf[3]) elif buf[0] == "Format:": self.mode = buf[1] elif buf[0] == "Scale:": self.di.scale = float(buf[1]) elif buf[0] == "Source": # like Source 0: 0.000000 0.000000 0.000000 self.source = Source(float(buf[2]), float(buf[3]), float(buf[4])) datafile = self.fname.parent.joinpath(self.rawfile_name) if datafile.namebase != self.fname.namebase: datafile = self.fname.parent.joinpath(self.fname.namebase + '.raw') print 'adapting datafile from to %s' % (datafile.name) assert datafile.exists(), '%s does not exists' % datafile assert self.mode == "FLOAT" shape = (self.di.x, self.di.y, self.di.z) self.data = np.reshape(np.fromfile(datafile, np.float32, -1), shape, 'F') if offset > 0: # TODO: Why is everything shifted by one unit self.data = np.roll(np.roll(np.roll(self.data,offset,axis=0),offset,axis=1),offset,axis=2) #make first offset row(s) and collum(s) zero self.data[0:,0:,:offset] = self.data[0:,:offset,0:] = self.data[:offset,0:,0:] = 0