Merge pull request #205 from scipion-em/devel

3.0.0b2

README.rst

@@ -38,16 +38,16 @@ b) Developer's version
 
       scipion installp -p path_to_scipion-em-relion --devel
 
-RELION sources will be downloaded and compiled automatically with the plugin, but you can also link an existing installation. Default installation path assumed is ``software/em/relion-3.1``, if you want to change it, set *RELION_HOME* in ``scipion.conf`` file to the folder where the RELION is installed. If you need to use CUDA different from the one used during Scipion installation (defined by *CUDA_LIB*), you can add *RELION_CUDA_LIB* variable to the config file. Moreover, if you have to use a MPI for Relion different from Scipion MPI, you can set *RELION_MPI_BIN* and *RELION_MPI_LIB* variables in your shell environment - they will be recognized by Scipion.
+RELION sources will be downloaded and compiled automatically with the plugin, but you can also link an existing installation. Default installation path assumed is ``software/em/relion-3.1.0``, if you want to change it, set *RELION_HOME* in ``scipion.conf`` file to the folder where the RELION is installed. If you need to use CUDA different from the one used during Scipion installation (defined by *CUDA_LIB*), you can add *RELION_CUDA_LIB* variable to the config file. Moreover, if you have to use a MPI for Relion different from Scipion MPI, you can set *RELION_MPI_BIN* and *RELION_MPI_LIB* variables in the config file.
 
 To check the installation, simply run one of the tests. A complete list of tests can be displayed by executing ``scipion test --show --grep relion``
 
 Supported versions
 ------------------
 
-3.0, 3.1
+3.0, 3.1.0
 
-In 2020 the plugin was updated to support the latest RELION 3.1. We discontinued any support for 2.x versions. We are still working towards stable plugin release, so some things might not work yet.
+In 2020 the plugin was updated to support the latest RELION 3.1.0. We discontinued any support for 2.x versions. We are still working towards stable plugin release, so some things might not work yet. Handling optics groups needs more testing, you are welcome to be a volunteer.
 
 Protocols
 ---------

relion/convert/__init__.py

@@ -165,7 +165,7 @@ def _updateClass(self, item):
             if self._alignType == pwem.ALIGN_PROJ:
                 fn += ':mrc'  # mark reference as a MRC volume
             item.getRepresentative().setLocation(index, fn)
-            item._rlnclassDistribution = Float(row.rlnClassDistribution)
+            item._rlnClassDistribution = Float(row.rlnClassDistribution)
             item._rlnAccuracyRotations = Float(row.rlnAccuracyRotations)
             if Plugin.IS_GT30():
                 item._rlnAccuracyTranslationsAngst = Float(row.rlnAccuracyTranslationsAngst)

relion/convert/convert30.py

@@ -26,6 +26,11 @@
 # *
 # **************************************************************************
 
+import numpy as np
+
+import pwem
+import pwem.convert.transformations as tfs
+
 from ..constants import *
 from .convert_base import WriterBase, ReaderBase
 from .convert_deprecated import rowToAlignment, readSetOfParticles
@@ -93,6 +98,15 @@ def _micToRow(self, mic, row):
 
 class Reader(ReaderBase):
 
