# **************************************************************************
# *
# * Authors: Javier Vargas (jvargas@cnb.csic.es)
# * Adrian Quintana (aquintana@cnb.csic.es)
# *
# * 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
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * 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 'scipion@cnb.csic.es'
# *
# **************************************************************************
from os.path import exists
from pyworkflow.object import String
from pyworkflow.protocol.constants import LEVEL_ADVANCED
import pyworkflow.protocol.params as params
from pyworkflow.utils.path import cleanPath
from pwem.protocols import ProtInitialVolume
from pwem.objects import Volume, SetOfParticles
from pwem.constants import ALIGN_2D
from pwem import emlib
from xmipp3.convert import getImageLocation, alignmentToRow, convertToMrc
[docs]class XmippProtRCT(ProtInitialVolume):
"""Creates initial volumes by using a set of projections/classes
from a tilted-pair picking process and using RCT algorithm. """
_label = 'random conical tilt'
def __init__(self, **args):
ProtInitialVolume.__init__(self, **args)
#--------------------------- DEFINE param functions ------------------------
def _defineParams(self, form):
form.addSection(label='Input')
# TODO: Input can be either a SetOfParticles or
# a SetOfClasses2D
form.addParam('inputParticlesTiltPair', params.PointerParam,
label="Input particles tilt pair",
important=True,
pointerClass='ParticlesTiltPair',
help='Select the input particles tilt pair file '
'that will be used. file. This file is used '
'to associate each micrograph with its tilted '
'equivalent.')
form.addParam('inputParticles', params.PointerParam,
label="Input classes", important=True,
pointerClass='SetOfParticles,SetOfClasses',
help='Select the input images or classes from '
'the project.')
form.addSection(label='Alignment parameters')
form.addParam('thinObject', params.BooleanParam,
default=False, label='Thin Object',
help='If the object is thin, then the tilted '
'projections can be stretched to match the '
'untilted projections')
form.addParam('maxShift', params.IntParam,
default="10", expertLevel=LEVEL_ADVANCED,
label="Maximum allowed shift for tilted particles (pixels)",
help='Particles that shift more will be discarded. '
'A value larger than the image size will not '
'discard any particle.')
form.addParam('skipTranslation', params.BooleanParam,
default=False, expertLevel=LEVEL_ADVANCED,
label='Skip tilted translation alignment',
help='If the tilted image quality is very low, '
'then this alignment might result in poor '
'estimates.')
form.addSection(label='Reconstruction')
form.addParam('additionalParams', params.StringParam,
default="-n 5 -l 0.01", expertLevel=LEVEL_ADVANCED,
label='Additional reconstruction parameters',
help='See: http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/Reconstruct_art_v31')
form.addParam('doFilter', params.BooleanParam, default=True,
label='Filter reconstructed volumes?',
help='Filtering may be useful to remove noise, '
'especially when few particles '
'contribute to the reconstruction.')
form.addParam('resoLowPassFilter', params.FloatParam, default=0.2,
label='Resolution of the low-pass filter (dig.freq)',
help='Resolution of the low-pass filter (dig.freq)')
form.addParallelSection(threads=4, mpi=1)
#--------------------------- STEPS functions -------------------------------
def _insertAllSteps(self):
self.inputSet = self.inputParticles.get()
self.rctClassesFn = self._getExtraPath('rct_classes.xmd')
firstStepId = self._insertFunctionStep('createRctImagesStep')
reconSteps = []
if isinstance(self.inputSet, SetOfParticles):
reconStep = self._reconstructImages(self.inputSet, firstStepId)
reconSteps.append(reconStep)
else:
for class2D in self.inputSet:
reconStep = self._reconstructImages(class2D, firstStepId)
reconSteps.append(reconStep)
self._insertFunctionStep('createOutputStep', prerequisites=reconSteps)
[docs] def createRctImagesStep(self):
""" Function to create the Xmipp metadata needed to run the
Xmipp protocol """
if isinstance(self.inputSet, SetOfParticles):
self._appendRctImages(self.inputSet)
else:
for class2D in self.inputSet:
self._appendRctImages(class2D)
def _appendRctImages(self, particles):
blockMd = "class%06d_images@%s" % (particles.getObjId(),
self.rctClassesFn)
classMd = emlib.MetaData()
partPairs = self.inputParticlesTiltPair.get()
uImages = partPairs.getUntilted()
tImages = partPairs.getTilted()
sangles = partPairs.getCoordsPair().getAngles()
micPairs = partPairs.getCoordsPair().getMicsPair()
uMics = micPairs.getUntilted()
tMics = micPairs.getTilted()
scaleFactor = uImages.getSamplingRate() / particles.getSamplingRate()
for img in particles:
imgId = img.getObjId()
uImg = uImages[imgId]
tImg = tImages[imgId]
if uImg is None or tImg is None:
print(">>> Warning, for id %d, tilted or untilted particle "
"was not found. Ignored." % imgId)
else:
objId = classMd.addObject()
pairRow = emlib.metadata.Row()
pairRow.setValue(emlib.MDL_IMAGE, getImageLocation(uImg))
uCoord = uImg.getCoordinate()
micId = uCoord.getMicId()
uMic = uMics[micId]
angles = sangles[micId]
pairRow.setValue(emlib.MDL_MICROGRAPH, uMic.getFileName())
pairRow.setValue(emlib.MDL_XCOOR, uCoord.getX())
pairRow.setValue(emlib.MDL_YCOOR, uCoord.getY())
pairRow.setValue(emlib.MDL_ENABLED, 1)
pairRow.setValue(emlib.MDL_ITEM_ID, int(imgId))
pairRow.setValue(emlib.MDL_REF, 1)
alignment = img.getTransform()
# Scale alignment by scaleFactor
alignment.scale(scaleFactor)
alignmentToRow(alignment, pairRow, alignType=ALIGN_2D)
pairRow.setValue(emlib.MDL_IMAGE_TILTED, getImageLocation(tImg))
tMic = tMics[micId]
pairRow.setValue(emlib.MDL_MICROGRAPH_TILTED, tMic.getFileName())
(angleY, angleY2, angleTilt) = angles.getAngles()
pairRow.setValue(emlib.MDL_ANGLE_Y, float(angleY))
pairRow.setValue(emlib.MDL_ANGLE_Y2, float(angleY2))
pairRow.setValue(emlib.MDL_ANGLE_TILT, float(angleTilt))
pairRow.writeToMd(classMd, objId)
classMd.write(blockMd, emlib.MD_APPEND)
def _reconstructImages(self, particles, deps):
""" Function to insert the step needed to reconstruct a
class (or setOfParticles) """
classNo = particles.getObjId()
blockMd = "class%06d_images@%s" % (classNo, self.rctClassesFn)
classNameIn = blockMd
classNameOut = self._getExtraPath("rct_images_%06d.xmd" % classNo)
classVolumeOut = self._getPath("rct_%06d.vol" % classNo)
if particles.hasRepresentative():
classImage = getImageLocation(particles.getRepresentative())
else:
classImage = None
reconStep = self._insertFunctionStep('reconstructClass',
classNameIn, classNameOut,
classImage, classVolumeOut,
prerequisites=[deps])
return reconStep
[docs] def reconstructClass(self, classIn, classOut, classImage, classVolumeOut):
# If class image doesn't exists, generate it by averaging
if classImage is None:
classRootOut = classOut.replace(".xmd", "") + "_"
statsFn = self._getExtraPath('stats.xmd')
args = "-i %(classIn)s --save_image_stats %(classRootOut)s -o %(statsFn)s"
self.runJob("xmipp_image_statistics", args % locals(),
numberOfMpi=1, numberOfThreads=1)
classImage = classRootOut + "average.xmp"
centerMaxShift = self.maxShift.get()
args = "-i %(classIn)s -o %(classOut)s --ref %(classImage)s " % locals()
args += " --max_shift %(centerMaxShift)d" % locals()
if self.thinObject.get():
args += " --do_stretch"
if self.skipTranslation.get():
args += " --do_not_align_tilted"
self.runJob("xmipp_image_align_tilt_pairs", args,
numberOfMpi=1, numberOfThreads=1)
reconstructAdditionalParams = self.additionalParams.get()
args = "-i %(classOut)s -o %(classVolumeOut)s " % locals()
args += " %(reconstructAdditionalParams)s" % locals()
if self.numberOfMpi >1:
# threading is not supported in mpi version
self.runJob("xmipp_reconstruct_art", args, numberOfThreads=1)
else:
args += " --thr %d" % self.numberOfThreads.get()
self.runJob("xmipp_reconstruct_art", args)
if exists(classVolumeOut):
mdFn = self._getPath('volumes.xmd')
md = emlib.MetaData()
if exists(mdFn):
md.read(mdFn)
objId = md.addObject()
md.setValue(emlib.MDL_IMAGE, classVolumeOut, objId)
if self.doFilter.get():
filteredVolume = classVolumeOut.replace('.vol', '_filtered.vol')
lowPassFilter = self.resoLowPassFilter.get()
args = "-i %(classVolumeOut)s -o %(filteredVolume)s " % locals()
args += " --fourier low_pass %(lowPassFilter)f " % locals()
args += " --thr %d" % self.numberOfThreads.get()
self.runJob("xmipp_transform_filter", args, numberOfMpi=1)
objId = md.addObject()
md.setValue(emlib.MDL_IMAGE, filteredVolume, objId)
md.write(mdFn)
[docs] def createOutputStep(self):
# TODO: Refactor following code if possible
self.volumesSet = self._createSetOfVolumes()
self.volumesSet.setStore(False)
self.sampling = self.inputParticlesTiltPair.get().getUntilted().getSamplingRate()
self.volumesSet.setSamplingRate(self.sampling)
if self.doFilter.get():
self.volumesFilterSet = self._createSetOfVolumes('filtered')
self.volumesFilterSet.setStore(False)
self.volumesFilterSet.setSamplingRate(self.sampling)
if isinstance(self.inputSet, SetOfParticles):
volumeOut = self._getPath("rct_%06d.vol" % self.inputSet.getObjId())
self._appendOutputVolume(volumeOut)
else:
for class2D in self.inputSet:
volumeOut = self._getPath("rct_%06d.vol" % class2D.getObjId())
self._appendOutputVolume(volumeOut)
self._defineOutputs(outputVolumes=self.volumesSet)
if self.doFilter.get():
self._defineOutputs(outputFilteredVolumes=self.volumesFilterSet)
self._defineSourceRelation(self.inputParticlesTiltPair, self.volumesSet)
def _appendOutputVolume(self, volumeOut):
fnMrc = convertToMrc(self, volumeOut, self.sampling, True)
vol = Volume()
vol.setFileName(fnMrc)
vol.setSamplingRate(self.sampling)
self.volumesSet.append(vol)
if self.doFilter.get():
volumeFilterOut = volumeOut.replace('.vol', '_filtered.vol')
fnMrc = convertToMrc(self, volumeFilterOut, self.sampling, True)
volf = Volume()
volf.setFileName(fnMrc)
volf.setSamplingRate(self.sampling)
self.volumesFilterSet.append(volf)
#--------------------------- INFO functions --------------------------------
def _validate(self):
errors = []
if isinstance(self.inputParticles.get(), SetOfParticles):
if not self.inputParticles.get().hasAlignment():
errors.append("The input particles should have "
"alignment information.")
else:
for class2D in self.inputParticles.get():
if not class2D.hasAlignment():
errors.append("The input classes should have "
"alignment information.")
return errors
def _summary(self):
summary = []
if not hasattr(self, 'outputVolumes'):
summary.append("Output volumes not ready yet.")
else:
if isinstance(self.inputParticles.get(), SetOfParticles):
summary.append("Input particles: %d"
% self.inputParticles.get().getSize())
else:
summary.append("Input classes: %d"
% self.inputParticles.get().getSize())
summary.append("Output volumes: %d"
% self.outputVolumes.getSize())
if self.doFilter.get():
summary.append("Output filtered volumes: %d"
% self.outputFilteredVolumes.getSize())
return summary
def _methods(self):
methods = []
if not hasattr(self, 'outputVolumes'):
methods.append("Output volumes not ready yet.")
else:
if isinstance(self.inputParticles.get(), SetOfParticles):
methods.append('Set of %d particles %s was employed to create '
'an initial volume using RCT method.'
% (len(self.inputParticles.get()),
self.getObjectTag('inputParticles')))
else:
particlesArray = [len(s) for s in self.inputParticles.get()]
particlesArrayString = String(particlesArray)
methods.append('Set of %d classes %s was employed to create %d '
'initial volumes using RCT method. '
% (len(self.inputParticles.get()),
self.getObjectTag('inputParticles'),
len(self.inputParticles.get())))
methods.append('For each initial volume were used respectively '
'%s particles' % particlesArrayString)
methods.append("Output volumes: %s" %
self.getObjectTag('outputVolumes'))
if self.doFilter.get():
methods.append("Output filtered volumes: %s"
% self.getObjectTag('outputFilteredVolumes'))
return methods
def _citations(self):
return ['Sorzano2015b']