Source code for xmipp3.protocols.protocol_generate_reprojections

# **************************************************************************
# *
# * Authors:     Amaya Jimenez (
# *              Javier Mota (
# *
# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 2 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307  USA
# *
# *  All comments concerning this program package may be sent to the
# *  e-mail address ''
# *
# **************************************************************************

import os

from pyworkflow import VERSION_2_0
from pyworkflow.protocol.params import PointerParam
from pyworkflow.utils.path import cleanPath
from pwem.protocols import ProtAnalysis3D
from pwem.objects import Image
import pwem.emlib.metadata as md

from pwem import emlib
from xmipp3.convert import setXmippAttributes, xmippToLocation

[docs]class XmippProtGenerateReprojections(ProtAnalysis3D): """Compares a set of classes or averages with the corresponding projections of a reference volume. The set of images must have a 3D angular assignment and the protocol computes the residues (the difference between the experimental images and the reprojections). The zscore of the mean and variance of the residues are computed. Large values of these scores may indicate outliers. The protocol also analyze the covariance matrix of the residual and computes the logarithm of its determinant [Cherian2013]. The extremes of this score (called zScoreResCov), that is values particularly low or high, may indicate outliers.""" _label = 'generate reprojections' _lastUpdateVersion = VERSION_2_0 def __init__(self, **args): ProtAnalysis3D.__init__(self, **args) #--------------------------- DEFINE param functions -------------------------------------------- def _defineParams(self, form): form.addSection(label='Input') form.addParam('inputSet', PointerParam, label="Input images", important=True, pointerClass='SetOfParticles') form.addParam('inputVolume', PointerParam, label="Volume to compare images to", important=True, pointerClass='Volume', help='Volume to be used for class comparison') form.addParallelSection(threads=0, mpi=8) #--------------------------- INSERT steps functions -------------------------------------------- def _insertAllSteps(self): # Convert input images if necessary self.imgsFn = self._getExtraPath('input_imgs.xmd') vol = self.inputVolume.get() self._insertFunctionStep("convertStep", self.imgsFn) imgSet = self.inputSet.get() anglesFn=self.imgsFn self._insertFunctionStep("produceProjections", vol.getFileName(), anglesFn, vol.getSamplingRate()) self._insertFunctionStep("createOutputStep") #--------------------------- STEPS functions ---------------------------------------------------
[docs] def convertStep(self, imgsFn): from xmipp3.convert import writeSetOfParticles imgSet = self.inputSet.get() writeSetOfParticles(imgSet, self.imgsFn) from pwem.emlib.image import ImageHandler img = ImageHandler() fnVol = self._getTmpPath("volume.vol") img.convert(self.inputVolume.get(), fnVol) xdim=self.inputVolume.get().getDim()[0] imgXdim=imgSet.getDim()[0] if xdim!=imgXdim: self.runJob("xmipp_image_resize","-i %s --dim %d"%(fnVol,imgXdim),numberOfMpi=1)
[docs] def produceProjections(self, fnVol, fnAngles, Ts): fnVol = self._getTmpPath("volume.vol") anglesOutFn=self._getExtraPath("anglesCont.stk") projectionsOutFn=self._getExtraPath("projections.stk") args="-i %s -o %s --ref %s --oprojections %s --sampling %f " \ "--max_angular_change 90"%(fnAngles,anglesOutFn,fnVol, projectionsOutFn,Ts) self.runJob("xmipp_angular_continuous_assign2", args) fnNewParticles=self._getExtraPath("images.stk") if os.path.exists(fnNewParticles): cleanPath(fnNewParticles)
[docs] def createOutputStep(self): fnImgs = self._getExtraPath('images.stk') if os.path.exists(fnImgs): cleanPath(fnImgs) imgSet = self.inputSet.get() imgFn = self._getExtraPath("anglesCont.xmd") self.newAssignmentPerformed = os.path.exists(self._getExtraPath("angles.xmd")) self.samplingRate = self.inputSet.get().getSamplingRate() outputSet = self._createSetOfParticles() outputSet.copyInfo(imgSet) if not self.newAssignmentPerformed: outputSet.setAlignmentProj() self.iterMd = md.iterRows(imgFn, md.MDL_ITEM_ID) self.lastRow = next(self.iterMd) outputSet.copyItems(imgSet, updateItemCallback=self._updateItem ) self._defineOutputs(outputParticles=outputSet) self._defineSourceRelation(self.inputSet, outputSet) imgSet = self.inputSet.get() outputSet2 = self._createSetOfParticles('2') outputSet2.copyInfo(imgSet) if not self.newAssignmentPerformed: outputSet2.setAlignmentProj() self.iterMd2 = md.iterRows(imgFn, md.MDL_ITEM_ID) self.lastRow2 = next(self.iterMd2) outputSet2.copyItems(imgSet, updateItemCallback=self._updateItem2, ) self._defineOutputs(outputProjections=outputSet2) self._defineSourceRelation(self.inputSet, outputSet2)
def _updateItem(self, particle, row): count = 0 while self.lastRow and particle.getObjId() == self.lastRow.getValue( md.MDL_ITEM_ID): count += 1 if count: self._processRow(particle, self.lastRow) try: self.lastRow = next(self.iterMd) except StopIteration: self.lastRow = None particle._appendItem = count > 0 def _updateItem2(self, particle2, row): count = 0 while self.lastRow2 and particle2.getObjId() == self.lastRow2.getValue( md.MDL_ITEM_ID): count += 1 if count: self._processRow2(particle2, self.lastRow2) try: self.lastRow2 = next(self.iterMd2) except StopIteration: self.lastRow2 = None particle2._appendItem = count > 0 def _processRow(self, particle, row): def __setXmippImage(label): attr = '_xmipp_' + emlib.label2Str(label) if not hasattr(particle, attr): img = Image() setattr(particle, attr, img) img.setSamplingRate(particle.getSamplingRate()) else: img = getattr(particle, attr) img.setLocation(xmippToLocation(row.getValue(label))) particle.setLocation(xmippToLocation(row.getValue(emlib.MDL_IMAGE))) #__setXmippImage(emlib.MDL_IMAGE) #__setXmippImage(emlib.MDL_IMAGE_REF) setXmippAttributes(particle, row, emlib.MDL_IMAGE_ORIGINAL, emlib.MDL_IMAGE_REF) def _processRow2(self, particle, row): def __setXmippImage(label): attr = '_xmipp_' + emlib.label2Str(label) if not hasattr(particle, attr): img = Image() setattr(particle, attr, img) img.setSamplingRate(particle.getSamplingRate()) else: img = getattr(particle, attr) img.setLocation(xmippToLocation(row.getValue(label))) particle.setLocation(xmippToLocation(row.getValue( emlib.MDL_IMAGE_REF))) #__setXmippImage(emlib.MDL_IMAGE) #__setXmippImage(emlib.MDL_IMAGE_REF) setXmippAttributes(particle, row, emlib.MDL_IMAGE_ORIGINAL, emlib.MDL_IMAGE_REF) #--------------------------- INFO functions -------------------------------------------- def _summary(self): summary = [] summary.append("Images evaluated: %i" % self.inputSet.get().getSize()) summary.append("Volume: %s" % self.inputVolume.getNameId()) return summary def _methods(self): methods = [] if hasattr(self, 'outputParticles'): methods.append( "We projected the volume %s following the directions in %i input images %s." \ % ( self.getObjectTag('inputVolume'), self.inputSet.get().getSize(), self.getObjectTag('inputSet'))) return methods