+    ALIGNMENT_LABELS = [
+        "rlnOriginX",
+        "rlnOriginY",
+        "rlnOriginZ",
+        "rlnAngleRot",
+        "rlnAngleTilt",
+        "rlnAnglePsi",
+    ]
+
     def readSetOfParticles(self, starFile, partsSet, **kwargs):
         """ Convert a star file into a set of particles.
 
@@ -110,4 +124,68 @@ def readSetOfParticles(self, starFile, partsSet, **kwargs):
 
     def setParticleTransform(self, particle, row):
         """ Set the transform values from the row. """
-        particle.setTransform(rowToAlignment(row, self._alignType))
+
+        if ((self._alignType == pwem.ALIGN_NONE) or
+                not row.hasAnyColumn(self.ALIGNMENT_LABELS)):
+            self.setParticleTransform = self.__setParticleTransformNone
+        else:
+            # Ensure the Transform object exists
+            self._angles = np.zeros(3)
+            self._shifts = np.zeros(3)
+
+            particle.setTransform(pwem.objects.Transform())
+
+            if self._alignType == pwem.ALIGN_2D:
+                self.setParticleTransform = self.__setParticleTransform2D
+            elif self._alignType == pwem.ALIGN_PROJ:
+                self.setParticleTransform = self.__setParticleTransformProj
+            else:
+                raise Exception("Unexpected alignment type: %s"
+                                % self._alignType)
+
+        # Call again the modified function
+        self.setParticleTransform(particle, row)
+
+    def __setParticleTransformNone(self, particle, row):
+        particle.setTransform(None)
+
+    def __setParticleTransform2D(self, particle, row):
+        angles = self._angles
+        shifts = self._shifts
+        ips = self._invPixelSize
+
+        def _get(label):
+            return float(getattr(row, label, 0.))
+
+        shifts[0] = _get('rlnOriginX') * ips
+        shifts[1] = _get('rlnOriginY') * ips
+        angles[2] = -_get('rlnAnglePsi')
+        radAngles = -np.deg2rad(angles)
+        M = tfs.euler_matrix(radAngles[0], radAngles[1], radAngles[2], 'szyz')
+        M[:3, 3] = shifts[:3]
+        particle.getTransform().setMatrix(M)
+
+    def __setParticleTransformProj(self, particle, row):
+        angles = self._angles
+        shifts = self._shifts
+        ips = self._invPixelSize
+
+        def _get(label):
+            return float(getattr(row, label, 0.))
+
+        shifts[0] = _get('rlnOriginX') * ips
+        shifts[1] = _get('rlnOriginY') * ips
+        shifts[2] = _get('rlnOriginZ') * ips
+
+        angles[0] = _get('rlnAngleRot')
+        angles[1] = _get('rlnAngleTilt')
+        angles[2] = _get('rlnAnglePsi')
+
+        radAngles = -np.deg2rad(angles)
+
+        # TODO: jmrt: Maybe we should test performance and consider if keeping
+        # TODO: the matrix and not creating one everytime will make things faster
+        M = tfs.euler_matrix(radAngles[0], radAngles[1], radAngles[2], 'szyz')
+        M[:3, 3] = -shifts[:3]
+        M = np.linalg.inv(M)
+        particle.getTransform().setMatrix(M)

relion/convert/dataimport.py

@@ -135,7 +135,10 @@ def _updateClass(self, item):
             item.getRepresentative().setLocation(index, fn)
             item._rlnclassDistribution = Float(row.get('rlnClassDistribution'))
             item._rlnAccuracyRotations = Float(row.get('rlnAccuracyRotations'))
-            item._rlnAccuracyTranslations = Float(row.get('rlnAccuracyTranslations'))
+            if self.version30:
+                item._rlnAccuracyTranslations = Float(row.get('rlnAccuracyTranslations'))
+            else:
+                item._rlnAccuracyTranslations = Float(row.get('rlnAccuracyTranslationsAngst'))
 
     def _createClasses(self, partSet):
         self._classesDict = {}  # store classes info, indexed by class id

relion/protocols/_legacy/protocol30_ctf_refinement.py

@@ -42,8 +42,8 @@ class ProtRelionCtfRefinement(ProtParticles):
 
     def _defineParams(self, form):
         form.addSection(label='Input')
-        # TODO: conditions on particles?
         form.addParam('inputParticles', params.PointerParam,
+                      pointerCondition='hasAlignmentProj',
                       important=True,
                       label='Input particles',
                       pointerClass='SetOfParticles',
@@ -118,7 +118,7 @@ def _defineParams(self, form):
                            'and good signal-to-noise ratios per micrograph.')
         form.addParam('doBeamtiltEstimation', params.BooleanParam,
                       default=False,
-                      label='Perform beamtilt estimation?',
+                      label='Perform beam tilt estimation?',
                       help='If set to Yes, then relion_ctf_refine will also '
                            'estimate the beamtilt over the entire data set. '
                            'This option is only recommended for '
@@ -213,7 +213,6 @@ def createGlobalInfoStep(self):
 
     def _updateItemCtfBeamTilt(self, particle, row):
         particle.setCTF(convert.rowToCtfModel(row))
-        # TODO: Add other field from the .star file when other options?
         # check if beamtilt is available and save it
         if row.hasLabel('rlnBeamTiltX'):
             particle._rlnBeamTiltX = Float(row.getValue('rlnBeamTiltX', 0))

relion/protocols/protocol_assign_optic_groups.py

@@ -220,6 +220,9 @@ def updateItem(item, row):
     def _validate(self):
         validateMsgs = []
 
+        if self.mtfFile.hasValue() and not pwutils.exists(self.mtfFile.get()):
+            validateMsgs.append("MTF file %s does not exist!" % self.mtfFile.get())
+
         return validateMsgs
 
     def _summary(self):

relion/protocols/protocol_ctf_refinement.py

@@ -42,6 +42,33 @@ class ProtRelionCtfRefinement(ProtParticles):
     """ Wrapper protocol for the Relion's CTF refinement. """
     _label = 'ctf refinement'
 
