Source code for xmipp3.viewers.viewer_resolution_fso

# **************************************************************************
# *
# * Authors:     J.L. 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'
# *
# **************************************************************************

from pyworkflow.gui.plotter import Plotter
from pyworkflow.protocol.params import LabelParam, StringParam, EnumParam, IntParam, LEVEL_ADVANCED
from pyworkflow.viewer import ProtocolViewer, DESKTOP_TKINTER
from pyworkflow.utils import getExt, removeExt

from pwem.emlib.metadata import MetaData
from pwem.wizards import ColorScaleWizardBase
from pwem.emlib.image import ImageHandler
from pwem.viewers import ChimeraView, DataView, LocalResolutionViewer
from pwem import emlib


import numpy as np

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import FuncFormatter

from xmipp3.protocols.protocol_resolution_fso import \
    XmippProtFSO, OUTPUT_3DFSC, OUTPUT_DIRECTIONAL_FILTER
from xmipp3.viewers.plotter import XmippPlotter
import xmippLib

# Axis code
AX_X = 0
AX_Y = 1
AX_Z = 2


[docs]class XmippProtFSOViewer(LocalResolutionViewer): """ Visualization tools for the FSO, FSC, and 3DFSC. """ _label = 'viewer FSO' _targets = [XmippProtFSO] _environments = [DESKTOP_TKINTER]
[docs] @staticmethod def getColorMapChoices(): return plt.colormaps()
def _defineParams(self, form): form.addSection(label='Visualization') form.addParam('doShowOriginalVolumeSlices', LabelParam, label="Original Half Maps slices") form.addParam('doShowDirectionalFilter', LabelParam, label="Directionally filtered map") groupDirFSC = form.addGroup('FSO and Resolution Analysis') groupDirFSC.addParam('doShowAnisotropyCurve', LabelParam, label="Show FSO curve") groupDirFSC.addParam('doShowFSC', LabelParam, label="Global FSC") groupDirFSC.addParam('doShow3DFSC', LabelParam, label="3DFSC map") groupDirFSC.addParam('doShow3DFSCcolorSlices', LabelParam, label="Show 3DFSC Color slices") groupDirFSC.addParam('doShowDirectionalResolution', LabelParam, label="Show Directional Resolution on sphere") group = form.addGroup('Choose a Color Map') group.addParam('sliceAxis', EnumParam, default=AX_Z, choices=['x', 'y', 'z'], display=EnumParam.DISPLAY_HLIST, label='Slice axis') group.addParam('doShowOneColorslice', LabelParam, expertLevel=LEVEL_ADVANCED, label='Show selected slice') group.addParam('sliceNumber', IntParam, default=-1, expertLevel=LEVEL_ADVANCED, label='Show slice number') ColorScaleWizardBase.defineColorScaleParams(group, defaultLowest=0.0, defaultHighest=1.0) def _getVisualizeDict(self): self.protocol._createFilenameTemplates() return {'doShowOriginalVolumeSlices': self._showOriginalVolumeSlices, 'doShow3DFSC': self._show3DFSC, 'doShowFSC': self._showFSCCurve, 'doShow3DFSCcolorSlices': self._show3DFSCcolorSlices, 'doShowOneColorslice': self._showOneColorslice, 'doShowAnisotropyCurve': self._showAnisotropyCurve, 'doShowDirectionalFilter': self._showDirectionalFilter, 'doShowDirectionalResolution': self._showDirectionalResolution, } def _showOriginalVolumeSlices(self, param=None): """ This function opens the two half maps to visualize the slices """ cm = DataView(self.protocol.half1.get().getFileName()) cm2 = DataView(self.protocol.half2.get().getFileName()) return [cm, cm2] def _show3DFSC(self, param=None): """ This function opens the slices of the 3DFSC """ cm = DataView(self.protocol._getExtraPath(OUTPUT_3DFSC)) return [cm] def _showFSCCurve(self, paramName=None): """ It shows the FSC curve in terms of the resolution The horizontal axis is linear in this plot. Note That this is not the normal representation of hte FSC """ fnmd = self.protocol._getExtraPath('GlobalFSC.xmd') title = 'Global FSC' xTitle = 'Resolution (1/A)' yTitle = 'FSC (a.u.)' mdLabelX = emlib.MDL_RESOLUTION_FREQ mdLabelY = emlib.MDL_RESOLUTION_FRC self._plotCurveFSC(fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY) def _show3DFSCcolorSlices(self, param=None): """ It opens 4 colores slices of the 3DFSC """ img = ImageHandler().read(self.protocol._getExtraPath(OUTPUT_3DFSC)) imgData = img.getData() xplotter = XmippPlotter(x=2, y=2, mainTitle="3DFSC Color Slices" "along %s-axis." % self._getAxis()) # The slices to be shown are close to the center. Volume size is divided in # 9 segments, the fouth central ones are selected i.e. 3,4,5,6 for i in range(3, 7): sliceNumber = self.getSlice(i, imgData) a = xplotter.createSubPlot("Slice %s" % (sliceNumber + 1), '', '') matrix = self.getSliceImage(imgData, sliceNumber, self._getAxis()) plot = xplotter.plotMatrix(a, matrix, 0, 1, cmap=self.getColorMap(), interpolation="nearest") xplotter.getColorBar(plot) return [xplotter] def _showAnisotropyCurve(self, paramName=None): """ It shows the FSO curve in terms of the resolution The horizontal axis is linear in this plot. """ fnmd = self.protocol._getExtraPath('fso.xmd') title = 'Anisotropy Curve' xTitle = 'Resolution (1/A)' yTitle = 'Anisotropy (a.u.)' mdLabelX = emlib.MDL_RESOLUTION_FREQ mdLabelY = emlib.MDL_RESOLUTION_FSO self._plotCurveAnisotropy(fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY) def _showDirectionalFilter(self, param=None): """ The directionally filtered map using the 3DFSC as low pass filter """ cm = DataView(self.protocol._getExtraPath(OUTPUT_DIRECTIONAL_FILTER)) return [cm] def _showOneColorslice(self, param=None): """ It shows a single colored slice of the 3DFSC """ imageFile = self.protocol._getExtraPath(OUTPUT_3DFSC) img = ImageHandler().read(imageFile) imgData = img.getData() xplotter = XmippPlotter(x=1, y=1, mainTitle="3DFSC Slices " "along %s-axis." % self._getAxis()) sliceNumber = self.sliceNumber.get() if sliceNumber < 0: x, _, _, _ = ImageHandler().getDimensions(imageFile) sliceNumber = int(x / 2) else: sliceNumber -= 1 # sliceNumber has no sense to start in zero a = xplotter.createSubPlot("Slice %s" % (sliceNumber + 1), '', '') matrix = self.getSliceImage(imgData, sliceNumber, self._getAxis()) plot = xplotter.plotMatrix(a, matrix, 0, 1, cmap=self.getColorMap(), interpolation="nearest") xplotter.getColorBar(plot) return [plt.show(xplotter)] def _formatFreq(self, value, pos): """ Format function for Matplotlib formatter. """ inv = 999. if value: inv = 1 / value return "1/%0.2f" % inv def _plotCurveFSC(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY): """ This function is called by _showFSCCurve. It shows the FSC curve in terms of the resolution That this is not the normal representation of the FSC """ md = xmippLib.MetaData(fnmd) xplotter = XmippPlotter(figure=None) xplotter.plot_title_fontsize = 11 a = xplotter.createSubPlot(title, xTitle, yTitle, 1, 1) xplotter.plotMdFile(md, mdLabelX, mdLabelY, 'g') a.xaxis.set_major_formatter(FuncFormatter(self._formatFreq)) xx, yy = self._prepareDataForPlot(md, mdLabelX, mdLabelY) a.hlines(0.143, xx[0], xx[-1], colors='k', linestyles='dashed') a.grid(True) return plt.show()
[docs] def interpolRes(self, thr, x, y): """ This function is called by _showAnisotropyCurve. It provides the cut point of the curve defined by the points (x,y) with a threshold thr. The flag okToPlot shows if there is no intersection points """ idx = np.arange(0, len(x)) aux = np.array(y) <= thr idx_x = idx[aux] okToPlot = True resInterp = [] if not idx_x.any(): okToPlot = False else: if len(idx_x) > 1: idx_2 = idx_x[0] idx_1 = idx_2 - 1 if idx_1 < 0: idx_2 = idx_x[1] idx_1 = idx_2 - 1 y2 = x[idx_2] y1 = x[idx_1] x2 = y[idx_2] x1 = y[idx_1] slope = (y2 - y1) / (x2 - x1) ny = y2 - slope * x2 resInterp = 1.0 / (slope * thr + ny) else: okToPlot = False return resInterp, okToPlot
def _plotCurveAnisotropy(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY): """ This function is called by _showAnisotropyCurve It shows the FSO curve in terms of the resolution The horizontal axis is linear in this plot. """ md = xmippLib.MetaData(fnmd) xplotter = XmippPlotter(figure=None) xplotter.plot_title_fontsize = 11 a = xplotter.createSubPlot(title, xTitle, yTitle, 1, 1) xplotter.plotMdFile(md, mdLabelX, mdLabelY, 'g') xx, yy = self._prepareDataForPlot(md, mdLabelX, mdLabelY) from matplotlib.ticker import FuncFormatter a.axes.xaxis.set_major_formatter(FuncFormatter(self._formatFreq)) a.axes.set_ylim([-0.1, 1.1]) a.axes.plot(xx, yy, 'g') a.axes.set_xlabel('Resolution (A)') a.axes.set_ylabel('FSO (a.u)') hthresholds = [0.1, 0.5, 0.9] a.axes.hlines(hthresholds, xx[0], xx[-1], colors='k', linestyles='dashed') a.axes.grid(True) textstr = '' res_01, okToPlot_01 = self.interpolRes(0.1, xx, yy) res_05, okToPlot_05 = self.interpolRes(0.5, xx, yy) res_09, okToPlot_09 = self.interpolRes(0.9, xx, yy) if (okToPlot_01 and okToPlot_05 and okToPlot_09): textstr = str(0.9) + ' --> ' + str("{:.2f}".format(res_09)) + 'A\n' + str(0.5) + ' --> ' + str( "{:.2f}".format(res_05)) + 'A\n' + str(0.1) + ' --> ' + str("{:.2f}".format(res_01)) + 'A' a.axes.axvspan(1.0 / res_09, 1.0 / res_01, alpha=0.3, color='green') props = dict(boxstyle='round', facecolor='white') a.axes.text(0.0, 0.0, textstr, fontsize=12, ha="left", va="bottom", bbox=props) return plt.show() def _prepareDataForPlot(self, md, mdLabelX, mdLabelY): """ plot metadata columns mdLabelX and mdLabelY if nbins is in args then and histogram over y data is made """ xx = md.getColumnValues(mdLabelX) yy = md.getColumnValues(mdLabelY) return xx, yy def _showDirectionalResolution(self, zparam=None): """ This function shows the angular distribution of the resolution """ fnmd = self.protocol._getExtraPath('Resolution_Distribution.xmd') self._showPolarPlot(fnmd) def _showPolarPlot(self, fnmd): """ It is called by _showDirectionalResolution This function shows the angular distribution of the resolution """ md = emlib.MetaData(fnmd) radius = md.getColumnValues(emlib.MDL_ANGLE_ROT) azimuth = md.getColumnValues(emlib.MDL_ANGLE_TILT) counts = md.getColumnValues(emlib.MDL_RESOLUTION_FRC) # define binning azimuths = np.radians(np.linspace(0, 360, 360)) zeniths = np.arange(0, 91, 1) r, theta = np.meshgrid(zeniths, azimuths) values = np.zeros((len(azimuths), len(zeniths))) for i in range(0, len(azimuth)): values[int(radius[i]), int(azimuth[i])] = counts[i] # ------ Plot ------ stp = 0.1 lowlim = max(0.0, values.min()) highlim = values.max() + stp fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) pc = plt.contourf(theta, r, values, np.arange(lowlim, highlim, stp)) plt.colorbar(pc) plt.show() def _getAxis(self): return self.getEnumText('sliceAxis')
[docs] def getColorMap(self): cmap = cm.get_cmap(self.colorMap.get()) if cmap is None: cmap = cm.jet return cmap