# -*- coding: utf-8 -*-
# **************************************************************************
# *
# * Authors: Jose Luis Vilas (jlvilas@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'
# *
# **************************************************************************
import numpy as np
from pyworkflow import VERSION_2_0
from pyworkflow.protocol.params import (PointerParam, BooleanParam, FloatParam,
LEVEL_ADVANCED)
from pwem.protocols import ProtAnalysis3D
from pyworkflow.object import Float
from pwem.emlib.image import ImageHandler
from pyworkflow.utils import getExt
from pwem.objects import Volume
import pwem.emlib.metadata as md
MONORES_METHOD_URL = 'http://github.com/I2PC/scipion/wiki/XmippProtMonoRes'
OUTPUT_RESOLUTION_FILE = 'resolutionMap'
FN_FILTERED_MAP = 'filteredMap'
FN_MEAN_VOL = 'meanvol'
FN_MASK_WEDGE = 'maskwedge'
METADATA_MASK_FILE = 'metadataresolutions'
FN_METADATA_HISTOGRAM = 'mdhist'
BINARY_MASK = 'binarymask'
FN_GAUSSIAN_MAP = 'gaussianfilter'
[docs]class XmippProtMonoTomo(ProtAnalysis3D):
"""
Given a map the protocol assigns local resolutions to each voxel of the map.
"""
_label = 'local MonoTomo'
_lastUpdateVersion = VERSION_2_0
def __init__(self, **args):
ProtAnalysis3D.__init__(self, **args)
self.min_res_init = Float()
self.max_res_init = Float()
# --------------------------- DEFINE param functions ----------------------
def _defineParams(self, form):
form.addSection(label='Input')
form.addParam('inputVolume', PointerParam, pointerClass='Volume',
label="Odd tomogram", important=True,
help='Select a volume for determining its '
'local resolution.')
form.addParam('inputVolume2', PointerParam, pointerClass='Volume',
label="Even Tomogram", important=True,
help='Select a second volume for determining a '
'local resolution.')
form.addParam('useMask', BooleanParam, default=False,
label="Use mask?",
help='The mask determines which points are specimen'
' and which are not.')
form.addParam('Mask', PointerParam, pointerClass='VolumeMask',
condition='useMask', allowsNull=True,
label="Binary Mask",
help='The mask determines which points are specimen'
' and which are not')
group = form.addGroup('Extra parameters')
line = group.addLine('Resolution Range (Å)',
help="Resolution range (and step in expert mode) "
"to evaluate the local resolution.")
group.addParam('significance', FloatParam, default=0.95,
expertLevel=LEVEL_ADVANCED,
label="Significance",
help='Relution is computed using hipothesis tests, '
'this value determines the significance of that test')
line.addParam('minRes', FloatParam, default=0, label='High')
line.addParam('maxRes', FloatParam, allowsNull=True, label='Low')
line.addParam('stepSize', FloatParam, allowsNull=True, default=0.25,
expertLevel=LEVEL_ADVANCED, label='Step')
form.addParallelSection(threads=4, mpi=0)
# --------------------------- INSERT steps functions --------------------------------------------
def _createFilenameTemplates(self):
""" Centralize how files are called """
myDict = {
FN_MEAN_VOL: self._getExtraPath('mean_volume.mrc'),
FN_MASK_WEDGE: self._getExtraPath('wedgeMask.mrc'),
FN_FILTERED_MAP: self._getExtraPath('filteredMap.vol'),
OUTPUT_RESOLUTION_FILE: self._getExtraPath('mgresolution.mrc'),
METADATA_MASK_FILE: self._getExtraPath('mask_data.xmd'),
FN_METADATA_HISTOGRAM: self._getExtraPath('hist.xmd'),
BINARY_MASK: self._getExtraPath('binarized_mask.mrc'),
FN_GAUSSIAN_MAP: self._getExtraPath('gaussianfilted.mrc'),
}
self._updateFilenamesDict(myDict)
def _insertAllSteps(self):
# Convert input into xmipp Metadata format
self._createFilenameTemplates()
self._insertFunctionStep('convertInputStep')
self._insertFunctionStep('resolutionMonoTomoStep')
self._insertFunctionStep('createOutputStep')
self._insertFunctionStep("createHistrogram")
[docs] def ifNomask(self, fnVol):
xdim, _ydim, _zdim = self.inputVolume.get().getDim()
params = ' -i %s' % fnVol
params += ' -o %s' % self._getFileName(FN_GAUSSIAN_MAP)
setsize = 0.02 * xdim
params += ' --fourier real_gaussian %f' % (setsize)
self.runJob('xmipp_transform_filter', params)
img = ImageHandler().read(self._getFileName(FN_GAUSSIAN_MAP))
imgData = img.getData()
max_val = np.amax(imgData) * 0.05
params = ' -i %s' % self._getFileName(FN_GAUSSIAN_MAP)
params += ' --select below %f' % max_val
params += ' --substitute binarize'
params += ' -o %s' % self._getFileName(BINARY_MASK)
self.runJob('xmipp_transform_threshold', params)
self.maskFn = self._getFileName(BINARY_MASK)
[docs] def maskRadius(self):
xdim, _ydim, _zdim = self.inputVolume.get().getDim()
xdim = xdim * 0.5
return xdim
[docs] def resolutionMonoTomoStep(self):
# Number of frequencies
max_ = self.maxRes.get()
if self.stepSize.hasValue():
freq_step = self.stepSize.get()
else:
freq_step = 0.25
xdim = self.maskRadius()
params = ' --vol %s' % self.vol1Fn
params += ' --vol2 %s' % self.vol2Fn
params += ' --meanVol %s' % self._getFileName(FN_MEAN_VOL)
if self.useMask.get():
params += ' --mask %s' % self._getFileName(BINARY_MASK)
params += ' --sampling_rate %f' % self.inputVolume.get().getSamplingRate()
params += ' --minRes %f' % self.minRes.get()
params += ' --maxRes %f' % max_
params += ' --step %f' % freq_step
params += ' -o %s' % self._getFileName(OUTPUT_RESOLUTION_FILE)
params += ' --significance %f' % self.significance.get()
self.runJob('xmipp_resolution_monotomo', params)
[docs] def createHistrogram(self):
freq_step = self.stepSize.get() if self.stepSize.hasValue() else 10
M = float(self.max_res_init)
m = float(self.min_res_init)
range_res = round((M - m) / freq_step)
params = ' -i %s' % self._getFileName(OUTPUT_RESOLUTION_FILE)
if self.useMask.get():
params += ' --mask binary_file %s' % self._getFileName(BINARY_MASK)
params += ' --steps %f' % range_res
params += ' --range %f %f' % (self.min_res_init.get(),
(float(self.max_res_init.get()) - float(freq_step)))
params += ' -o %s' % self._getFileName(FN_METADATA_HISTOGRAM)
self.runJob('xmipp_image_histogram', params)
[docs] def getMinMax(self, imageFile):
img = ImageHandler().read(imageFile)
imgData = img.getData()
min_res = round(np.amin(imgData) * 100) / 100
max_res = round(np.amax(imgData) * 100) / 100
return min_res, max_res
[docs] def createOutputStep(self):
volume = Volume()
volume.setFileName(self._getFileName(OUTPUT_RESOLUTION_FILE))
volume.setSamplingRate(self.inputVolume.get().getSamplingRate())
self._defineOutputs(resolution_Volume=volume)
self._defineSourceRelation(self.inputVolume, volume)
# Setting the min max for the summary
imageFile = self._getFileName(OUTPUT_RESOLUTION_FILE)
min_, max_ = self.getMinMax(imageFile)
self.min_res_init.set(round(min_ * 100) / 100)
self.max_res_init.set(round(max_ * 100) / 100)
self._store(self.min_res_init)
self._store(self.max_res_init)
# --------------------------- INFO functions ------------------------------
def _methods(self):
messages = []
if hasattr(self, 'resolution_Volume'):
messages.append(
'Information about the method/article in ' + MONORES_METHOD_URL)
return messages
def _summary(self):
summary = []
summary.append("Highest resolution %.2f Å, "
"Lowest resolution %.2f Å. \n" % (self.min_res_init,
self.max_res_init))
return summary
def _citations(self):
return ['Vilas2018']