+    def _initialize(self):
+        self._createFilenameTemplates()
+
+    def _createFilenameTemplates(self):
+        """ Centralize how files are called. """
+        myDict = {
+            'output_star': self._getExtraPath("particles_ctf_refine.star"),
+            'ctf_sqlite': self._getExtraPath("ctf_analyze.sqlite"),
+            'mag_obs_x': self._getExtraPath("mag_disp_x_optics-group_%(og)d.mrc:mrc"),
+            'mag_obs_y': self._getExtraPath("mag_disp_y_optics-group_%(og)d.mrc:mrc"),
+            'mag_fit_x': self._getExtraPath("mag_disp_x_fit_optics-group_%(og)d.mrc:mrc"),
+            'mag_fit_y': self._getExtraPath("mag_disp_y_fit_optics-group_%(og)d.mrc:mrc"),
+            'tetrafoil_it_fit': self._getExtraPath("aberr_delta-phase_iter-fit_optics-group_%(og)d_N-4.mrc:mrc"),
+            'tetrafoil_fit': self._getExtraPath("aberr_delta-phase_lin-fit_optics-group_%(og)d_N-4.mrc:mrc"),
+            'tetrafoil_residual_fit': self._getExtraPath("aberr_delta-phase_lin-fit_optics-group_%(og)d_N-4_residual.mrc:mrc"),
+            'tetrafoil_obs': self._getExtraPath("aberr_delta-phase_per-pixel_optics-group_%(og)d.mrc:mrc"),
+            'beamtilt_it_fit': self._getExtraPath("beamtilt_delta-phase_iter-fit_optics-group_%(og)d.mrc:mrc"),
+            'beamtilt_fit': self._getExtraPath("beamtilt_delta-phase_lin-fit_optics-group_%(og)d.mrc:mrc"),
+            'trefoil_it_fit': self._getExtraPath("beamtilt_delta-phase_iter-fit_optics-group_%(og)d_N-3.mrc:mrc"),
+            'trefoil_fit': self._getExtraPath("beamtilt_delta-phase_lin-fit_optics-group_%(og)d_N-3.mrc:mrc"),
+            'trefoil_residual_fit': self._getExtraPath("beamtilt_delta-phase_lin-fit_optics-group_%(og)d_N-3_residual.mrc:mrc"),
+            'beamtilt_obs': self._getExtraPath("beamtilt_delta-phase_per-pixel_optics-group_%(og)d.mrc:mrc")
+        }
+
+        self._updateFilenamesDict(myDict)
+
+    # -------------------------- DEFINE param functions -----------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputParticles', params.PointerParam,
@@ -164,10 +191,13 @@ def _defineParams(self, form):
 
     # -------------------------- STEPS functions ------------------------------
     def _insertAllSteps(self):
