home/diploma/lws_root/utils/volumeimage.py

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