add dynamic tissue to AtlasViewer

atlasviewer/atlasviewer/dtissue.py

+import cPickle
+import numpy as np
+from atlasviewer.avtiff import voxelDimensionsFromTiffFile, avimread
+
+
+class DTissue(object):
+    def __init__(self, timePoints = None, background = 1):
+        self.timePoints = timePoints
+        self.background = 1
+        if self.timePoints is not None:
+            self.dtissue = {tp : {} for tp in self.timePoints}
+        else:
+            self.dtissue = {}
+            
+    def setSegmentationFiles(self, timePointSegmentationFiles):
+        for tp in self.dtissue:
+            self.dtissue[tp]["segmentation_file"] = timePointSegmentationFiles[tp]
+    
+    def setResolutions(self, resolutions = None):
+        if resolutions is None:
+            for tp in self.dtissue:
+                fn = self.dtissue[tp]["segmentation_file"]
+                res = voxelDimensionsFromTiffFile(fn)
+                self.dtissue[tp]["resolution"] = res
+        else:
+            for tp in self.dtissue:
+                self.dtissue[tp]["resolution"] = resolutions[tp]
+    
+    def extractLabels(self):
+        print "Extraing labels ... "
+        for tp in self.dtissue:
+            image, info = avimread(self.dtissue[tp]["segmentation_file"])
+            labels = list(np.unique(image))
+            self.dtissue[tp]["labels"] = labels
+            print tp, 
+        print "\n Done."
+        
+    def computeVolumes(self):
+        volumesDict = {}
+        print "Computing volumes ..."
+        for tp in self.dtissue:
+            print tp
+            image, tags = avimread(self.dtissue[tp]["segmentation_file"])
+            bc = np.bincount(image.flatten())
+            voxelS = np.prod(self.dtissue[tp]["resolution"])
+            vols = bc * voxelS
+            self.dtissue[tp]["volumes"] = {cid: vols[cid] for cid in self.dtissue[tp]["labels"] if cid != self.background}
+        print " Done."
+            
+    def setRealTimepoints(self, realTPD):
+        for tp in self.dtissue:
+            self.dtissue[tp]["real_timepoint"] = realTPD[tp]
+            
+    def getDaughters(self, motherTp, motherId):
+        try:
+            return self.dtissue[motherTp]["daughters"][motherId]
+        except AttributeError, e:
+            print "Daughters property is not set."
+            
+    def getMother(self, daughterTp, daughterId):
+        try:
+            return self.dtissue[daughterTp]["mother"][daughterId]
+        except AttributeError, e:
+            print "Daughters property is not set."       
+                    
+    def __getitem__(self, index): 
+        return self.dtissue[index]
+        
+    def extractDescendants(self, motherTp, motherCellId, descendantsTp):
+        family = {motherTp: [motherCellId]}
+        motherTpIndex, desesendantsTpIndex = self.timePoints.index(motherTp), self.timePoints.index(descendantsTp)
+        for ind in xrange(motherTpIndex, desesendantsTpIndex):
+            tp = self.timePoints[ind]
+            family[self.timePoints[ind + 1]] = []
+            for cid in family[tp]:
+                family[self.timePoints[ind + 1]].extend(self.dtissue[tp]["daughters"][cid])
+        return family[descendantsTp]    
+    
+
+    def setTrackingFiles(self, trackingFileNames):
+        for (tp0, tp1), fname in trackingFileNames.iteritems():
+            tp = tp0
+            self.dtissue[tp]["daughters"] = {}
+            fobj = file(fname)
+            tracking = cPickle.load(fobj)[0]
+            fobj.close()
+            for (daughter, mother) in tracking:
+                self.dtissue[tp]["daughters"].setdefault(mother, []).append(daughter)
+            for cid in  (set(self.dtissue[tp]["labels"]) - set(self.dtissue[tp]["daughters"].keys())):
+                self.dtissue[tp]["daughters"][cid] = []
+        
+        for (tp0, tp1), fname in trackingFileNames.iteritems():
+            tp = tp1
+            self.dtissue[tp]["mother"] = {}
+            fobj = file(fname)
+            tracking = cPickle.load(fobj)[0]
+            fobj.close()
+            for (daughter, mother) in tracking:
+                self.dtissue[tp]["mother"][daughter] =  mother
+            for cid in  (set(self.dtissue[tp]["labels"]) - set(self.dtissue[tp]["mother"].keys())):
+                self.dtissue[tp]["mother"][cid] = -1
+    
+        
+    def setPatternFiles(self, patternFNames):
+        for tp, fn in patternFNames.iteritems():
+            fobj = file(fn)
+            patterns = cPickle.load(fobj)
+            fobj.close()
+            self.dtissue[tp]["patterns"] = patterns
+            print tp, "L1: ", patterns["L1"]
+            
+    def computeGrowthRate(self, sourceTP, targetTP):
+        growthRate = {}
+        withDes = []
+        growth = []
+        deltaT = (self.dtissue[targetTP]["real_timepoint"] - self.dtissue[sourceTP]["real_timepoint"]) 
+        for cid in self.dtissue[sourceTP]["labels"]:
+            if cid != self.background:
+                descendants = self.extractDescendants(sourceTP, cid, targetTP)
+                if len(descendants) > 0:# and (cid in self.dtissue[sourceTP]["patterns"]["L1"]):
+                    dVolume = sum([self.dtissue[targetTP]["volumes"][did] for did in descendants])
+                    mVolume = self.dtissue[sourceTP]["volumes"][cid]
+                    growthRate[cid] = (dVolume - mVolume) / (deltaT * mVolume)
+        return growthRate
+                    
+    def computeRelativeGrowthRate(self, sourceTP, targetTP):
+        relativeGrowthRate = {}
+        withDes = []
+        growth = []
+        deltaT = (self.dtissue[targetTP]["real_timepoint"] - self.dtissue[sourceTP]["real_timepoint"]) 
+        for cid in self.dtissue[sourceTP]["labels"]:
+            if cid != self.background:
+                descendants = self.extractDescendants(sourceTP, cid, targetTP)
+                if len(descendants) > 0:# and (cid in self.dtissue[sourceTP]["patterns"]["L1"]):
+                    dVolume = sum([self.dtissue[targetTP]["volumes"][did] for did in descendants])
+                    mVolume = self.dtissue[sourceTP]["volumes"][cid]
+                    relativeGrowthRate[cid] = (dVolume - mVolume) / (deltaT * mVolume)
+        return relativeGrowthRate
+            
+