+        self._initialize()
         self._insertFunctionStep('convertInputStep')
         self._insertFunctionStep('refineCtfStep')
         self._insertFunctionStep('createOutputStep')
-        self._insertFunctionStep('createGlobalInfoStep')
+
+        if self.doCtfFitting:
+            self._insertFunctionStep('createGlobalInfoStep')
 
     def convertInputStep(self):
         inputParts = self.inputParticles.get()
@@ -225,7 +255,7 @@ def createOutputStep(self):
         imgSet = self.inputParticles.get()
         outImgSet = self._createSetOfParticles()
         outImgSet.copyInfo(imgSet)
-        outImgsFn = self.fileWithRefinedCTFName()
+        outImgsFn = self._getFileName("output_star")
         imgSet.setAlignmentProj()
 
         if self.IS_GT30():
@@ -255,7 +285,7 @@ def createGlobalInfo(self, filename):
         return ctfInfo
 
     def createGlobalInfoStep(self):
-        self.createGlobalInfo(self.fileWithAnalyzeInfo())
+        self.createGlobalInfo(self._getFileName("ctf_sqlite"))
 
     def _updateItem30(self, particle, row):
         particle.setCTF(convert.rowToCtfModel(row))
@@ -279,19 +309,5 @@ def _validate(self):
         errors = []
         return errors
 
-    def fileWithRefinedCTFName(self):
-        return self._getExtraPath('particles_ctf_refine.star')
-
-    def fileWithPhaseDifferenceName(self):
-        return self._getExtraPath(
-            'beamtilt_delta-phase_per-pixel_class_0.mrc:mrc')
-
-    def fileWithModelFitterName(self):
-        return self._getExtraPath(
-            'beamtilt_delta-phase_lin-fit_class_0.mrc:mrc')
-
-    def fileWithAnalyzeInfo(self):
-        return self._getExtraPath('ctf_analyze.sqlite')
-
     def IS_GT30(self):
         return relion.Plugin.IS_GT30()

relion/protocols/protocol_extract_particles.py

@@ -360,7 +360,7 @@ def _updateOutputSet(self, outputName, outputSet,
         if first:
             og = OpticsGroups.fromImages(outputSet)
             og.updateAll(rlnImagePixelSize=self._getNewSampling(),
-                         rlnImageSize=self.getBoxSize())
+                         rlnImageSize=self.getNewImgSize())
             og.toImages(outputSet)
 
         ProtExtractParticles._updateOutputSet(self, outputName, outputSet,
@@ -430,9 +430,8 @@ def getBoxScale(self):
         f = self.getScaleFactor()
         return f / self._getDownFactor() if self._doDownsample() else f
 
-    def getBoxSize(self):
-        # This function is needed by the wizard
-        return int(self.getCoords().getBoxSize() * self.getBoxScale())
+    def getNewImgSize(self):
+        return int(self.rescaledSize if self._doDownsample() else self.boxSize)
 
     def _getOutputImgMd(self):
         return self._getPath('images.xmd')

relion/tests/test_convert_relion.py

@@ -168,7 +168,8 @@ def launchTest(self, fileKey, mList, alignType=None, **kwargs):
                           sphericalAberration=2,
                           amplitudeContrast=0.1,
                           magnification=60000)
-        og = OpticsGroups.create(rlnMtfFileName="mtfFile1.star")
+        og = OpticsGroups.create(rlnMtfFileName="mtfFile1.star",
+                                 rlnImageSize=128)
         partSet.setSamplingRate(1.0)
         partSet.setAcquisition(acq)
         og.toImages(partSet)
@@ -640,7 +641,7 @@ def test_reconstRotandShiftFlip(self):
         a3 -> flip a2
         a4 -> flip a2 (identical to a3)
         """
-        # in practice this is irrelevant since no converson with |mat|==-1
+        # in practice this is irrelevant since no conversion with |mat|==-1
         mList = [
                  [[1., 0., 0., 0.],  # a1
                   [0., 1., 0., 0.],