diff --git a/atlasviewer/atlasviewer/dtissue.py b/atlasviewer/atlasviewer/dtissue.py new file mode 100755 index 0000000000000000000000000000000000000000..7e1993a9b79bfab757350dcbda524c13f1fdc2ed --- /dev/null +++ b/atlasviewer/atlasviewer/dtissue.py @@ -0,0 +1,140 @@ +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 + +