Repository: I2PC/scipion-em-xmipp (language: python)

Issue:
FlexAlign: fixing default number of control points



Discussion / hints:

Fix patch:
diff --git a/setup.py b/setup.py
index 9a10569c9..25bd031ea 100644
--- a/setup.py
+++ b/setup.py
@@ -46,7 +46,7 @@
     # For a discussion on single-sourcing the version across setup.py and the
     # project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='21.06.1',  # Required
+    version='21.9.0',  # Required
 
     # This is a one-line description or tagline of what your project does. This
     # corresponds to the "Summary" metadata field:
diff --git a/xmipp3/__init__.py b/xmipp3/__init__.py
index d983fc59d..905252184 100644
--- a/xmipp3/__init__.py
+++ b/xmipp3/__init__.py
@@ -39,7 +39,7 @@
 
 _logo = "xmipp_logo.png"
 _references = ['delaRosaTrevin2013', 'Sorzano2013']
-_currentVersion = '3.21.06.1'
+_currentVersion = '3.21.09.0'
 
 
 class Plugin(pwem.Plugin):
diff --git a/xmipp3/bibtex.py b/xmipp3/bibtex.py
index af51b9672..2cef143f9 100644
--- a/xmipp3/bibtex.py
+++ b/xmipp3/bibtex.py
@@ -591,5 +591,18 @@
 ARTICLE-NUMBER = {xxx},
 URL = {https://github.com/I2PC/xmipp/wiki/FSO---Fourier-Shell-Occupancy},
 }
+
+@article{Fernandez-Gimenez2021,
+title = {Cryo-EM density maps adjustment for subtraction, consensus and sharpening},
+journal = {Journal of Structural Biology},
+volume = {213},
+number = {4},
+year = {2021},
+issn = {1047-8477},
+doi = {https://doi.org/10.1016/j.jsb.2021.107780},
+url = {https://www.sciencedirect.com/science/article/pii/S104784772100085X},
+author = {E. Fernández-Giménez and M. Martínez and R. Sánchez-García and R. Marabini and E. Ramírez-Aportela and
+P. Conesa and J.M. Carazo and C.O.S. Sorzano},
+}
 """
 
diff --git a/xmipp3/protocols/__init__.py b/xmipp3/protocols/__init__.py
index 3cc3957cc..c753c8be9 100644
--- a/xmipp3/protocols/__init__.py
+++ b/xmipp3/protocols/__init__.py
@@ -117,14 +117,18 @@
 from .protocol_metaprotocol_create_subset import XmippMetaProtCreateSubset
 from .protocol_metaprotocol_create_output import XmippMetaProtCreateOutput
 from .protocol_metaprotocol_discrete_heterogeneity_scheduler import XmippMetaProtDiscreteHeterogeneityScheduler # keep this after XmippMetaProtCreateSubset and XmippMetaProtCreateOutput
-# from .protocol_angular_alignment_sph import XmippProtAngularAlignmentSPH
-# from .protocol_volume_deform_sph import XmippProtVolumeDeformSPH
-# from .protocol_structure_map_sph import XmippProtStructureMapSPH
+from .protocol_angular_alignment_zernike3d import XmippProtAngularAlignmentZernike3D
+from .protocol_volume_deform_zernike3d import XmippProtVolumeDeformZernike3D
+from .protocol_structure_map_zernike3d import XmippProtStructureMapZernike3D
 from .protocol_align_volume_and_particles import XmippProtAlignVolumeParticles
 from .protocol_local_ctf import XmippProtLocalCTF
 from .protocol_analyze_local_ctf import XmippProtAnalyzeLocalCTF
 from .protocol_consensus_local_ctf import XmippProtConsensusLocalCTF
 from .protocol_particle_pick_remove_duplicates import XmippProtPickingRemoveDuplicates
 from .protocol_metaprotocol_golden_highres import XmippMetaProtGoldenHighRes
+# from .protocol_apply_deformation_zernike3d import XmippProtApplyZernike3D
+# from .protocol_kmeans_clustering import XmippProtKmeansSPH
+from .protocol_structure_map import XmippProtStructureMap
+from .protocol_pdb_deform_zernike3d import XmippProtDeformPDBZernike3D
 from .protocol_volume_adjust_sub import XmippProtVolAdjust, XmippProtVolSubtraction
-from .protocol_volume_consensus import XmippProtVolConsensus
+from .protocol_volume_consensus import XmippProtVolConsensus
\ No newline at end of file
diff --git a/xmipp3/protocols/protocol_align_volume.py b/xmipp3/protocols/protocol_align_volume.py
index 6ea0ae89d..858b1d418 100644
--- a/xmipp3/protocols/protocol_align_volume.py
+++ b/xmipp3/protocols/protocol_align_volume.py
@@ -318,7 +318,7 @@ def _getMaskArgs(self):
             if self.maskType == ALIGN_MASK_CIRCULAR:
                 maskArgs+=" --mask circular -%d" % self.maskRadius
             else:
-                maskArgs+=" --mask binary_file %s" % self.volMask
+                maskArgs+=" --mask binary_file %s" % self.maskFile.get().getFileName()
         return maskArgs
     
     def _getAlignArgs(self):
diff --git a/xmipp3/protocols/protocol_angular_alignment_sph.py b/xmipp3/protocols/protocol_angular_alignment_sph.py
deleted file mode 100644
index 8b6f9770e..000000000
--- a/xmipp3/protocols/protocol_angular_alignment_sph.py
+++ /dev/null
@@ -1,196 +0,0 @@
-########
-
-# This protocol has been removed from Scipion as there is a newer version being developed that
-# will solve some of the problems of this current implementation. The code will be available
-# again in a newer release.
-
-########
-
-
-# **************************************************************************
-# *
-# * Authors:     Amaya Jimenez Moreno (ajimenez@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 pyworkflow.protocol.params as params
-# from pyworkflow.utils.path import moveFile
-#
-# from pwem.protocols import ProtAnalysis3D
-# import pwem.emlib.metadata as md
-# from pwem.emlib.image import ImageHandler
-# from pwem.constants import ALIGN_PROJ
-#
-# from xmipp3.convert import (writeSetOfParticles, createItemMatrix,
-#                             setXmippAttributes)
-# from xmipp3.base import writeInfoField, readInfoField
-# import numpy as np
-#
-# from pyworkflow import VERSION_2_0
-#
-#
-# class XmippProtAngularAlignmentSPH(ProtAnalysis3D):
-#     """ Protocol for flexible angular alignment based on spherical harmonics. """
-#     _label = 'sph angular align'
-#     _lastUpdateVersion = VERSION_2_0
-#
-#     # --------------------------- DEFINE param functions --------------------------------------------
-#     def _defineParams(self, form):
-#         form.addSection(label='Input')
-#         form.addParam('inputParticles', params.PointerParam, label="Input particles", pointerClass='SetOfParticles')
-#         form.addParam('inputVolume', params.PointerParam, label="Input volume", pointerClass='Volume')
-#         form.addParam('targetResolution', params.FloatParam, label="Target resolution (A)", default=8.0,
-#                       help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
-#                            "2/3 of the Fourier spectrum.")
-#         form.addParam('depth', params.IntParam, default=3,
-#                       label='Harmonical depth', expertLevel=params.LEVEL_ADVANCED,
-#                       help='Harmonical depth of the deformation=1,2,3,...')
-#         form.addParam('maxShift', params.FloatParam, default=-1,
-#                       label='Maximum shift (px)', expertLevel=params.LEVEL_ADVANCED,
-#                       help='Maximum shift allowed in pixels')
-#         form.addParam('maxAngular', params.FloatParam, default=5,
-#                       label='Maximum angular change (degrees)', expertLevel=params.LEVEL_ADVANCED,
-#                       help='Maximum angular change allowed (in degrees)')
-#         form.addParam('maxResolution', params.FloatParam, default=4.0,
-#                       label='Maximum resolution (A)', expertLevel=params.LEVEL_ADVANCED,
-#                       help='Maximum resolution (A)')
-#         form.addParallelSection(threads=1, mpi=8)
-#
-#     def _createFilenameTemplates(self):
-#         """ Centralize how files are called """
-#         myDict = {
-#             'imgsFn': self._getExtraPath('input_particles.xmd'),
-#             'fnVol': self._getExtraPath('input_volume.vol'),
-#             'fnOut': self._getExtraPath('output_particles.xmd'),
-#             'fnOutDir': self._getExtraPath()
-#                  }
-#         self._updateFilenamesDict(myDict)
-#
-#     # --------------------------- INSERT steps functions --------------------------------------------
-#     def _insertAllSteps(self):
-#         self._createFilenameTemplates()
-#         self._insertFunctionStep("convertStep")
-#         self._insertFunctionStep("alignmentStep")
-#         self._insertFunctionStep("createOutputStep")
-#
-#     # --------------------------- STEPS functions ---------------------------------------------------
-#     def convertStep(self):
-#         imgsFn = self._getFileName('imgsFn')
-#         fnVol = self._getFileName('fnVol')
-#
-#         inputParticles = self.inputParticles.get()
-#         writeSetOfParticles(inputParticles, imgsFn)
-#         Xdim = inputParticles.getXDim()
-#         Ts = inputParticles.getSamplingRate()
-#         newTs = self.targetResolution.get() * 1.0 /3.0
-#         newTs = max(Ts, newTs)
-#         self.newXdim = int(Xdim * Ts / newTs)
-#         writeInfoField(self._getExtraPath(), "sampling", md.MDL_SAMPLINGRATE, newTs)
-#         writeInfoField(self._getExtraPath(), "size", md.MDL_XSIZE, self.newXdim)
-#         if self.newXdim != Xdim:
-#             self.runJob("xmipp_image_resize",
-#                         "-i %s -o %s --save_metadata_stack %s --fourier %d" %
-#                         (imgsFn,
-#                          self._getExtraPath('scaled_particles.stk'),
-#                          self._getExtraPath('scaled_particles.xmd'),
-#                          self.newXdim), numberOfMpi=1)
-#             moveFile(self._getExtraPath('scaled_particles.xmd'), imgsFn)
-#
-#         ih = ImageHandler()
-#         ih.convert(self.inputVolume.get(), fnVol)
-#         Xdim = self.inputVolume.get().getDim()[0]
-#         if Xdim != self.newXdim:
-#             self.runJob("xmipp_image_resize"
-#                         ,"-i %s --dim %d " %(fnVol, self.newXdim), numberOfMpi=1)
-#
-#
-#     def alignmentStep(self):
-#         imgsFn = self._getFileName('imgsFn')
-#         fnVol = self._getFileName('fnVol')
-#         fnOut =  self._getFileName('fnOut')
-#         fnOutDir = self._getFileName('fnOutDir')
-#         Ts = readInfoField(self._getExtraPath(), "sampling", md.MDL_SAMPLINGRATE)
-#         params = ' -i %s --ref %s -o %s --optimizeAlignment --optimizeDeformation ' \
-#                  '--depth %d --max_shift %f --max_angular_change %f --sampling %f ' \
-#                  ' --max_resolution %f --odir %s --resume' %\
-#                  (imgsFn, fnVol, fnOut, self.depth, self.maxShift, self.maxAngular,
-#                   Ts, self.maxResolution, fnOutDir)
-#         if self.inputParticles.get().isPhaseFlipped(): #preguntar
-#             params += ' --phaseFlipped'
-#         self.runJob("xmipp_angular_sph_alignment", params, numberOfMpi=self.numberOfMpi.get())
-#
-#
-#     def createOutputStep(self):
-#         from sklearn.manifold import TSNE
-#         fnOut = self._getFileName('fnOut')
-#         mdOut = md.MetaData(fnOut)
-#         i = 0
-#         for row in md.iterRows(mdOut):
-#             coeffs = mdOut.getValue(md.MDL_SPH_COEFFICIENTS, row.getObjId())
-#             if i==0:
-#                 coeffMatrix = coeffs
-#             else:
-#                 coeffMatrix = np.vstack((coeffMatrix, coeffs))
-#             i+=1
-#         X_tsne_1d = TSNE(n_components=1).fit_transform(coeffMatrix)
-#         X_tsne_2d = TSNE(n_components=2).fit_transform(coeffMatrix)
-#
-#         newMdOut = md.MetaData()
-#         i=0
-#         for row in md.iterRows(mdOut):
-#             newRow = row
-#             newRow.setValue(md.MDL_SPH_TSNE_COEFF1D, float(X_tsne_1d[i,0]))
-#             newRow.setValue(md.MDL_SPH_TSNE_COEFF2D, [float(X_tsne_2d[i, 0]),  float(X_tsne_2d[i, 1])])
-#             newRow.addToMd(newMdOut)
-#             i+=1
-#             newMdOut.write(fnOut)
-#
-#         inputSet = self.inputParticles.get()
-#         partSet = self._createSetOfParticles()
-#
-#         partSet.copyInfo(inputSet)
-#         partSet.setAlignmentProj()
-#         partSet.copyItems(inputSet,
-#                           updateItemCallback=self._updateParticle,
-#                           itemDataIterator=md.iterRows(fnOut, sortByLabel=md.MDL_ITEM_ID))
-#
-#         self._defineOutputs(outputParticles=partSet)
-#         self._defineTransformRelation(self.inputParticles, partSet)
-#
-#
-# # --------------------------- UTILS functions --------------------------------------------
-#     def _updateParticle(self, item, row):
-#         setXmippAttributes(item, row, md.MDL_ANGLE_ROT, md.MDL_ANGLE_TILT,
-#                            md.MDL_ANGLE_PSI, md.MDL_SHIFT_X, md.MDL_SHIFT_Y,
-#                            md.MDL_FLIP, md.MDL_SPH_DEFORMATION,
-#                            md.MDL_SPH_COEFFICIENTS, md.MDL_SPH_TSNE_COEFF1D,
-#                            md.MDL_SPH_TSNE_COEFF2D)
-#         createItemMatrix(item, row, align=ALIGN_PROJ)
-#
-#     def getInputParticles(self):
-#         return self.inputParticles.get()
-
-
-
-
-
diff --git a/xmipp3/protocols/protocol_angular_alignment_zernike3d.py b/xmipp3/protocols/protocol_angular_alignment_zernike3d.py
new file mode 100644
index 000000000..345d7a353
--- /dev/null
+++ b/xmipp3/protocols/protocol_angular_alignment_zernike3d.py
@@ -0,0 +1,251 @@
+
+# **************************************************************************
+# *
+# * Authors:     Amaya Jimenez Moreno (ajimenez@cnb.csic.es)
+# *              David Herreros Calero (dherreos@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
+
+import pyworkflow.protocol.params as params
+from pyworkflow.object import Integer
+from pyworkflow.utils.path import moveFile
+from pyworkflow import VERSION_2_0
+
+from pwem.protocols import ProtAnalysis3D
+import pwem.emlib.metadata as md
+from pwem.emlib.image import ImageHandler
+from pwem.constants import ALIGN_PROJ
+
+from xmipp3.convert import (writeSetOfParticles, createItemMatrix,
+                            setXmippAttributes)
+from xmipp3.base import writeInfoField, readInfoField
+
+
+class XmippProtAngularAlignmentZernike3D(ProtAnalysis3D):
+    """ Protocol for flexible angular alignment based on Zernike3D basis. """
+    _label = 'angular align - Zernike3D'
+    _lastUpdateVersion = VERSION_2_0
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addHidden(params.USE_GPU, params.BooleanParam, default=True,
+                       label="Use GPU for execution",
+                       help="This protocol has both CPU and GPU implementation.\
+                           Select the one you want to use.")
+        form.addHidden(params.GPU_LIST, params.StringParam, default='0',
+                       expertLevel=params.LEVEL_ADVANCED,
+                       label="Choose GPU IDs",
+                       help="Add a list of GPU devices that can be used")
+        form.addParam('inputParticles', params.PointerParam, label="Input particles", pointerClass='SetOfParticles')
+        form.addParam('inputVolume', params.PointerParam, label="Input volume", pointerClass='Volume')
+        form.addParam('inputVolumeMask', params.PointerParam, label="Input volume mask", pointerClass='VolumeMask',
+                      allowsNull=True)
+        form.addParam('targetResolution', params.FloatParam, label="Target resolution (A)", default=8.0,
+                      help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
+                           "2/3 of the Fourier spectrum.")
+        form.addParam('l1', params.IntParam, default=3,
+                      label='Zernike Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Zernike Polynomials of the deformation=1,2,3,...')
+        form.addParam('l2', params.IntParam, default=2,
+                      label='Harmonical Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Spherical Harmonics of the deformation=1,2,3,...')
+        form.addParam('maxShift', params.FloatParam, default=-1,
+                      label='Maximum shift (px)', expertLevel=params.LEVEL_ADVANCED,
+                      help='Maximum shift allowed in pixels')
+        form.addParam('maxAngular', params.FloatParam, default=5,
+                      label='Maximum angular change (degrees)', expertLevel=params.LEVEL_ADVANCED,
+                      help='Maximum angular change allowed (in degrees)')
+        form.addParam('maxResolution', params.FloatParam, default=4.0,
+                      label='Maximum resolution (A)', expertLevel=params.LEVEL_ADVANCED,
+                      help='Maximum resolution (A)')
+        form.addParam('regularization', params.FloatParam, default=0.005, label='Regularization',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Penalization to deformations (higher values penalize more the deformation).')
+        form.addParam('ignoreCTF', params.BooleanParam, default=False, label='Ignore CTF?',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help="If true, volume projection won't be subjected to CTF corrections")
+        form.addParam('optimizeAlignment', params.BooleanParam, default=True, label='Optimize alignment?',
+                     expertLevel=params.LEVEL_ADVANCED)
+        form.addParallelSection(threads=1, mpi=8)
+
+    def _createFilenameTemplates(self):
+        """ Centralize how files are called """
+        myDict = {
+            'imgsFn': self._getExtraPath('input_particles.xmd'),
+            'fnVol': self._getExtraPath('input_volume.vol'),
+            'fnVolMask': self._getExtraPath('input_volume_mask.vol'),
+            'fnOut': self._getExtraPath('output_particles.xmd'),
+            'fnOutDir': self._getExtraPath()
+                 }
+        self._updateFilenamesDict(myDict)
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._createFilenameTemplates()
+        self._insertFunctionStep("convertStep")
+        self._insertFunctionStep("alignmentStep")
+        self._insertFunctionStep("createOutputStep")
+
+    # --------------------------- STEPS functions ---------------------------------------------------
+    def convertStep(self):
+        imgsFn = self._getFileName('imgsFn')
+        fnVol = self._getFileName('fnVol')
+        fnVolMask = self._getFileName('fnVolMask')
+
+        inputParticles = self.inputParticles.get()
+        writeSetOfParticles(inputParticles, imgsFn)
+        Xdim = inputParticles.getXDim()
+        self.Ts = inputParticles.getSamplingRate()
+        newTs = self.targetResolution.get() * 1.0 /3.0
+        self.newTs = max(self.Ts, newTs)
+        self.newXdim = int(Xdim * self.Ts / newTs)
+        writeInfoField(self._getExtraPath(), "sampling", md.MDL_SAMPLINGRATE, newTs)
+        writeInfoField(self._getExtraPath(), "size", md.MDL_XSIZE, self.newXdim)
+        if self.newXdim != Xdim:
+            self.runJob("xmipp_image_resize",
+                        "-i %s -o %s --save_metadata_stack %s --fourier %d" %
+                        (imgsFn,
+                         self._getExtraPath('scaled_particles.stk'),
+                         self._getExtraPath('scaled_particles.xmd'),
+                         self.newXdim), numberOfMpi=1)
+            moveFile(self._getExtraPath('scaled_particles.xmd'), imgsFn)
+
+        ih = ImageHandler()
+        ih.convert(self.inputVolume.get(), fnVol)
+        Xdim = self.inputVolume.get().getDim()[0]
+        if Xdim != self.newXdim:
+            self.runJob("xmipp_image_resize",
+                        "-i %s --dim %d " % (fnVol, self.newXdim), numberOfMpi=1)
+        if self.inputVolumeMask.get():
+            ih.convert(self.inputVolumeMask.get(), fnVolMask)
+            if Xdim != self.newXdim:
+                self.runJob("xmipp_image_resize",
+                            "-i %s --dim %d --interp nearest" % (fnVolMask, self.newXdim), numberOfMpi=1)
+
+    def alignmentStep(self):
+        imgsFn = self._getFileName('imgsFn')
+        fnVol = self._getFileName('fnVol')
+        fnOut = self._getFileName('fnOut')
+        fnVolMask = self._getFileName('fnVolMask')
+        fnOutDir = self._getFileName('fnOutDir')
+        Ts = readInfoField(self._getExtraPath(), "sampling", md.MDL_SAMPLINGRATE)
+        params = ' -i %s --ref %s -o %s --optimizeDeformation --optimizeDefocus ' \
+                 '--l1 %d --l2 %d --max_shift %f --max_angular_change %f --sampling %f ' \
+                 ' --max_resolution %f --odir %s --resume --regularization %f' %\
+                 (imgsFn, fnVol, fnOut, self.l1.get(), self.l2.get(), self.maxShift,
+                  self.maxAngular, Ts, self.maxResolution, fnOutDir, self.regularization.get())
+        if self.optimizeAlignment.get():
+            params += ' --optimizeAlignment'
+        if self.ignoreCTF.get():
+            params += ' --ignoreCTF'
+        if self.inputParticles.get().isPhaseFlipped():
+            params += ' --phaseFlipped'
+        if self.inputVolumeMask.get():
+            params += ' --mask %s' % fnVolMask
+
+        if self.useGpu.get():
+            params += ' --device %d' % self.getGpuList()[0]
+            program = 'xmipp_cuda_angular_sph_alignment'
+            self.runJob(program, params)
+        else:
+            program = 'xmipp_angular_sph_alignment'
+            self.runJob(program, params, numberOfMpi=self.numberOfMpi.get())
+
+
+    def createOutputStep(self):
+        from sklearn.manifold import TSNE
+        Xdim = self.inputParticles.get().getXDim()
+        self.Ts = self.inputParticles.get().getSamplingRate()
+        newTs = self.targetResolution.get() * 1.0 /3.0
+        self.newTs = max(self.Ts, newTs)
+        self.newXdim = int(Xdim * self.Ts / newTs)
+        fnOut = self._getFileName('fnOut')
+        mdOut = md.MetaData(fnOut)
+        coeffMatrix = np.vstack(mdOut.getColumnValues(md.MDL_SPH_COEFFICIENTS))
+        X_tsne_1d = TSNE(n_components=1).fit_transform(coeffMatrix)
+        X_tsne_2d = TSNE(n_components=2).fit_transform(coeffMatrix)
+
+        newMdOut = md.MetaData()
+        i = 0
+        for row in md.iterRows(mdOut):
+            newRow = row
+            newRow.setValue(md.MDL_SPH_TSNE_COEFF1D, float(X_tsne_1d[i, 0]))
+            newRow.setValue(md.MDL_SPH_TSNE_COEFF2D, [float(X_tsne_2d[i, 0]),  float(X_tsne_2d[i, 1])])
+            if self.newTs != self.Ts:
+                coeffs = mdOut.getValue(md.MDL_SPH_COEFFICIENTS, row.getObjId())
+                correctionFactor = self.inputVolume.get().getDim()[0] / self.newXdim
+                coeffs = [correctionFactor*coeff for coeff in coeffs]
+                newRow.setValue(md.MDL_SPH_COEFFICIENTS, coeffs)
+            newRow.addToMd(newMdOut)
+            i += 1
+        newMdOut.write(fnOut)
+
+        inputSet = self.inputParticles.get()
+        partSet = self._createSetOfParticles()
+
+        partSet.copyInfo(inputSet)
+        partSet.setAlignmentProj()
+        partSet.copyItems(inputSet,
+                          updateItemCallback=self._updateParticle,
+                          itemDataIterator=md.iterRows(fnOut, sortByLabel=md.MDL_ITEM_ID))
+        partSet.L1 = Integer(self.l1.get())
+        partSet.L2 = Integer(self.l2.get())
+        partSet.Rmax = Integer(self.inputVolume.get().getDim()[0] / 2)
+
+        self._defineOutputs(outputParticles=partSet)
+        self._defineTransformRelation(self.inputParticles, partSet)
+
+
+# --------------------------- UTILS functions --------------------------------------------
+    def _updateParticle(self, item, row):
+        setXmippAttributes(item, row, md.MDL_ANGLE_ROT, md.MDL_ANGLE_TILT,
+                           md.MDL_ANGLE_PSI, md.MDL_SHIFT_X, md.MDL_SHIFT_Y,
+                           md.MDL_FLIP, md.MDL_SPH_DEFORMATION,
+                           md.MDL_SPH_COEFFICIENTS, md.MDL_SPH_TSNE_COEFF1D,
+                           md.MDL_SPH_TSNE_COEFF2D)
+        createItemMatrix(item, row, align=ALIGN_PROJ)
+
+    def getInputParticles(self):
+        return self.inputParticles.get()
+
+    # ----------------------- VALIDATE functions ----------------------------------------
+    def validate(self):
+        """ Try to find errors on define params. """
+        errors = []
+        l1 = self.l1.get()
+        l2 = self.l2.get()
+        if (l1 - l2) < 0:
+            errors.append('Zernike degree must be higher than '
+                          'SPH degree.')
+        return errors
+
+
+
+
+
diff --git a/xmipp3/protocols/protocol_apply_deformation_zernike3d.py b/xmipp3/protocols/protocol_apply_deformation_zernike3d.py
new file mode 100644
index 000000000..fbe719722
--- /dev/null
+++ b/xmipp3/protocols/protocol_apply_deformation_zernike3d.py
@@ -0,0 +1,112 @@
+# # *************************************************
+# # * This protocol will be offered in future releases (more testing is needed)
+# # *************************************************
+
+# # **************************************************************************
+# # *
+# # * Authors:     David Herreros Calero (dherreros@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 os, glob
+# import numpy as np
+#
+# from pwem.protocols import ProtAnalysis3D
+# from pwem.objects import Volume, Integer
+# import pwem.emlib.metadata as md
+# import pyworkflow.protocol.params as params
+# import pyworkflow.utils as pwutils
+#
+#
+# class XmippProtApplyZernike3D(ProtAnalysis3D):
+#     """ Protocol to apply the deformation computed through spherical harmonics to
+#     EM maps. """
+#     _label = 'apply deformation - Zernike3D'
+#
+#     # --------------------------- DEFINE param functions --------------------------------------------
+#     def _defineParams(self, form):
+#         form.addSection(label='Input')
+#         form.addParam('inputVol', params.PointerParam, label="Input volume",
+#                       pointerClass='Volume', important=True,
+#                       help='Select a Volume to be deformed.')
+#         form.addParam('inputClasses', params.PointerParam, pointerClass='SetOfClasses2D',
+#                       label="Input Coefficients", important=True,
+#                       help='Specify a path to the deformation coefficients metadata.')
+#
+#     # --------------------------- INSERT steps functions -----------------------
+#     def _insertAllSteps(self):
+#         self._insertFunctionStep("deformStep")
+#         self._insertFunctionStep("createOutputStep")
+#
+#     # --------------------------- STEPS functions ------------------------------
+#     def deformStep(self):
+#         self.samplingRate_Volume = self.inputVol.get().getSamplingRate()
+#         self.outParams = []
+#         classes = self.inputClasses.get()
+#         self.samplingRate_Coefficients = classes.getSamplingRate()
+#         correctionFactor = self.samplingRate_Coefficients / self.samplingRate_Volume
+#         for rep in classes.iterItems():
+#             basisParams = [rep.L1, rep.L2, rep.Rmax]
+#             coeffs = rep.getRepresentative().get()
+#             coeffs = np.fromstring(coeffs, sep=',')
+#             coeffs = correctionFactor * coeffs
+#             idx = rep.getObjId() + 1
+#             file = self._getTmpPath('coeffs.txt')
+#             # np.savetxt(file, coeffs)
+#             with open(file, 'w') as fid:
+#                 fid.write(' '.join(map(str, basisParams)) + "\n")
+#                 fid.write(' '.join(map(str, coeffs)) + "\n")
+#             outFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.vol' % idx
+#             params = ' -i %s --clnm %s -o %s' % \
+#                      (self.inputVol.get().getFileName(), file, self._getExtraPath(outFile))
+#             self.runJob("xmipp_volume_apply_deform_sph", params)
+#             params = ' -i %s --sampling_rate %f' % (self._getExtraPath(outFile), self.samplingRate_Volume)
+#             self.runJob("xmipp_image_header", params)
+#             self.outParams.append((self._getExtraPath(outFile), rep.getSize()))
+#
+#     def createOutputStep(self):
+#         classes3DSet = self._createSetOfClasses3D(self.inputClasses.get().getImages())
+#         classes3DSet.classifyItems(updateItemCallback=self._updateParticle,
+#                                    updateClassCallback=self._updateClass,
+#                                    itemDataIterator=self.iterClassesId())
+#         result = {'outputClasses': classes3DSet}
+#         self._defineOutputs(**result)
+#         self._defineSourceRelation(self.inputClasses, classes3DSet)
+#         self._defineOutputs(**result)
+#         self._defineSourceRelation(self.inputVol, classes3DSet)
+#
+#     # ------------------------------- UTILS functions -------------------------------
+#     def _updateParticle(self, item, idc):
+#         item.setClassId(idc)
+#
+#     def _updateClass(self, item):
+#         representative = item.getRepresentative()
+#         volumeFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.vol' \
+#                      % (item.getObjId() + 1)
+#         volumeFile = self._getExtraPath(volumeFile)
+#         representative.setSamplingRate(self.samplingRate_Volume)
+#         representative.setLocation(volumeFile)
+#
+#     def iterClassesId(self):
+#         for img in self.inputClasses.get().iterClassItems():
+#             yield img.getClassId()
diff --git a/xmipp3/protocols/protocol_deep_align.py b/xmipp3/protocols/protocol_deep_align.py
index 7a9fb0bb9..0cee544df 100644
--- a/xmipp3/protocols/protocol_deep_align.py
+++ b/xmipp3/protocols/protocol_deep_align.py
@@ -53,6 +53,14 @@ class XmippProtDeepAlign(ProtRefine3D, xmipp3.XmippProtocol):
     _lastUpdateVersion = VERSION_3_0
     _conda_env = 'xmipp_DLTK_v0.3'
     _ih = ImageHandler()
+    _cond_modelPretrainTrue = 'modelPretrain==True'
+    _cond_modelPretrainFalse = 'modelPretrain==False'
+    _fnCorrectedParticlesStk = 'corrected_particles.stk'    
+    _fnCorrectedParticlesXmd = 'corrected_particles.xmd'
+    _fnVolumeVol = 'volume.vol'
+    _fnConeCenterDoc = 'coneCenters.doc'
+    _tempNumXmd = '%d.xmd'
+
 
     def __init__(self, **args):
         ProtRefine3D.__init__(self, **args)
@@ -78,7 +86,7 @@ def _defineParams(self, form):
                            'If you choose "no" new models will be trained.')
         form.addParam('pretrainedModels', PointerParam,
                       pointerClass=self.getClassName(),
-                      condition='modelPretrain==True',
+                      condition=self._cond_modelPretrainTrue,
                       label='Set pretrained models',
                       help='Choose the protocol where your models were trained. '
                            'Be careful with using proper models for your new prediction.')
@@ -87,7 +95,7 @@ def _defineParams(self, form):
                       pointerClass='SetOfParticles',
                       pointerCondition='hasAlignmentProj',
                       help='The set of particles previously aligned to be used as training set',
-                      condition='modelPretrain==False')
+                      condition=self._cond_modelPretrainFalse)
         form.addParam('targetResolution', FloatParam, label="Target resolution",
                       default=3.0,
                       help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
@@ -99,17 +107,17 @@ def _defineParams(self, form):
                       label="Number of epochs for training",
                       default=10,
                       help="Number of epochs for training.",
-                      condition='modelPretrain==False')
+                      condition=self._cond_modelPretrainFalse)
         form.addParam('batchSize', IntParam,
                       label="Batch size for training",
                       default=128,
                       help="Batch size for training.",
-                      condition='modelPretrain==False')
+                      condition=self._cond_modelPretrainFalse)
         form.addParam('spanConesTilt', FloatParam,
                       label="Distance between region centers",
                       default=30,
                       help="Distance in degrees between region centers.",
-                      condition='modelPretrain==False')
+                      condition=self._cond_modelPretrainFalse)
         form.addParam('numConesSelected', IntParam,
                       label="Number of selected regions per image",
                       default=2,
@@ -170,8 +178,44 @@ def _insertAllSteps(self):
         self._insertFunctionStep("createOutputStep", prerequisites=[stepId])
 
     # --------------------------- STEPS functions ---------------------------------------------------
-    def convertStep(self):
+    
+    def _getProjectionsExp(self, num):
+        return self._getExtraPath('projectionsExp%d.xmd' % num)
+
+    def _getConePrediction(self):
+        return self._getExtraPath('conePrediction.txt')
+
+    def _getOutCone(self, num):
+        return 'outCone%d.xmd' % num
+
+    def _resize(self, Xdim, fnCorrected, prefix, fnTarget):
+        if self.newXdim != Xdim:
+            self.runJob("xmipp_image_resize",
+                        "-i %s -o %s --save_metadata_stack %s --fourier %d" %
+                        (fnCorrected,
+                            self._getExtraPath(prefix + '.stk'),
+                            self._getExtraPath(prefix + '.xmd'),
+                            self.newXdim), numberOfMpi=self.myMPI.get())
+            moveFile(self._getExtraPath(prefix + '.xmd'), fnTarget)
+
+    def _correctWiener(self, hasCTF, fn, Ts):
+        if hasCTF and self.applyCTF.get():
+            fnCorrectedStk = self._getExtraPath(self._fnCorrectedParticlesStk)
+            fnCorrected = self._getExtraPath(self._fnCorrectedParticlesXmd)
+            args = "-i %s -o %s --save_metadata_stack %s --keep_input_columns" % (
+                fn, fnCorrectedStk, fnCorrected)
+            args += " --sampling_rate %f --correct_envelope" % Ts
+            if self.inputSet.get().isPhaseFlipped():
+                args += " --phase_flipped"
+            self.runJob("xmipp_ctf_correct_wiener2d",
+                        args, numberOfMpi=self.myMPI.get())
+
+    def _removeCorrectedParticles(self):
+        if exists(self._getExtraPath(self._fnCorrectedParticlesStk)):
+            remove(self._getExtraPath(self._fnCorrectedParticlesStk))
+            remove(self._getExtraPath(self._fnCorrectedParticlesXmd))
 
+    def convertStep(self):
         if self.modelPretrain.get() is True:
             fnPreProtocol = self.pretrainedModels.get()._getExtraPath()
             preXDim = readInfoField(fnPreProtocol, "size", emlib.MDL_XSIZE)
@@ -184,46 +228,30 @@ def convertStep(self):
         Ts = inputParticles.getSamplingRate()
         row = getFirstRow(self.imgsFn)
         hasCTF = row.containsLabel(emlib.MDL_CTF_DEFOCUSU) or emlib.containsLabel(
-                    emlib.MDL_CTF_MODEL)
+            emlib.MDL_CTF_MODEL)
 
         fnCorrected = self.imgsFn
-        
-        if hasCTF and self.applyCTF.get() is True:
-            fnCorrectedStk = self._getExtraPath('corrected_particles.stk')
-            fnCorrected = self._getExtraPath('corrected_particles.xmd')
-            args = "-i %s -o %s --save_metadata_stack %s --keep_input_columns" % (
-                    self.imgsFn, fnCorrectedStk, fnCorrected)
-            args += " --sampling_rate %f --correct_envelope" % Ts
-            if self.inputSet.get().isPhaseFlipped():
-                args += " --phase_flipped"
-            self.runJob("xmipp_ctf_correct_wiener2d", args,numberOfMpi=self.myMPI.get())
+
+        self._correctWiener(hasCTF, self.imgsFn, Ts)
 
         Xdim = inputParticles.getXDim()
         newTs = self.targetResolution.get() * 1.0 / 3.0
         newTs = max(Ts, newTs)
-        if self.modelPretrain.get() is False:
-            self.newXdim = int(float(Xdim * Ts / newTs))
-        else:
-            self.newXdim = int(preXDim)
+
+        self.newXdim = int(preXDim) if self.modelPretrain.get() else int(
+            float(Xdim * Ts / newTs))
+
         self.firstMaxShift = int(round(self.newXdim / 10))
         writeInfoField(self._getExtraPath(), "sampling",
                        emlib.MDL_SAMPLINGRATE, newTs)
         writeInfoField(self._getExtraPath(), "size", emlib.MDL_XSIZE,
                        self.newXdim)
-        if self.newXdim != Xdim:
-            self.runJob("xmipp_image_resize",
-                        "-i %s -o %s --save_metadata_stack %s --fourier %d" %
-                        (fnCorrected,
-                         self._getExtraPath('scaled_particles.stk'),
-                         self._getExtraPath('scaled_particles.xmd'),
-                         self.newXdim), numberOfMpi=self.myMPI.get())
-            moveFile(self._getExtraPath('scaled_particles.xmd'), self.imgsFn)
 
-        if exists(self._getExtraPath('corrected_particles.stk')):
-            remove(self._getExtraPath('corrected_particles.stk'))
-            remove(self._getExtraPath('corrected_particles.xmd'))
+        self._resize(Xdim, fnCorrected, 'scaled_particles', self.imgsFn)
 
-        fnVol = self._getTmpPath("volume.vol")
+        self._removeCorrectedParticles()
+
+        fnVol = self._getTmpPath(self._fnVolumeVol)
         self._ih.convert(self.inputVolume.get(), fnVol)
         Xdim = self.inputVolume.get().getDim()[0]
         if Xdim != self.newXdim:
@@ -235,37 +263,19 @@ def convertStep(self):
         writeSetOfParticles(inputTrain, self.trainImgsFn)
         row = getFirstRow(self.trainImgsFn)
         hasCTF = row.containsLabel(emlib.MDL_CTF_DEFOCUSU) or emlib.containsLabel(
-                    emlib.MDL_CTF_MODEL)
+            emlib.MDL_CTF_MODEL)
 
         fnCorrected = self.trainImgsFn
 
-        if hasCTF and self.applyCTF.get() is True:
-            fnCorrectedStk = self._getExtraPath('corrected_particles.stk')
-            fnCorrected = self._getExtraPath('corrected_particles.xmd')
-            args = "-i %s -o %s --save_metadata_stack %s --keep_input_columns" % (
-                    self.trainImgsFn, fnCorrectedStk, fnCorrected)
-            args += " --sampling_rate %f --correct_envelope" % Ts
-            if self.inputSet.get().isPhaseFlipped():
-                args += " --phase_flipped"
-            self.runJob("xmipp_ctf_correct_wiener2d", args,numberOfMpi=self.myMPI.get())
+        self._correctWiener(hasCTF, self.trainImgsFn, Ts)
 
-        Xdim = inputTrain.getXDim()
-        if self.newXdim != Xdim:
-            self.runJob("xmipp_image_resize",
-                        "-i %s -o %s --save_metadata_stack %s --fourier %d" %
-                        (fnCorrected,
-                         self._getExtraPath('scaled_train_particles.stk'),
-                         self._getExtraPath('scaled_train_particles.xmd'),
-                         self.newXdim), numberOfMpi=self.myMPI.get())
-            moveFile(self._getExtraPath('scaled_train_particles.xmd'),
-                     self.trainImgsFn)
+        self._resize(Xdim, fnCorrected,
+                     'scaled_train_particles', self.trainImgsFn)
 
-        if exists(self._getExtraPath('corrected_particles.stk')):
-            remove(self._getExtraPath('corrected_particles.stk'))
-            remove(self._getExtraPath('corrected_particles.xmd'))
+        self._removeCorrectedParticles()
 
     def generateConeCenters(self, fn):
-        fnVol = self._getTmpPath("volume.vol")
+        fnVol = self._getTmpPath(self._fnVolumeVol)
         fnCenters = self._getExtraPath(fn + ".stk")
         fnCentersMd = self._getExtraPath(fn + ".doc")
         # if self.spanConesTilt.get()>30:
@@ -316,7 +326,7 @@ def computeTrainingSet(self, nameTrain):
 
         totalCones = self.generateConeCenters('coneCenters')
         self.numCones = totalCones
-        fnCentersMd = self._getExtraPath("coneCenters.doc")
+        fnCentersMd = self._getExtraPath(self._fnConeCenterDoc)
         mdCones = emlib.MetaData(fnCentersMd)
 
         auxList = []
@@ -340,13 +350,13 @@ def computeTrainingSet(self, nameTrain):
             auxList.append(numCone - 1)
             row.addToMd(mdCone)
         for i in range(totalCones):
-            fnTrain = self._getExtraPath(nameTrain + "%d.xmd" % (i + 1))
+            fnTrain = self._getExtraPath(nameTrain + self._tempNumXmd % (i + 1))
             mdList[i].write(fnTrain)
 
 
     def prepareImagesForTraining(self):
 
-        fnCentersMd = self._getExtraPath("coneCenters.doc")
+        fnCentersMd = self._getExtraPath(self._fnConeCenterDoc)
         mdCones = emlib.MetaData(fnCentersMd)
         span = self.spanConesTilt.get()
         counterCones = 0
@@ -384,7 +394,7 @@ def prepareImagesForTraining(self):
             counterCones = counterCones + 1
 
         if self.modelPretrain.get() is False:
-            fnToFilter = self._getExtraPath('projectionsExp%d.xmd' % (lastLabel))
+            fnToFilter = self._getProjectionsExp(lastLabel)
             self.runJob("xmipp_transform_filter", " -i %s --fourier low_pass %f" %
                         (fnToFilter, 0.15), numberOfMpi=self.myMPI.get())
 
@@ -392,7 +402,7 @@ def projectStep(self, numProj, iniRot, endRot, iniTilt, endTilt, fn, idx):
 
         newXdim = readInfoField(self._getExtraPath(), "size",
                                 emlib.MDL_XSIZE)
-        fnVol = self._getTmpPath("volume.vol")
+        fnVol = self._getTmpPath(self._fnVolumeVol)
 
         uniformProjectionsStr = """
 # XMIPP_STAR_1 *
@@ -427,17 +437,11 @@ def projectStep(self, numProj, iniRot, endRot, iniTilt, endTilt, fn, idx):
         cleanPattern(self._getExtraPath('uniformProjections*'))
 
     def generateExpImagesStep(self, Nimgs, nameProj, nameExp, label):
-
-        newXdim = readInfoField(self._getExtraPath(), "size",
-                                emlib.MDL_XSIZE)
-        fnProj = self._getExtraPath(nameProj + "%d.xmd" % label)
-        fnExp = self._getExtraPath(nameExp + "%d.xmd" % label)
+        fnProj = self._getExtraPath(nameProj + self._tempNumXmd % label)
+        fnExp = self._getExtraPath(nameExp + self._tempNumXmd % label)
         fnLabels = self._getExtraPath('labels.txt')
         mdIn = emlib.MetaData(fnProj)
         mdExp = emlib.MetaData()
-        maxPsi = 180
-        maxShift = round(newXdim / 10)
-        idx = 1
         NimgsMd = mdIn.size()
         Nrepeats = int(Nimgs / NimgsMd)
         # if Nrepeats<10:
@@ -446,6 +450,28 @@ def generateExpImagesStep(self, Nimgs, nameProj, nameExp, label):
         if (label == 1 and exists(fnLabels)):
             remove(fnLabels)
         fileLabels = open(fnLabels, "a")
+        self._processRows(label, fnExp, mdIn, mdExp, Nrepeats, fileLabels)
+        mdExp.write(fnExp)
+        fileLabels.close()
+        if (label - 1) > 0:
+            labelPrev = -1
+            for n in range(1, label):
+                if exists(self._getExtraPath(nameExp + self._tempNumXmd % (label - n))):
+                    labelPrev = label - n
+                    break
+            if labelPrev != -1:
+                lastFnExp = self._getExtraPath(
+                    nameExp + self._tempNumXmd % (labelPrev))
+                self.runJob("xmipp_metadata_utilities",
+                            " -i %s --set union %s -o %s " %
+                            (lastFnExp, fnExp, fnExp), numberOfMpi=1)
+        remove(fnProj)
+
+    def _processRows(self, label, fnExp, mdIn, mdExp, Nrepeats, fileLabels):
+        newXdim = readInfoField(self._getExtraPath(), "size",
+                                emlib.MDL_XSIZE)
+        maxPsi = 180
+        maxShift = round(newXdim / 10)
         for row in iterRows(mdIn):
             fnImg = row.getValue(emlib.MDL_IMAGE)
             myRow = row
@@ -453,7 +479,7 @@ def generateExpImagesStep(self, Nimgs, nameProj, nameExp, label):
             Xdim, Ydim, _, _ = I.getDimensions()
             Xdim2 = Xdim / 2
             Ydim2 = Ydim / 2
-            if Nrepeats==0:
+            if Nrepeats == 0:
                 myRow.addToMd(mdExp)
                 idx += 1
                 fileLabels.write(str(label - 1) + '\n')
@@ -468,7 +494,8 @@ def generateExpImagesStep(self, Nimgs, nameProj, nameExp, label):
                     M = np.float32([[c, s, (1 - c) * Xdim2 - s * Ydim2 + deltaX],
                                     [-s, c, s * Xdim2 + (1 - c) * Ydim2 + deltaY]])
                     newFn = ('%06d@' % idx) + fnExp[:-3] + 'stk'
-                    self._ih.applyTransform(fnImg, newFn, M, (Ydim, Xdim), doWrap=True)
+                    self._ih.applyTransform(
+                        fnImg, newFn, M, (Ydim, Xdim), doWrap=True)
 
                     myRow.setValue(emlib.MDL_IMAGE, newFn)
                     myRow.setValue(emlib.MDL_ANGLE_PSI, psiDeg)
@@ -477,24 +504,10 @@ def generateExpImagesStep(self, Nimgs, nameProj, nameExp, label):
                     myRow.addToMd(mdExp)
                     idx += 1
                     fileLabels.write(str(label - 1) + '\n')
-        mdExp.write(fnExp)
-        fileLabels.close()
-        if (label - 1) > 0:
-            labelPrev = -1
-            for n in range(1, label):
-                if exists(self._getExtraPath(nameExp + "%d.xmd" % (label - n))):
-                    labelPrev = label - n
-                    break
-            if labelPrev is not -1:
-                lastFnExp = self._getExtraPath(nameExp + "%d.xmd" % (labelPrev))
-                self.runJob("xmipp_metadata_utilities",
-                            " -i %s --set union %s -o %s " %
-                            (lastFnExp, fnExp, fnExp), numberOfMpi=1)
-        remove(fnProj)
 
     def trainNClassifiers2ClassesStep(self, thIdx, gpuId, totalGpu):
 
-        mdNumCones = emlib.MetaData(self._getExtraPath("coneCenters.doc"))
+        mdNumCones = emlib.MetaData(self._getExtraPath(self._fnConeCenterDoc))
         self.numCones = mdNumCones.size()
 
         for i in range(self.numCones):
@@ -510,53 +523,55 @@ def trainNClassifiers2ClassesStep(self, thIdx, gpuId, totalGpu):
                     copy(self.pretrainedModels.get()._getExtraPath(modelFn + '.h5'),
                          self._getExtraPath(modelFn + '.h5'))
 
-            expCheck = self._getExtraPath('projectionsExp%d.xmd' % idx)
-            if exists(expCheck):
-                if not exists(self._getExtraPath(modelFn + '.h5')):
-                    fnLabels = self._getExtraPath('labels.txt')
-                    fileLabels = open(fnLabels, "r")
-                    expSet = self._getExtraPath('projectionsExp%d.xmd' % self.numCones)
-                    if not exists(expSet):
-                        for n in range(1, self.numCones):
-                            if exists(self._getExtraPath('projectionsExp%d.xmd' % (self.numCones - n))):
-                                expSet = self._getExtraPath('projectionsExp%d.xmd' % (self.numCones - n))
-                                break
-                    newFnLabels = self._getExtraPath('labels%d.txt' % idx)
-                    newFileLabels = open(newFnLabels, "w")
-                    lines = fileLabels.readlines()
-                    for line in lines:
-                        if line == str(idx - 1) + '\n':
-                            newFileLabels.write('1\n')
-                        else:
-                            newFileLabels.write('0\n')
-                    newFileLabels.close()
-                    fileLabels.close()
-
-                    newXdim = readInfoField(self._getExtraPath(), "size",
-                                            emlib.MDL_XSIZE)
-                    fnLabels = self._getExtraPath('labels%d.txt' % idx)
-
-                    print("Training region ", idx, " in GPU ", gpuId)
-                    sys.stdout.flush()
-
-                    try:
-                        args = "%s %s %s %s %d %d %d %d " % (
-                        expSet, fnLabels, self._getExtraPath(),
-                        modelFn+'_aux', self.numEpochs, newXdim, 2, self.batchSize.get())
-                        #args += " %(GPU)s"
-                        args += " %s " % (gpuId)
-                        #args += " %s " %(int(idx % totalGpu))
-                        self.runJob("xmipp_cone_deepalign", args, numberOfMpi=1, env=self.getCondaEnv())
-                    except Exception as e:
-                        raise Exception(
-                            "ERROR: Please, if you are suffering memory problems, "
-                            "check the target resolution to work with lower dimensions.")
-
-                    moveFile(self._getExtraPath(modelFn + '_aux.h5'), self._getExtraPath(modelFn + '.h5'))
+            expCheck = self._getProjectionsExp(idx)
+            if exists(expCheck) and not exists(self._getExtraPath(modelFn + '.h5')):
+                self._coneStep(gpuId, idx, modelFn)
+
+    def _coneStep(self, gpuId, idx, modelFn):
+        fnLabels = self._getExtraPath('labels.txt')
+        fileLabels = open(fnLabels, "r")
+        expSet = self._getProjectionsExp(self.numCones)
+        if not exists(expSet):
+            for n in range(1, self.numCones):
+                if exists(self._getProjectionsExp(self.numCones - n)):
+                    expSet = self._getProjectionsExp(self.numCones - n)
+                    break
+        newFnLabels = self._getExtraPath('labels%d.txt' % idx)
+        newFileLabels = open(newFnLabels, "w")
+        lines = fileLabels.readlines()
+        for line in lines:
+            if line == str(idx - 1) + '\n':
+                newFileLabels.write('1\n')
+            else:
+                newFileLabels.write('0\n')
+        newFileLabels.close()
+        fileLabels.close()
+
+        newXdim = readInfoField(self._getExtraPath(), "size",
+                                        emlib.MDL_XSIZE)
+        fnLabels = self._getExtraPath('labels%d.txt' % idx)
+
+        print("Training region ", idx, " in GPU ", gpuId)
+        sys.stdout.flush()
+
+        try:
+            args = "%s %s %s %s %d %d %d %d " % (
+                    expSet, fnLabels, self._getExtraPath(),
+                    modelFn+'_aux', self.numEpochs, newXdim, 2, self.batchSize.get())
+                    #args += " %(GPU)s"
+            args += " %s " % (gpuId)
+                    #args += " %s " %(int(idx % totalGpu))
+            self.runJob("xmipp_cone_deepalign", args, numberOfMpi=1, env=self.getCondaEnv())
+        except Exception as e:
+            raise Exception(
+                        "ERROR: Please, if you are suffering memory problems, "
+                        "check the target resolution to work with lower dimensions.")
+
+        moveFile(self._getExtraPath(modelFn + '_aux.h5'), self._getExtraPath(modelFn + '.h5'))
 
     def predictStep(self, gpuId):
 
-        if not exists(self._getExtraPath('conePrediction.txt')):
+        if not exists(self._getConePrediction()):
             # if self.useQueueForSteps() or self.useQueue():
             #     myStr = os.environ["CUDA_VISIBLE_DEVICES"]
             # else:
@@ -566,7 +581,7 @@ def predictStep(self, gpuId):
             # print("Predict", myStr, numGPU)
             # sys.stdout.flush()
 
-            mdNumCones = emlib.MetaData(self._getExtraPath("coneCenters.doc"))
+            mdNumCones = emlib.MetaData(self._getExtraPath(self._fnConeCenterDoc))
             self.numCones = mdNumCones.size()
 
             imgsOutStk = self._getExtraPath('images_out_filtered.stk')
@@ -588,11 +603,11 @@ def predictStep(self, gpuId):
 
     def correlationCudaStep(self, thIdx, gpuId, totalGpu):
 
-        mdNumCones = emlib.MetaData(self._getExtraPath("coneCenters.doc"))
+        mdNumCones = emlib.MetaData(self._getExtraPath(self._fnConeCenterDoc))
         self.numCones = mdNumCones.size()
 
         # Cuda Correlation step - creating the metadata
-        predCones = np.loadtxt(self._getExtraPath('conePrediction.txt'))
+        predCones = np.loadtxt(self._getConePrediction())
         mdConeList = []
         numMax = int(self.numConesSelected)
         for i in range(self.numCones):
@@ -633,7 +648,7 @@ def correlationCudaStep(self, thIdx, gpuId, totalGpu):
 
                 self.runJob("xmipp_metadata_utilities", "-i %s --fill ref lineal 1 1 " % (fnProjCone), numberOfMpi=1)
 
-                fnOutCone = 'outCone%d.xmd' % (i + 1)
+                fnOutCone = self._getOutCone(i + 1)
 
                 if not exists(self._getExtraPath(fnOutCone)):
                     params = '  -i %s' % fnExpCone
@@ -645,11 +660,11 @@ def correlationCudaStep(self, thIdx, gpuId, totalGpu):
 
     def correlationSignificantStep(self):
 
-        mdNumCones = emlib.MetaData(self._getExtraPath("coneCenters.doc"))
+        mdNumCones = emlib.MetaData(self._getExtraPath(self._fnConeCenterDoc))
         self.numCones = mdNumCones.size()
 
         # Cuda Correlation step - creating the metadata
-        predCones = np.loadtxt(self._getExtraPath('conePrediction.txt'))
+        predCones = np.loadtxt(self._getConePrediction())
         mdConeList = []
         numMax = int(self.numConesSelected)
         for i in range(self.numCones):
@@ -674,7 +689,7 @@ def correlationSignificantStep(self):
                 mdConeList[i].write(fnExpCone)
 
                 fnProjCone = self._getExtraPath('projectionsCudaCorr%d.xmd' % (i + 1))
-                fnOutCone = 'outCone%d.xmd' % (i + 1)
+                fnOutCone = self._getOutCone(i + 1)
 
                 if not exists(self._getExtraPath(fnOutCone)):
                     # Correlation step - calling significant program
@@ -689,7 +704,7 @@ def correlationSignificantStep(self):
 
     def createOutputMetadataStep(self):
 
-        mdNumCones = emlib.MetaData(self._getExtraPath("coneCenters.doc"))
+        mdNumCones = emlib.MetaData(self._getExtraPath(self._fnConeCenterDoc))
         self.numCones = mdNumCones.size()
         numMax = int(self.numConesSelected)
         mdIn = emlib.MetaData(self.imgsFn)
@@ -703,27 +718,11 @@ def createOutputMetadataStep(self):
         fnFinal = self._getExtraPath('outConesParticles.xmd')
         for i in range(self.numCones):
 
-            fnOutCone = 'outCone%d.xmd' % (i + 1)
+            fnOutCone = self._getOutCone(i + 1)
 
             if exists(self._getExtraPath(fnOutCone)):             
-
-                if numMax == 1:
-                    if not exists(fnFinal):
-                        copy(self._getExtraPath(fnOutCone), fnFinal)
-                    else:
-                        params = ' -i %s --set union %s -o %s' % (fnFinal, self._getExtraPath(fnOutCone),
-                                                                  fnFinal)
-                        self.runJob("xmipp_metadata_utilities", params,
-                                    numberOfMpi=1)
-                else:
-                    mdCones.append(emlib.MetaData(self._getExtraPath(fnOutCone)))
-                    coneFns.append(mdCones[i].getColumnValues(emlib.MDL_IMAGE))
-                    shiftX.append(mdCones[i].getColumnValues(emlib.MDL_SHIFT_X))
-                    shiftY.append(mdCones[i].getColumnValues(emlib.MDL_SHIFT_Y))
-                    coneCCs.append(mdCones[i].getColumnValues(emlib.MDL_MAXCC))
-
-            else:
-                if numMax > 1:
+                self._updateCone(numMax, coneFns, coneCCs, mdCones, shiftX, shiftY, fnFinal, i, fnOutCone)
+            elif numMax > 1:
                     mdCones.append(None)
                     coneFns.append([])
                     coneCCs.append([])
@@ -731,29 +730,7 @@ def createOutputMetadataStep(self):
                     shiftY.append([])
 
         if numMax > 1:
-            mdFinal = emlib.MetaData()
-            row = md.Row()
-            for myFn in allInFns:
-                myCCs = []
-                myCones = []
-                myPos = []
-                for n in range(self.numCones):
-                    if myFn in coneFns[n]:
-                        pos = coneFns[n].index(myFn)
-                        myPos.append(pos)
-                        if abs(shiftX[n][pos])<30 and abs(shiftY[n][pos])<30:
-                            myCCs.append(coneCCs[n][pos])
-                        else:
-                            myCCs.append(0)
-                        myCones.append(n + 1)
-                if len(myPos) > 0:
-                    if max(myCCs)==0:
-                        continue
-                    coneMax = myCones[myCCs.index(max(myCCs))]
-                    objId = myPos[myCCs.index(max(myCCs))] + 1
-                    row.readFromMd(mdCones[coneMax - 1], objId)
-                    row.addToMd(mdFinal)
-            mdFinal.write(fnFinal)
+            self._writeMDFinal(allInFns, coneFns, coneCCs, mdCones, shiftX, shiftY, fnFinal)
         else:
             mdCones = emlib.MetaData(fnFinal)
             mdCones.removeObjects(emlib.MDValueGT(emlib.MDL_SHIFT_X, 30.))
@@ -762,6 +739,47 @@ def createOutputMetadataStep(self):
             mdCones.removeObjects(emlib.MDValueLT(emlib.MDL_SHIFT_Y, -30.))
             mdCones.write(fnFinal)
 
+    def _updateCone(self, numMax, coneFns, coneCCs, mdCones, shiftX, shiftY, fnFinal, i, fnOutCone):
+        if numMax == 1:
+            if not exists(fnFinal):
+                copy(self._getExtraPath(fnOutCone), fnFinal)
+            else:
+                params = ' -i %s --set union %s -o %s' % (fnFinal, self._getExtraPath(fnOutCone),
+                                                                  fnFinal)
+                self.runJob("xmipp_metadata_utilities", params,
+                                    numberOfMpi=1)
+        else:
+            mdCones.append(emlib.MetaData(self._getExtraPath(fnOutCone)))
+            coneFns.append(mdCones[i].getColumnValues(emlib.MDL_IMAGE))
+            shiftX.append(mdCones[i].getColumnValues(emlib.MDL_SHIFT_X))
+            shiftY.append(mdCones[i].getColumnValues(emlib.MDL_SHIFT_Y))
+            coneCCs.append(mdCones[i].getColumnValues(emlib.MDL_MAXCC))
+
+    def _writeMDFinal(self, allInFns, coneFns, coneCCs, mdCones, shiftX, shiftY, fnFinal):
+        mdFinal = emlib.MetaData()
+        row = md.Row()
+        for myFn in allInFns:
+            myCCs = []
+            myCones = []
+            myPos = []
+            for n in range(self.numCones):
+                if myFn in coneFns[n]:
+                    pos = coneFns[n].index(myFn)
+                    myPos.append(pos)
+                    if abs(shiftX[n][pos])<30 and abs(shiftY[n][pos])<30:
+                        myCCs.append(coneCCs[n][pos])
+                    else:
+                        myCCs.append(0)
+                    myCones.append(n + 1)
+            if len(myPos) > 0:
+                if max(myCCs)==0:
+                    continue
+                coneMax = myCones[myCCs.index(max(myCCs))]
+                objId = myPos[myCCs.index(max(myCCs))] + 1
+                row.readFromMd(mdCones[coneMax - 1], objId)
+                row.addToMd(mdFinal)
+        mdFinal.write(fnFinal)
+
     def createOutputStep(self):
 
         cleanPattern(self._getExtraPath('*.stk'))
diff --git a/xmipp3/protocols/protocol_kmeans_clustering.py b/xmipp3/protocols/protocol_kmeans_clustering.py
new file mode 100644
index 000000000..344643c29
--- /dev/null
+++ b/xmipp3/protocols/protocol_kmeans_clustering.py
@@ -0,0 +1,122 @@
+# # *************************************************
+# # * This protocol will be offered in future releases (more testing is needed)
+# # *************************************************
+
+# # *****************************************************************************
+# # *
+# # * Authors:     David Herreros Calero     (dherreros@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 os
+# from sklearn.cluster import KMeans
+# import numpy as np
+#
+# from pyworkflow.object import Integer
+# import pyworkflow.protocol.params as param
+# from pyworkflow import VERSION_2_0
+# import pwem.emlib.metadata as md
+# from pyworkflow.utils.properties import Message
+# from pwem.protocols import ProtClassify2D, pwobj
+#
+# from xmipp3.convert import getXmippAttribute
+#
+#
+# class XmippProtKmeansSPH(ProtClassify2D):
+#     """ Reduce the number of deformation coefficients based on k-means vectors """
+#
+#     _label = 'kmeans clustering'
+#     _lastUpdateVersion = VERSION_2_0
+#
+#     def __init__(self, **args):
+#         ProtClassify2D.__init__(self, **args)
+#         self.stepsExecutionMode = param.STEPS_PARALLEL
+#
+#     # --------------------------- DEFINE param functions -----------------------
+#     def _defineParams(self, form):
+#         form.addSection(label=Message.LABEL_INPUT)
+#         form.addParam('inputParts', param.MultiPointerParam, pointerClass='SetOfParticles',
+#                       label="Input particles",
+#                       important=True)
+#         form.addParam('column', param.EnumParam, choices=['SPH', 'NMA'],
+#                       label='Column to cluster', default=0,
+#                       help='Selecto from the type of coefficients to be used for the '
+#                            'clustering.')
+#
+#     # --------------------------- INSERT steps functions -----------------------
+#     def _insertAllSteps(self):
+#         self._insertFunctionStep('findClustersStep')
+#         self._insertFunctionStep('createOutputStep')
+#
+#     def findClustersStep(self):
+#         self.kDict = []
+#         coeffs = []
+#         for inputParts in self.inputParts:
+#             particles = inputParts.get()
+#             colData = self.column.get()
+#             if colData == 0:
+#                 mdLabel = md.MDL_SPH_COEFFICIENTS
+#             elif colData == 1:
+#                 mdLabel = md.MDL_NMA
+#             for particle in particles.iterItems():
+#                 coeffs.append(np.fromstring(getXmippAttribute(particle, mdLabel).get(), sep=','))
+#         maxClusters = 101 if len(coeffs) > 101 else len(coeffs) + 1
+#         rmse = np.zeros((maxClusters - 1))
+#         coeffs = np.asarray(coeffs)[:, :-8]
+#         for nClusters in range(1, maxClusters):
+#             kmeans = KMeans(n_clusters=nClusters).fit(coeffs)
+#             rmse[nClusters-1] = np.sqrt(kmeans.inertia_ / len(coeffs))
+#         p1 = np.array((1, rmse[0]))
+#         p2 = np.array((maxClusters - 1, rmse[-1]))
+#         d = np.zeros((maxClusters - 1))
+#         for nClusters in range(2, maxClusters - 1):
+#             p3 = np.array((nClusters, rmse[nClusters - 1]))
+#             d[nClusters-1] = np.linalg.norm(np.cross(p2 - p1, p1 - p3)) / np.linalg.norm(p2 - p1)
+#         nClusters = np.argmax(d) + 1
+#         self.kmeans = KMeans(n_clusters=nClusters).fit(np.asarray(coeffs))
+#
+#     def createOutputStep(self):
+#         classes2DSet = self._createSetOfClasses2D(self.inputParts[0].get())
+#         classes2DSet.classifyItems(updateItemCallback=self._updateParticle,
+#                                    updateClassCallback=self._updateClass,
+#                                    itemDataIterator=iter(self.kmeans.labels_))
+#         result = {'outputClasses': classes2DSet}
+#         classes2DSet.L1 = Integer(self.inputParts[0].get().L1)
+#         classes2DSet.L2 = Integer(self.inputParts[0].get().L2)
+#         classes2DSet.Rmax = Integer(self.inputParts[0].get().Rmax)
+#         self._defineOutputs(**result)
+#         for inputParts in self.inputParts:
+#             self._defineSourceRelation(inputParts, classes2DSet)
+#
+#     # ------------------------------- UTILS functions -------------------------------
+#     def _updateParticle(self, item, idc):
+#         item.setClassId(idc)
+#
+#     def _updateClass(self, item):
+#         coeff = self.kmeans.cluster_centers_[item.getObjId()-1]
+#         coeff = ','.join(['%f' % num for num in coeff])
+#         representative = pwobj.CsvList()
+#         representative.set(coeff)
+#         # setattr(representative, '_classID', pwobj.Integer(item.getObjId()))
+#         # setattr(representative, '_objDoStore', True)
+#         item.setRepresentative(representative)
diff --git a/xmipp3/protocols/protocol_metaprotocol_golden_highres.py b/xmipp3/protocols/protocol_metaprotocol_golden_highres.py
index 3f90f0b50..394826c3e 100644
--- a/xmipp3/protocols/protocol_metaprotocol_golden_highres.py
+++ b/xmipp3/protocols/protocol_metaprotocol_golden_highres.py
@@ -247,7 +247,7 @@ def monitorStep(self):
             newHighRes.inputVolumes.set(previousProtVol)
             newHighRes.inputVolumes.setExtended(namePreviousVol)
 
-            project.launchProtocol(newHighRes)
+            project.scheduleProtocol(newHighRes)
 
             # Next schedule will be after this one
             self._runPrerequisites.append(newHighRes.getObjId())
diff --git a/xmipp3/protocols/protocol_movie_correlation.py b/xmipp3/protocols/protocol_movie_correlation.py
index 6d6caa6a8..139d4e6ea 100644
--- a/xmipp3/protocols/protocol_movie_correlation.py
+++ b/xmipp3/protocols/protocol_movie_correlation.py
@@ -352,8 +352,8 @@ def _setControlPoints(self):
         x, y, frames = self.inputMovies.get().getDim()
         Ts = self.inputMovies.get().getSamplingRate()
         # one control point each 1000 A
-        self.controlPointX.set(max([int(x / Ts) / 1000 + 2, 3]))
-        self.controlPointY.set(max([int(y / Ts) / 1000 + 2, 3]))
+        self.controlPointX.set(max([int(x * Ts) / 1000 + 2, 3]))
+        self.controlPointY.set(max([int(y * Ts) / 1000 + 2, 3]))
         self.controlPointT.set(max([ceil(frames/7.) + 2, 3]))
 
     def _getMovieShifts(self, movie):
diff --git a/xmipp3/protocols/protocol_pdb_deform_zernike3d.py b/xmipp3/protocols/protocol_pdb_deform_zernike3d.py
new file mode 100644
index 000000000..f4896d539
--- /dev/null
+++ b/xmipp3/protocols/protocol_pdb_deform_zernike3d.py
@@ -0,0 +1,64 @@
+
+# **************************************************************************
+# *
+# * Authors:     David Herreros Calero (dherreros@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 pwem.protocols import ProtAnalysis3D
+from pwem.objects import AtomStruct
+import pyworkflow.protocol.params as params
+import pyworkflow.utils as pwutils
+
+
+class XmippProtDeformPDBZernike3D(ProtAnalysis3D):
+    """ Protocol for PDB deformation based on Zernike3D basis. """
+    _label = 'deform pdb - Zernike3D'
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('inputPDB', params.PointerParam, label="Input PDB",
+                      pointerClass='AtomStruct', important=True,
+                      help='Select a PDB to be deformed.')
+        form.addParam('inputCoeff', params.PathParam, label="Input Coefficients",
+                      important=True,
+                      help='Specify a path to the deformation coefficients file.')
+
+    # --------------------------- INSERT steps functions -----------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep("deformStep")
+        self._insertFunctionStep("createOutputStep")
+
+    # --------------------------- STEPS functions ------------------------------
+    def deformStep(self):
+        outFile = pwutils.removeBaseExt(self.inputPDB.get().getFileName()) + '_deformed.pdb'
+        params = ' --pdb %s --clnm %s -o %s' % \
+                 (self.inputPDB.get().getFileName(), self.inputCoeff.get(), self._getExtraPath(outFile))
+        self.runJob("xmipp_pdb_sph_deform", params)
+
+    def createOutputStep(self):
+        outFile = pwutils.removeBaseExt(self.inputPDB.get().getFileName()) + '_deformed.pdb'
+        pdb = AtomStruct(self._getExtraPath(outFile))
+        self._defineOutputs(outputPDB=pdb)
+        self._defineSourceRelation(self.inputPDB, pdb)
diff --git a/xmipp3/protocols/protocol_screen_deepConsensus.py b/xmipp3/protocols/protocol_screen_deepConsensus.py
index a502e2623..c95cafb2e 100644
--- a/xmipp3/protocols/protocol_screen_deepConsensus.py
+++ b/xmipp3/protocols/protocol_screen_deepConsensus.py
@@ -28,7 +28,7 @@
 """
 Deep Consensus picking protocol
 """
-import os, sys
+import os, sys, time
 from glob import glob
 import six
 import json, shutil, pickle
@@ -78,10 +78,14 @@ class XmippProtScreenDeepConsensus(ProtParticlePicking, XmippProtocol):
     _conda_env = 'xmipp_DLTK_v0.3'
     _stepsCheckSecs = 5              # time in seconds to check the steps
 
+    USING_INPUT_COORDS = False
     CONSENSUS_COOR_PATH_TEMPLATE="consensus_coords_%s"
     CONSENSUS_PARTS_PATH_TEMPLATE="consensus_parts_%s"
     PRE_PROC_MICs_PATH="preProcMics"
 
+    PARTICLES_TEMPLATE = "particles{}.xmd"
+    NET_TEMPLATE = "nnetData{}"
+
     ADD_DATA_TRAIN_TYPES = ["None", "Precompiled", "Custom"]
     ADD_DATA_TRAIN_NONE = 0
     ADD_DATA_TRAIN_PRECOMP = 1
@@ -450,6 +454,9 @@ def _doContinue(self):
 
     def justPredict(self):
       return self.skipTraining.get()==True
+    
+    def _useNewModel(self):
+      return self.modelInitialization.get() == self.ADD_MODEL_TRAIN_NEW
 
     def _usePretrainedModel(self):
         return self.modelInitialization.get()== self.ADD_MODEL_TRAIN_PRETRAIN
@@ -465,81 +472,52 @@ def _insertAllSteps(self):
 
     def _stepsCheck(self):
         '''Checks if new steps can be executed'''
-        newSteps = []
-        if not self.ENDED:
-          # Functions streamed. Input is processed as soon as it is generated
-          if not self.checkIfPrevRunIsCompatible("mics_"):
-            if len(self.readyToPreprocessMics(shared = False)) > 0 and not self.PREPROCESSING:
-              self.PREPROCESSING = True
-              self.lastDeps = [self._insertFunctionStep("preprocessMicsStep", prerequisites=self.initDeps)]
-
-          if len(self.readyToExtractMicFns('OR')) >= self.extractingBatch.get() and not self.EXTRACTING['OR']:
-            self.EXTRACTING['OR'] = True
-            newSteps += [self.insertCaculateConsensusSteps('OR', prerequisites=self.initDeps)]
-            newSteps += self.insertExtractPartSteps('OR', prerequisites=newSteps)
-
-
-          trainedParams = self.loadTrainedParams()
-          toTrainSize = self.toTrainDataSize.get() if self.toTrainDataSize.get() != -1 else 1e10
-          if not self.skipTraining.get() and trainedParams['posParticlesTrained'] < toTrainSize:
-            if len(self.readyToExtractMicFns('NOISE')) >= self.extractingBatch.get() and not self.EXTRACTING['NOISE']:
-              self.EXTRACTING['NOISE'] = True
-              depNoise = self._insertFunctionStep('pickNoise', prerequisites=self.initDeps)
-              depsNoise = self.insertExtractPartSteps('NOISE', prerequisites=[depNoise])
-              newSteps += depsNoise
-
-            if len(self.readyToExtractMicFns('AND')) >= self.extractingBatch.get() and not self.EXTRACTING['AND']:
-              self.EXTRACTING['AND'] = True
-              depsAnd = self.insertCaculateConsensusSteps('AND', prerequisites=self.initDeps)
-              depsAnd = self.insertExtractPartSteps('AND', prerequisites=[depsAnd])
-              newSteps += depsAnd
-
-            trainedParams = self.loadTrainedParams()
-            if self.addTrainingData.get() == self.ADD_DATA_TRAIN_CUST and\
-                    self.trainingDataType == self.ADD_DATA_TRAIN_CUSTOM_OPT_COORS:
-              if self.trainTrueSetOfCoords.get() is not None and \
-                      len(self.readyToExtractMicFns('ADDITIONAL_COORDS_TRUE')) >= self.extractingBatch.get():
-                self.TO_EXTRACT_MICFNS['ADDITIONAL_COORDS_TRUE'] = self.readyToExtractMicFns('ADDITIONAL_COORDS_TRUE')
-                newSteps += self.insertExtractPartSteps('ADDITIONAL_COORDS_TRUE', prerequisites=self.initDeps)
-              if self.trainFalseSetOfCoords.get() is not None and \
-                      len(self.readyToExtractMicFns('ADDITIONAL_COORDS_FALSE')) >= self.extractingBatch.get():
-                self.TO_EXTRACT_MICFNS['ADDITIONAL_COORDS_FALSE'] = self.readyToExtractMicFns('ADDITIONAL_COORDS_FALSE')
-                newSteps += self.insertExtractPartSteps('ADDITIONAL_COORDS_FALSE', prerequisites=self.initDeps)
-
-            if self.cnnFree():
-              self.TO_TRAIN_MICFNS = self.readyToTrainMicFns()
-              if len(self.TO_TRAIN_MICFNS) >= self.trainingBatch.get():
-                self.TRAINING = True
-                self.curTrainedParams = trainedParams
-                self.depsTrain = [self._insertFunctionStep('trainCNN', self.TO_TRAIN_MICFNS, prerequisites=self.initDeps)]
-                self.depsTrain = [self._insertFunctionStep('endTrainingStep', prerequisites=self.depsTrain)]
-                newSteps += self.depsTrain
-
-          elif trainedParams['posParticlesTrained'] >= toTrainSize and trainedParams['trainingPass'] != '':
-            lastTrainingPass = trainedParams['trainingPass']
-            self.retrievePreviousPassModel('', lastTrainingPass)
-            trainedParams['trainingPass'] = ''
-            self.saveTrainedParams(trainedParams)
-
-          if self.networkReadyToPredict() and self.cnnFree() and self.predictionsOn():
-              self.PREDICTING = True
-              depPredict = self._insertFunctionStep('predictCNN', prerequisites= newSteps)
-              newSteps += [self._insertFunctionStep('endPredictingStep', prerequisites=[depPredict])]
-              newSteps += [self._insertFunctionStep('createOutputStep', prerequisites=[depPredict])]
-
-          trainedParams = self.loadTrainedParams()
-          if self.allFree() and not self.LAST_ROUND and self.checkIfParentsFinished():
-            protLast = self._steps[self.lastStep - 1]
-            protLast.addPrerequisites(*newSteps)
-            protLast.setStatus(STATUS_NEW)
-            self.updateSteps()
-
-          protEnd = self._steps[self.endStep-1]
-          protEnd.addPrerequisites(*newSteps)
-          if self.LAST_ROUND and self.allFree():
-            protEnd.setStatus(STATUS_NEW)
+        self.newSteps = []
+        if self.ENDED:
+          return
+        # Preprocessing
+        if len(self.readyToPreprocessMics(shared = False)) > 0 and not self.PREPROCESSING:
+            self.PREPROCESSING = True
+            self.lastDeps = [self._insertFunctionStep("preprocessMicsStep", prerequisites=self.initDeps)]
+
+        #Particle extraction OR (for predictions)
+        if len(self.readyToExtractMicFns('OR')) >= self.extractingBatch.get() and not self.EXTRACTING['OR']:
+          self.EXTRACTING['OR'] = True
+          self.newSteps += self.insertCaculateConsensusSteps('OR', prerequisites=self.initDeps)
+          self.newSteps += self.insertExtractPartSteps('OR', prerequisites=self.newSteps)
+
+        #Particle extraction for training and training
+        trainedParams = self.loadTrainedParams()
+        toTrainSize = self.toTrainDataSize.get() if self.toTrainDataSize.get() != -1 else 1e10
+        if not self.skipTraining.get() and trainedParams['posParticlesTrained'] < toTrainSize:
+          self.doTraining()
+
+        elif trainedParams['posParticlesTrained'] >= toTrainSize and trainedParams['trainingPass'] != '':
+          lastTrainingPass = trainedParams['trainingPass']
+          self.retrievePreviousPassModel('', lastTrainingPass)
+          trainedParams['trainingPass'] = ''
+          self.saveTrainedParams(trainedParams)
+
+        #Prediction
+        if self.networkReadyToPredict() and self.cnnFree() and self.predictionsOn():
+            self.PREDICTING = True
+            depPredict = self._insertFunctionStep('predictCNN', prerequisites= self.newSteps)
+            self.newSteps += [self._insertFunctionStep('endPredictingStep', prerequisites=[depPredict])]
+            self.newSteps += [self._insertFunctionStep('createOutputStep', prerequisites=[depPredict])]
+
+        #Last round with batch size == 1 to include all input
+        if self.allFree() and not self.LAST_ROUND and self.checkIfParentsFinished():
+          protLast = self._steps[self.lastStep - 1]
+          protLast.addPrerequisites(*self.newSteps)
+          protLast.setStatus(STATUS_NEW)
           self.updateSteps()
-          sys.stdout.flush()
+
+        protEnd = self._steps[self.endStep-1]
+        protEnd.addPrerequisites(*self.newSteps)
+        #Ending the protocol when everything is done
+        if self.LAST_ROUND and self.allFree():
+          protEnd.setStatus(STATUS_NEW)
+        self.updateSteps()
 
 
     def endPredictingStep(self):
@@ -593,39 +571,11 @@ def initializeStep(self):
             writeSetOfParticles(self.retrieveTrainSets(),
                                 self._getTmpPath("addNegTrainParticles.xmd"))
 
-        if self.checkIfPrevRunIsCompatible(""):
-            for mode in ["OR", "AND", "NOISE"]:
-              print("copying particles %s"%(mode))
-              newDataPath= "particles_%s.xmd" % mode
-              createLink( self.continueRun.get()._getExtraPath(newDataPath),
-                        self._getExtraPath(newDataPath))
-              newDataPath= "parts_%s" % mode
-              createLink( self.continueRun.get()._getExtraPath(newDataPath),
-                        self._getExtraPath(newDataPath))
-            print("copying OR coordinates")
-            newDataPath= (self.CONSENSUS_COOR_PATH_TEMPLATE % "OR")
-            createLink( self.continueRun.get()._getExtraPath(newDataPath),
-                        self._getExtraPath(newDataPath))
 
-        else:
-          if self.checkIfPrevRunIsCompatible( "mics_"):
-            sourcePath= self.continueRun.get()._getExtraPath(self.PRE_PROC_MICs_PATH )
-            print("copying mics from %s"%(sourcePath))
-            createLink( sourcePath, self._getExtraPath(self.PRE_PROC_MICs_PATH ))
-
-          else:
-            makePath(self._getExtraPath(self.PRE_PROC_MICs_PATH))
-
-          if self.checkIfPrevRunIsCompatible( "coords_"):
-            for mode in ["OR", "AND", "NOISE"]:
-              print("copying coordinates %s"%(mode))
-              newDataPath= (self.CONSENSUS_COOR_PATH_TEMPLATE % mode)
-              createLink( self.continueRun.get()._getExtraPath(newDataPath),
-                        self._getExtraPath(newDataPath))
-          else:
-              for mode in ["AND", "OR", "NOISE"]:
-                consensusCoordsPath = self.CONSENSUS_COOR_PATH_TEMPLATE % mode
-                makePath(self._getExtraPath(consensusCoordsPath))
+        makePath(self._getExtraPath(self.PRE_PROC_MICs_PATH))
+        for mode in ["AND", "OR", "NOISE"]:
+            consensusCoordsPath = self.CONSENSUS_COOR_PATH_TEMPLATE % mode
+            makePath(self._getExtraPath(consensusCoordsPath))
 
         if self._doContinue():
           if self.skipTraining.get():
@@ -635,12 +585,20 @@ def initializeStep(self):
           self.retrievePreviousRunModel(self.continueRun.get(), trPass)
           self.uploadTrainedParam('trainingPass', trPass)
 
+        if self._usePretrainedModel() and self.skipTraining.get():
+          trPass=''
+          self.uploadTrainedParam('trainingPass', trPass)
+        
+        if self._useNewModel() and self.skipTraining.get():
+          self.skipTraining.set(False)
+
         preprocessParamsFname= self._getExtraPath("preprocess_params.json")
         preprocParams= self.getPreProcParamsFromForm()
         with open(preprocessParamsFname, "w") as f:
           json.dump(preprocParams, f)
 
         #Initializing outputs
+        self.inSamplingRate = self._getInputMicrographs().getSamplingRate()
         self.preCorrectedParSet, self.preCoordSet = [], []
 
 
@@ -648,7 +606,6 @@ def preprocessMicsStep(self):
         '''Step which preprocesses the input micrographs'''
         micIds = self.getMicsIds(filterOutNoCoords=False)
         if len(micIds) > 0:
-          samplingRate = self._getInputMicrographs().getSamplingRate()
           mics_ = self._getInputMicrographs()
           micsFnameSet = {mics_[micId].getMicName(): mics_[micId].getFileName() for micId in micIds
                           if mics_[micId] is not None}  # to skip failed mics
@@ -679,7 +636,7 @@ def preprocessMicsStep(self):
           with open(inputsFname, "w") as f:
             f.write(preproMicsContent)
           downFactor = self._getDownFactor()
-          args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, samplingRate, downFactor, ouputDir, nThrs)
+          args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, self.inSamplingRate, downFactor, ouputDir, nThrs)
           if not self.skipInvert.get():
             args+=" --invert_contrast"
 
@@ -697,7 +654,12 @@ def insertCaculateConsensusSteps(self, mode, prerequisites):
           makePath(outCoordsDataPath)
         newDep = self._insertFunctionStep('calculateCoorConsensusStep', outCoordsDataPath, mode, prerequisites=prerequisites)
         newDep = self._insertFunctionStep('loadCoords', outCoordsDataPath, mode, prerequisites=[newDep])
-        return newDep
+        return [newDep]
+    
+    def waitFreeInputCoords(self):
+      while self.USING_INPUT_COORDS:
+        time.sleep(1)
+      self.USING_INPUT_COORDS = True
 
     def calculateCoorConsensusStep(self, outCoordsDataPath, mode):
       '''Calculates the consensus coordinates from micrographs whose particles haven't been extracted yet in "mode"'''
@@ -707,24 +669,10 @@ def calculateCoorConsensusStep(self, outCoordsDataPath, mode):
         for micFn in self.getExtractedMicFns(mode):
           extractedSetOfCoordsFns.append(pwutils.path.replaceBaseExt(micFn,'pos'))
 
-        Tm = []
-        for coordinatesP in self.inputCoordinates:
-          mics = coordinatesP.get().getMicrographs()
-          Tm.append(mics.getSamplingRate())
-        nCoordsSets= len(Tm)
-
-        inputCoordsFnames={}
-        for coord_num, coordinatesP in enumerate(self.inputCoordinates):
-            tmpPosDir= self._getTmpPath("input_coords_%d_%s"%(coord_num, mode))
-            if not os.path.exists(tmpPosDir):
-              makePath(tmpPosDir)
-            writeSetOfCoordinates(tmpPosDir, coordinatesP.get(), scale=float(Tm[coord_num])/float(Tm[0]))
-            for posFname in os.listdir(tmpPosDir):
-                baseName, extension=os.path.splitext(os.path.basename(posFname))
-                if extension==".pos" and not posFname in extractedSetOfCoordsFns:
-                  if baseName not in inputCoordsFnames:
-                      inputCoordsFnames[baseName]=["None"]*nCoordsSets
-                  inputCoordsFnames[baseName][coord_num]= os.path.join(tmpPosDir, posFname)
+        self.waitFreeInputCoords()
+        inputCoordsFnames = self.getInpCoordsFns(mode, extractedSetOfCoordsFns)
+        self.USING_INPUT_COORDS = False
+          
         inputFileHeader="#pos_i\n"
         inputFileStr=inputFileHeader
         for baseName in inputCoordsFnames:
@@ -744,10 +692,7 @@ def calculateCoorConsensusStep(self, outCoordsDataPath, mode):
         self.TO_EXTRACT_MICFNS[mode] = self.readyToExtractMicFns(mode)
 
     def pickNoise(self):
-      '''Find noise coordinates from micrographs in order to use them as negatives in the training process'''
-      if self.checkIfPrevRunIsCompatible("coords_"):
-        print("using previous round noise particles")
-      else:
+        '''Find noise coordinates from micrographs in order to use them as negatives in the training process'''
         orPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "OR")
         outputPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "NOISE")
         if not "OR" in self.coordinatesDict:  # fill self.coordinatesDict['OR']
@@ -790,8 +735,7 @@ def pickNoise(self):
           args = (" -i %(toPickMicsDir)s -c %(toPickCoorsDir)s -o %(outputPosDir)s -s %(boxSize)s " +
                   "-n %(nToPick)s -t %(nThrs)s") % argsDict
 
-          if not self.checkIfPrevRunIsCompatible("coords_"):
-            self.runJob('xmipp_pick_noise', args, numberOfMpi=1)
+          self.runJob('xmipp_pick_noise', args, numberOfMpi=1)
 
           self.loadCoords(outputPosDir, 'NOISE', micSet=micNoiseSet)
           self.TO_EXTRACT_MICFNS['NOISE'] = toPickNoiseFns
@@ -815,9 +759,12 @@ def loadCoords(self, posCoorsPath, mode, micSet=[]):
           sqliteName = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode) + ".sqlite"
           if os.path.isfile(self._getExtraPath(sqliteName)):
             cleanPath(self._getExtraPath(sqliteName))
+          
+          self.waitFreeInputCoords()
           totalSetOfCoordinates = readSetOfCoordsFromPosFnames(posCoorsPath,
                                                                setOfInputCoords=self.inputCoordinates[0].get(),
                                                                sqliteOutName=sqliteName, write=True)
+          self.USING_INPUT_COORDS=False
           print("Coordinates %s size: %d" % (mode, totalSetOfCoordinates.getSize()))
           assert totalSetOfCoordinates.getSize() > MIN_NUM_CONSENSUS_COORDS, \
             ("Error, the consensus (%s) of your input coordinates was too small (%s). " +
@@ -828,11 +775,10 @@ def loadCoords(self, posCoorsPath, mode, micSet=[]):
 
     def insertExtractPartSteps(self, mode, prerequisites):
         '''Inserts the steps necessary for extracting the particles from the micrographs'''
-        newSteps = []
-        if not self.checkIfPrevRunIsCompatible(""):
-          newSteps.append(self._insertFunctionStep("extractParticles", mode, prerequisites= prerequisites))
-          newSteps.append(self._insertFunctionStep("joinSetOfParticlesStep", mode, prerequisites= newSteps))
-        return newSteps
+        self.newSteps = []
+        self.newSteps.append(self._insertFunctionStep("extractParticles", mode, prerequisites= prerequisites))
+        self.newSteps.append(self._insertFunctionStep("joinSetOfParticlesStep", mode, prerequisites= self.newSteps))
+        return self.newSteps
 
     def extractParticles(self, mode):
         '''Extract the particles from a set of micrographs with their corresponding coordinates'''
@@ -841,39 +787,41 @@ def extractParticles(self, mode):
         preprocMicsPath = self._getExtraPath(self.PRE_PROC_MICs_PATH)
         toExtractMicFns = self.TO_EXTRACT_MICFNS[mode]
         print('To extract in mode {}: {}'.format(mode, toExtractMicFns))
-        if len(toExtractMicFns) > 0:
-          for micFname in toExtractMicFns:
-            micFnameBase= pwutils.removeExt(micFname)
-            micFname= os.path.join(preprocMicsPath, micFname)
-            micsFnameSet[micFnameBase]= micFname
-          extractCoordsContent="#mics coords\n"
-
-          if mode.startswith("ADDITIONAL_COORDS"):
-            if not os.path.exists(posDir):
-              os.mkdir(posDir)
-            if mode.endswith("TRUE"):
-              coordSet= self.trainTrueSetOfCoords.get()
-            elif mode.endswith("FALSE"):
-              coordSet= self.trainFalseSetOfCoords.get()
-            writeSetOfCoordinates(posDir, coordSet)
-
-          for posFname in os.listdir(posDir):
-            posNameBase=  pwutils.removeExt(posFname)
-            posFname= os.path.join(posDir, posFname)
-            if posNameBase in micsFnameSet:
-              extractCoordsContent+= "%s particles@%s\n"%(micsFnameSet[posNameBase], posFname)
-
-          inputsFname= self._getTmpPath("extractParticles_inputs_%s.txt"%mode)
-          outputDir= self._getConsensusParticlesDir(mode)
-          if not os.path.exists(outputDir):
-            makePath(outputDir)
-          nThrs= self.numberOfThreads.get()
-          with open(inputsFname, "w") as f:
-            f.write(extractCoordsContent)
-          downFactor= self._getDownFactor()
-          args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, DEEP_PARTICLE_SIZE, downFactor, outputDir, nThrs)
-
-          self.runJob('xmipp_extract_particles', args, numberOfMpi=1)
+        if len(toExtractMicFns) <= 0:
+          return
+
+        for micFname in toExtractMicFns:
+          micFnameBase= pwutils.removeExt(micFname)
+          micFname= os.path.join(preprocMicsPath, micFname)
+          micsFnameSet[micFnameBase]= micFname
+        extractCoordsContent="#mics coords\n"
+
+        if mode.startswith("ADDITIONAL_COORDS"):
+          if not os.path.exists(posDir):
+            os.mkdir(posDir)
+          if mode.endswith("TRUE"):
+            coordSet= self.trainTrueSetOfCoords.get()
+          elif mode.endswith("FALSE"):
+            coordSet= self.trainFalseSetOfCoords.get()
+          writeSetOfCoordinates(posDir, coordSet)
+
+        for posFname in os.listdir(posDir):
+          posNameBase=  pwutils.removeExt(posFname)
+          posFname= os.path.join(posDir, posFname)
+          if posNameBase in micsFnameSet:
+            extractCoordsContent+= "%s particles@%s\n"%(micsFnameSet[posNameBase], posFname)
+
+        inputsFname= self._getTmpPath("extractParticles_inputs_%s.txt"%mode)
+        outputDir= self._getConsensusParticlesDir(mode)
+        if not os.path.exists(outputDir):
+          makePath(outputDir)
+        nThrs= self.numberOfThreads.get()
+        with open(inputsFname, "w") as f:
+          f.write(extractCoordsContent)
+        downFactor= self._getDownFactor()
+        args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, DEEP_PARTICLE_SIZE, downFactor, outputDir, nThrs)
+
+        self.runJob('xmipp_extract_particles', args, numberOfMpi=1)
 
     def joinSetOfParticlesStep(self, mode, micFns='', trainingPass='', clean=False):
         '''Stores the particles extracted from a set of micrographs in a images.xmd metadata file'''
@@ -905,6 +853,38 @@ def joinSetOfParticlesStep(self, mode, micFns='', trainingPass='', clean=False):
           imgsXmd.append(fnImages)
         self.EXTRACTING[mode] = False
 
+    def doTraining(self):
+      '''Prepares the positive (AND) and negative (NOISE) coordinates for the training and executes it'''
+      if len(self.readyToExtractMicFns('NOISE')) >= self.extractingBatch.get() and not self.EXTRACTING['NOISE']:
+        self.EXTRACTING['NOISE'] = True
+        depNoise = self._insertFunctionStep('pickNoise', prerequisites=self.initDeps)
+        self.newSteps += self.insertExtractPartSteps('NOISE', prerequisites=[depNoise])
+
+      if len(self.readyToExtractMicFns('AND')) >= self.extractingBatch.get() and not self.EXTRACTING['AND']:
+        self.EXTRACTING['AND'] = True
+        depsAnd = self.insertCaculateConsensusSteps('AND', prerequisites=self.initDeps)
+        self.newSteps += self.insertExtractPartSteps('AND', prerequisites=depsAnd)
+
+      trainedParams = self.loadTrainedParams()
+      if self.addTrainingData.get() == self.ADD_DATA_TRAIN_CUST and \
+              self.trainingDataType == self.ADD_DATA_TRAIN_CUSTOM_OPT_COORS:
+        if self.trainTrueSetOfCoords.get() is not None and \
+                len(self.readyToExtractMicFns('ADDITIONAL_COORDS_TRUE')) >= self.extractingBatch.get():
+          self.TO_EXTRACT_MICFNS['ADDITIONAL_COORDS_TRUE'] = self.readyToExtractMicFns('ADDITIONAL_COORDS_TRUE')
+          self.newSteps += self.insertExtractPartSteps('ADDITIONAL_COORDS_TRUE', prerequisites=self.initDeps)
+        if self.trainFalseSetOfCoords.get() is not None and \
+                len(self.readyToExtractMicFns('ADDITIONAL_COORDS_FALSE')) >= self.extractingBatch.get():
+          self.TO_EXTRACT_MICFNS['ADDITIONAL_COORDS_FALSE'] = self.readyToExtractMicFns('ADDITIONAL_COORDS_FALSE')
+          self.newSteps += self.insertExtractPartSteps('ADDITIONAL_COORDS_FALSE', prerequisites=self.initDeps)
+
+      if self.cnnFree():
+        self.TO_TRAIN_MICFNS = self.readyToTrainMicFns()
+        if len(self.TO_TRAIN_MICFNS) >= self.trainingBatch.get():
+          self.TRAINING = True
+          self.curTrainedParams = trainedParams
+          self.depsTrain = [self._insertFunctionStep('trainCNN', self.TO_TRAIN_MICFNS, prerequisites=self.initDeps)]
+          self.depsTrain = [self._insertFunctionStep('endTrainingStep', prerequisites=self.depsTrain)]
+          self.newSteps += self.depsTrain
 
     def trainCNN(self, toTrainMicFns):
         '''Trains the CNN with the particles from the ready to train micrographs'''
@@ -918,7 +898,7 @@ def trainCNN(self, toTrainMicFns):
         for mode in ['AND', 'NOISE']:
           self.joinSetOfParticlesStep(mode, toTrainMicFns, trPass, clean=True)
         #Creatting the training pass directory
-        netDataPath = self._getExtraPath("nnetData{}".format(trPass))
+        netDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass))
         if not os.path.exists(netDataPath):
           makePath(netDataPath)
         nEpochs = self.nEpochs.get()
@@ -1005,9 +985,9 @@ def predictCNN(self):
           self.joinSetOfParticlesStep(mode='OR', micFns=toPredictMicFns, trainingPass=predExten, clean=True)
           mdORPath = self._getExtraPath("particles_OR{}.xmd".format(predExten))
 
-        netDataPath = self._getExtraPath("nnetData{}".format(trPass))
+        netDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass))
         if not os.path.isdir(netDataPath) and self._doContinue():
-            prevRunPath = self.continueRun.get()._getExtraPath("nnetData{}".format(trPass))
+            prevRunPath = self.continueRun.get()._getExtraPath(self.NET_TEMPLATE.format(trPass))
             copyTree(prevRunPath, netDataPath)
         elif self.skipTraining.get() and self._usePretrainedModel():
           self.__retrievePreTrainedModel(netDataPath)
@@ -1032,7 +1012,7 @@ def predictCNN(self):
             posTestDict = None
             negTestDict = None
 
-        outParticlesPath = self._getPath("particles{}.xmd".format(predExten))
+        outParticlesPath = self._getPath(self.PARTICLES_TEMPLATE.format(predExten))
         fnamesPred, weightsPred= self.__dataDict_toStrs(predictDict)
 
         args= " -n %s --mode score -i %s -o %s "%(netDataPath, fnamesPred, outParticlesPath)
@@ -1069,9 +1049,8 @@ def createOutputStep(self, closeStream=False):
     def createFinalOutput(self, closeStream=False):
       predExten = '_partial'
       partSet = self._createSetOfParticles("outputParts_tmp{}".format(predExten))
-      readSetOfParticles(self._getPath("particles{}.xmd".format(predExten)), partSet)
-      inputSampling = self.inputCoordinates[0].get().getMicrographs().getSamplingRate()
-      partSet.setSamplingRate(self._getDownFactor() * inputSampling)
+      readSetOfParticles(self._getPath(self.PARTICLES_TEMPLATE.format(predExten)), partSet)
+      partSet.setSamplingRate(self._getDownFactor() * self.inSamplingRate)
 
       self.outputParticles, self.outputCoordinates = self.getParticlesOutput(partSet), self.getCoordinatesOutput()
 
@@ -1089,17 +1068,15 @@ def createFinalOutput(self, closeStream=False):
           self.outputCoordinates.append(coord)
           self.outputParticles.append(newPart)
 
-      cleanPattern(self._getPath("particles{}.xmd".format(predExten)))
+      cleanPattern(self._getPath(self.PARTICLES_TEMPLATE.format(predExten)))
       cleanPattern(self._getPath("*outputParts_tmp{}.sqlite".format(predExten)))
       writeSetOfParticles(self.outputParticles, self._getPath("particles.xmd"))
       self.updateOutput(closeStream)
 
     def createPreliminarOutput(self, trPass):
       partSet = self._createSetOfParticles("outputParts_tmp{}".format(trPass))
-      readSetOfParticles(self._getPath("particles{}.xmd".format(trPass)), partSet)
-      inputSampling = self.inputCoordinates[0].get().getMicrographs().getSamplingRate()
-      partSet.setSamplingRate(self._getDownFactor() * inputSampling)
-      boxSize = self._getBoxSize()
+      readSetOfParticles(self._getPath(self.PARTICLES_TEMPLATE.format(trPass)), partSet)
+      partSet.setSamplingRate(self._getDownFactor() * self.inSamplingRate)
 
       self.preliminarOutputParticles = self.getPreParticlesOutput(partSet)
       self.preliminarOutputCoordinates = self.getPreCoordinatesOutput()
@@ -1117,9 +1094,9 @@ def createPreliminarOutput(self, trPass):
           self.preliminarOutputCoordinates.append(coord)
           self.preliminarOutputParticles.append(part)
 
-      cleanPattern(self._getPath("particles{}.xmd".format(trPass)))
+      cleanPattern(self._getPath(self.PARTICLES_TEMPLATE.format(trPass)))
       cleanPattern(self._getPath("*outputParts_tmp{}.sqlite".format(trPass)))
-      writeSetOfParticles(self.preliminarOutputParticles, self._getPath("particles{}.xmd".format(trPass)))
+      writeSetOfParticles(self.preliminarOutputParticles, self._getPath(self.PARTICLES_TEMPLATE.format(trPass)))
       self.updatePreOutput(closeStream=True)
 
     def getPreCoordinatesOutput(self):
@@ -1129,8 +1106,10 @@ def getPreCoordinatesOutput(self):
       self.preliminarOutputCoordinates.setBoxSize(self._getBoxSize())
       self.preliminarOutputCoordinates.setStreamState(SetOfParticles.STREAM_OPEN)
       self._defineOutputs(preliminarOutputCoordinates=self.preliminarOutputCoordinates)
+      self.waitFreeInputCoords()
       for inSetOfCoords in self.inputCoordinates:
         self._defineSourceRelation(inSetOfCoords.get(), self.preliminarOutputCoordinates)
+      self.USING_INPUT_COORDS = False
       return self.preliminarOutputCoordinates
 
     def getCoordinatesOutput(self):
@@ -1141,8 +1120,10 @@ def getCoordinatesOutput(self):
         self.outputCoordinates.setBoxSize(self._getBoxSize())
         self.outputCoordinates.setStreamState(SetOfParticles.STREAM_OPEN)
         self._defineOutputs(outputCoordinates=self.outputCoordinates)
+        self.waitFreeInputCoords()
         for inSetOfCoords in self.inputCoordinates:
           self._defineSourceRelation(inSetOfCoords.get(), self.outputCoordinates)
+        self.USING_INPUT_COORDS = False
       else:
         # Micrographs of the set removed because there might be new ones in streaming
         self.outputCoordinates.setMicrographs(self.coordinatesDict['OR'].getMicrographs())
@@ -1243,26 +1224,29 @@ def getMicrographFnsWithCoordinates(self, shared=True):
       '''Return a list with the filenames of those microgrpahs which already have coordinates associated in the input
       sets. If shared, it must be in all the sets, if not shared, at least in one'''
       sharedMics = self.getAllCoordsInputMicrographs(shared)
+      self.waitFreeInputCoords()
       micPaths = []
-      for mic in sharedMics:
-        coordsInMic = []
+      for micFn in sharedMics:
+        coordsInMic, mic = [], sharedMics[micFn]
         for coordSet in self.inputCoordinates:
           for coord in coordSet.get().iterCoordinates(mic):
-            coordsInMic.append(True)
+            coordsInMic.append(coord)
             break
 
         if len(coordsInMic) == len(self.inputCoordinates):
           micPaths.append(mic.getFileName())
+      self.USING_INPUT_COORDS = False
 
       micFns = self.prunePaths(micPaths)
       return micFns
 
     def getAllCoordsInputMicrographs(self, shared=False):
-      '''Returns a set with the input micrographs present associated with all the input coordinates sets.
+      '''Returns a dic {micFn: mic} with the input micrographs present associated with all the input coordinates sets.
       If shared, the list contains only those micrographs present in all input coordinates sets, else the list contains
-      all microgrpah present in any set (Intersection vs Union)'''
-      micFns = set([])
-      micDict = {}
+      all microgrpah present in any set (Intersection vs Union)
+      Do not create a set, because of concurrency in the database'''
+      self.waitFreeInputCoords()
+      micDict, micFns = {}, set([])
       for inputCoord in self.inputCoordinates:
         newMics = inputCoord.get().getMicrographs()
         newMicFns = []
@@ -1275,17 +1259,17 @@ def getAllCoordsInputMicrographs(self, shared=False):
           micFns = micFns | set(newMicFns)
         else:
           micFns = micFns & set(newMicFns)
-
-      mics = self._createSetOfMicrographs()
-      mics.copyInfo(self.inputCoordinates[0].get().getMicrographs())
+      self.USING_INPUT_COORDS = False
+      sharedMicDict = {}
       for micFn in micFns:
-        mics.append(micDict[micFn])
+        sharedMicDict[micFn] = micDict[micFn]
 
-      return mics
+      return sharedMicDict
 
     def _getInputMicrographs(self):
       '''Return a list with the micrographs corresponding the input coordinates'''
       if not hasattr(self, "inputMicrographs") or not self.inputMicrographs:
+        self.waitFreeInputCoords()
         if len(self.inputCoordinates) == 0:
           print("WARNING. PROVIDE MICROGRAPHS FIRST")
         else:
@@ -1293,12 +1277,15 @@ def _getInputMicrographs(self):
           if inputMicrographs is None:
             raise ValueError("there are problems with your coordiantes, they do not have associated micrographs ")
           self.inputMicrographs = inputMicrographs
+        self.USING_INPUT_COORDS = False
       return self.inputMicrographs
 
     def _getBoxSize(self):
       '''Returns the box size of the input coordinates'''
       if not hasattr(self, "boxSize") or not self.boxSize:
+        self.waitFreeInputCoords()
         firstCoords = self.inputCoordinates[0].get()
+        self.USING_INPUT_COORDS = False
         self.boxSize = firstCoords.getBoxSize()
         self.downFactor = self.boxSize / float(DEEP_PARTICLE_SIZE)
       return self.boxSize
@@ -1340,55 +1327,16 @@ def getPreProcParamsFromForm(self):
                     "mics_skipInvert": self.skipInvert.get(),
                     "mics_pathToMics": pathToMics,
                     "mics_pathToCtfs": pathToCtfs}
+        self.waitFreeInputCoords()
         coordsNames=[]
         for inputCoords in self.inputCoordinates:
           coordsNames.append( inputCoords.get().getFileName() )
+        self.USING_INPUT_COORDS = False
         coordsNames= tuple(sorted(coordsNames))
         paramsInfo["coords_pathCoords"]= coordsNames
 
         return paramsInfo
 
-    def checkIfPrevRunIsCompatible(self, inputType=""):
-        '''
-        inputType can be mics_ or coords_ or ""
-        '''
-        def _makeTupleIfList(candidateList):
-          if isinstance(candidateList, list):
-            return tuple(candidateList)
-          elif isinstance(candidateList, six.string_types):
-            return tuple([candidateList])
-          else:
-            return candidateList
-        try:
-          preprocParams= self.getPreProcParamsFromForm()
-        except Exception:
-          return False
-        preprocParams= { k:_makeTupleIfList(preprocParams[k]) for k in preprocParams if k.startswith(inputType) }
-        if self._doContinue():
-            preprocessParamsFname = self.continueRun.get()._getExtraPath("preprocess_params.json")
-            with open(preprocessParamsFname) as f:
-              preprocParams_loaded = json.load(f)
-
-            preprocParams_loaded= {k: _makeTupleIfList(preprocParams_loaded[k])
-                                   for k in preprocParams_loaded
-                                   if k.startswith(inputType)}
-
-            for key in preprocParams_loaded:
-              if "path" in key:
-                for val in preprocParams_loaded[key]:
-                  if not val=="None" and not os.path.exists(val):
-                    return False
-
-            shared_items = {k: preprocParams_loaded[k]
-                            for k in preprocParams_loaded
-                            if k in preprocParams and preprocParams_loaded[k] == preprocParams[k]}
-
-            return (len(shared_items) == len(preprocParams) and
-                    len(shared_items) == len(preprocParams_loaded) and
-                    len(shared_items) > 0)
-
-        return False
-
     def __dataDict_toStrs(self, dataDict):
         fnamesStr=[]
         weightsStr=[]
@@ -1415,7 +1363,7 @@ def predictionsOn(self):
     def readyPreliminarPrediction(self):
       '''Return a boolean for whether to perform a preliminar predition. True if the user set it and the current
       trained network has not been used yet'''
-      if self.networkReadyToPredict:
+      if self.networkReadyToPredict():
         trainedParams = self.loadTrainedParams()
         if not trainedParams['trainingPass'] in trainedParams[
           'predictionPasses'] and self.doPreliminarPredictions.get():
@@ -1425,7 +1373,8 @@ def readyPreliminarPrediction(self):
     def networkReadyToPredict(self):
       '''Returns true if the CNN is trained or the user specified it does not need to be trained'''
       trainedParams = self.loadTrainedParams()
-      return self.skipTraining.get() or len(trainedParams['trainedMicFns']) > 0
+      return (self.skipTraining.get() and len(self.readyToPredictMicFns()) > 0) \
+             or len(trainedParams['trainedMicFns']) > 0
 
     def cnnFree(self):
       return not self.PREDICTING and not self.TRAINING
@@ -1441,12 +1390,16 @@ def allFree(self):
 
     def checkIfParentsFinished(self):
       '''Check the streamState of the coordinates input to check if the parent protocols are finsihed'''
+      self.waitFreeInputCoords()
+      finished=True
       for coords in self.inputCoordinates:
         coords = coords.get()
         coords.loadAllProperties()
         if coords.isStreamOpen():
-          return False
-      return True
+          finished = False
+          break
+      self.USING_INPUT_COORDS = False
+      return finished
 
     def checkIfNewMics(self, mode=''):
       '''Check if the are new micrographs ready for extracting particles'''
@@ -1464,21 +1417,18 @@ def getMicsIds(self, filterOutNoCoords=False):
         '''Returns the input micrographs Ids'''
         if not filterOutNoCoords:
           return self._getInputMicrographs().getIdSet()
-        micIds= set([])
-        micFnames= set([])
+        self.waitFreeInputCoords()
+        micFnames, micIds = set([]), set([])
         for coordinatesP in self.inputCoordinates:
             for coord in coordinatesP.get():
               micIds.add( coord.getMicId())
               micFnames.add( coord.getMicName() )
+        self.USING_INPUT_COORDS = False
         return sorted( micIds )
 
     def getInputMicsFns(self, shared):
       '''Returns the input micrographs filenames'''
-      micFnames = []
-      mics = self.getAllCoordsInputMicrographs(shared)
-      micPaths = mics.getFiles()
-      micFnames = self.prunePaths(micPaths)
-      return sorted(micFnames)
+      return sorted(self.getAllCoordsInputMicrographs(shared).keys())
 
     def prunePaths(self, paths):
       fns = []
@@ -1541,6 +1491,27 @@ def readyToPredictMicFns(self):
       readyToPredict = list(extractedMicFns - set(self.getPredictedMicFns()))
       return readyToPredict
 
+    def getInpCoordsFns(self, mode, extractedSetOfCoordsFns):
+      Tm = []
+      for coordinatesP in self.inputCoordinates:
+        mics = coordinatesP.get().getMicrographs()
+        Tm.append(mics.getSamplingRate())
+      nCoordsSets = len(Tm)
+
+      inputCoordsFnames = {}
+      for coord_num, coordinatesP in enumerate(self.inputCoordinates):
+        tmpPosDir = self._getTmpPath("input_coords_%d_%s" % (coord_num, mode))
+        if not os.path.exists(tmpPosDir):
+          makePath(tmpPosDir)
+        writeSetOfCoordinates(tmpPosDir, coordinatesP.get(), scale=float(Tm[coord_num]) / float(Tm[0]))
+        for posFname in os.listdir(tmpPosDir):
+          baseName, extension = os.path.splitext(os.path.basename(posFname))
+          if extension == ".pos" and not posFname in extractedSetOfCoordsFns:
+            if baseName not in inputCoordsFnames:
+              inputCoordsFnames[baseName] = ["None"] * nCoordsSets
+            inputCoordsFnames[baseName][coord_num] = os.path.join(tmpPosDir, posFname)
+      return inputCoordsFnames
+
     #Training params utils
     def loadTrainedParams(self):
       '''Load the dictionary stored in pickle format which stores the trained parameters.
@@ -1575,17 +1546,17 @@ def saveTrainedParams(self, params):
     #CNN models utils
     def retrievePreviousPassModel(self, trPass, lastTrPass=''):
       '''Retrieves a previous CNN model and copies its folders to used the network in a new location'''
-      curNetDataPath = self._getExtraPath("nnetData{}".format(trPass))
+      curNetDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass))
       if trPass == '':
-        prevNetDataPath = self._getExtraPath("nnetData{}".format(lastTrPass))
+        prevNetDataPath = self._getExtraPath(self.NET_TEMPLATE.format(lastTrPass))
       else:
-        prevNetDataPath = self._getExtraPath("nnetData{}".format(trPass - 1))
+        prevNetDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass - 1))
       if prevNetDataPath != curNetDataPath:
         copyTree(prevNetDataPath, curNetDataPath)
 
     def retrievePreviousRunModel(self, prevProt, trPass=''):
       '''Retrieves a CNN model from other protocol and copies its folders to used the network in a new location'''
-      curNetDataPath = self._getExtraPath("nnetData{}".format(trPass))
+      curNetDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass))
       prevNetDataPath = prevProt._getExtraPath("nnetData")
       if prevNetDataPath != curNetDataPath:
         copyTree(prevNetDataPath, curNetDataPath)
diff --git a/xmipp3/protocols/protocol_simulate_ctf.py b/xmipp3/protocols/protocol_simulate_ctf.py
index 3c0f1adf9..e5b9244fa 100644
--- a/xmipp3/protocols/protocol_simulate_ctf.py
+++ b/xmipp3/protocols/protocol_simulate_ctf.py
@@ -1,6 +1,7 @@
 # **************************************************************************
 # *
 # * Authors:     C.O.S. Sorzano (coss@cnb.csic.es)
+# *              Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -32,32 +33,27 @@
 from pwem.objects import CTFModel
 
 from xmippLib import createEmptyFile
-from xmipp3.convert import locationToXmipp
 
-from pwem import emlib
-from xmipp3.convert import (setXmippAttributes, writeSetOfParticles)
 
 class XmippProtSimulateCTF(Prot2D):
-    """    
+    """
     Simulate the effect of the CTF (no amplitude decay).
     A random defocus is chosen between the lower and upper defocus for each projection.
     """
 
     _label = 'simulate ctf'
     _lastUpdateVersion = VERSION_3_0
-    
+
     def __init__(self, *args, **kwargs):
         Prot2D.__init__(self, *args, **kwargs)
-        
-    #--------------------------- DEFINE param functions ------------------------
+
+    # --------------------------- DEFINE param functions ------------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
-        
         form.addParam('inputParticles', params.PointerParam,
                       allowsNull=False,
                       pointerClass='SetOfParticles',
                       label="Input particles")
-
         form.addParam('voltage', params.FloatParam, default=300,
                       label="Voltage (kV)")
         form.addParam('cs', params.FloatParam, default=2.7,
@@ -70,28 +66,20 @@ def _defineParams(self, form):
         form.addParam('DefocusF', params.FloatParam, default=25000,
                       label="Upper defocus (A)")
 
-    
-    #--------------------------- INSERT steps functions ------------------------
-
-    def _insertAllSteps(self):        
-        self._insertFunctionStep('convertInputStep',
-                                 self.inputParticles.get().getObjId())
-
-        self._insertFunctionStep('simulateStep',
-                                 self.inputParticles.get().getObjId())
-
-        self._insertFunctionStep('createOutputStep')
+    # --------------------------- INSERT steps functions ------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('convertInputStep')
+        self._insertFunctionStep('simulateStep')
 
-    #--------------------------- STEPS functions -------------------------------
-
-    def convertInputStep(self, particlesId):
+    # --------------------------- STEPS functions -------------------------------
+    def convertInputStep(self):
         x, y, _ = self.inputParticles.get().getDimensions()
         n = self.inputParticles.get().getSize()
-        createEmptyFile(self._getPath("images.stk"),x,y,1,n)
+        createEmptyFile(self._getPath("images.mrc"), x, y, 1, n)
 
-    def simulateStep(self, particlesId):
-        n=1
-        fnStk = self._getPath("images.stk")
+    def simulateStep(self):
+        n = 1
+        fnStk = self._getPath("images.mrc")
         Ts = self.inputParticles.get().getSamplingRate()
 
         imgSetOut = self._createSetOfParticles()
@@ -104,56 +92,38 @@ def simulateStep(self, particlesId):
         acquisition.setMagnification(1)
 
         for particle in self.inputParticles.get():
-
             location = particle.getLocation()
-            fnIn = str(location[0])+"@"+location[1]
-            fnOut = str(n)+"@"+fnStk
-            defocus = random.uniform(self.Defocus0.get(),self.DefocusF.get())
-            args = "-i %s -o %s"%(fnIn, fnOut)
-            args += " --fourier ctfdef %f %f %f %f --sampling %f -v 0"%(self.voltage, self.cs, self.Q0, defocus, Ts)
+            fnIn = str(location[0]) + "@" + location[1]
+            fnOut = str(n) + "@" + fnStk
+            defocus = random.uniform(self.Defocus0.get(), self.DefocusF.get())
+            args = "-i %s -o %s" % (fnIn, fnOut)
+            args += " --fourier ctfdef %f %f %f %f --sampling %f -v 0" % (self.voltage, self.cs, self.Q0, defocus, Ts)
             self.runJob("xmipp_transform_filter", args)
 
             newCTF = CTFModel()
             newCTF.setDefocusU(defocus)
             newCTF.setDefocusV(defocus)
             newCTF.setDefocusAngle(0.0)
-
             newParticle = particle.clone()
+            newParticle.setLocation((n, fnStk))
             acquisition = newParticle.getAcquisition()
             acquisition.setVoltage(self.voltage.get())
             acquisition.setAmplitudeContrast(self.Q0.get())
             acquisition.setSphericalAberration(self.cs.get())
             acquisition.setMagnification(1)
-
             newParticle.setCTF(newCTF)
-
             imgSetOut.append(newParticle)
-
-            n+=1
+            n += 1
 
         self._defineOutputs(outputParticles=imgSetOut)
         self._defineSourceRelation(self.inputParticles.get(), imgSetOut)
 
-    def createOutputStep(self):
-        # imgSet1 = self.inputParticles1.get()
-        # imgSetOut = self._createSetOfParticles()
-        # imgSetOut.copyInfo(imgSet1)
-        # imgSetOut.setAlignmentProj()
-        # self.iterMd = md.iterRows(self._getExtraPath("angles1_common.xmd"), md.MDL_ITEM_ID)
-        # self.lastRow = next(self.iterMd)
-        # imgSetOut.copyItems(imgSet1,
-        #                     updateItemCallback=self._updateItem)
-        # self._defineOutputs(outputParticles=imgSetOut)
-        # self._defineSourceRelation(self.inputParticles1, imgSetOut)
-        # self._defineSourceRelation(self.inputParticles2, imgSetOut)
-        pass
-
     # --------------------------- INFO functions -------------------------------
     def _summary(self):
         summary = []
-        summary.append("Voltage=%f kV"%self.voltage)
-        summary.append("Cs=%f mm"%self.cs)
-        summary.append("Q0=%f"%self.Q0)
-        summary.append("Defocus range=[%f,%f] A"%(self.Defocus0,self.DefocusF))
+        summary.append("Voltage=%f kV" % self.voltage)
+        summary.append("Cs=%f mm" % self.cs)
+        summary.append("Q0=%f" % self.Q0)
+        summary.append("Defocus range=[%f,%f] A" % (self.Defocus0, self.DefocusF))
         return summary
 
diff --git a/xmipp3/protocols/protocol_structure_map.py b/xmipp3/protocols/protocol_structure_map.py
new file mode 100644
index 000000000..72416717e
--- /dev/null
+++ b/xmipp3/protocols/protocol_structure_map.py
@@ -0,0 +1,201 @@
+
+# **************************************************************************
+# *
+# * Authors:     David Herreros Calero (dherreros@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
+import os
+
+import xmippLib
+
+from pwem.protocols import ProtAnalysis3D
+import pyworkflow.protocol.params as params
+from pwem.emlib.image import ImageHandler
+from pwem.objects import SetOfVolumes
+from pyworkflow import VERSION_2_0
+from pyworkflow.utils import cleanPath
+
+from ..protocols.protocol_structure_map_zernike3d import mds
+
+
+class XmippProtStructureMap(ProtAnalysis3D):
+    """ Protocol for structure mapping based on correlation distance. """
+    _label = 'struct map'
+    _lastUpdateVersion = VERSION_2_0
+    OUTPUT_SUFFIX = '_%d_crop.vol'
+
+    def __init__(self, **args):
+        ProtAnalysis3D.__init__(self, **args)
+        self.stepsExecutionMode = params.STEPS_PARALLEL
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('inputVolumes', params.PointerParam,
+                      pointerClass='SetOfClasses3D, SetOfVolumes',
+                      label="Input volume(s)", important=True,
+                      help='Select one or more volumes (SetOfClasses3D)\n'
+                           'for structure mapping.')
+        form.addParam('targetResolution', params.FloatParam, label="Target resolution",
+                      default=8.0,
+                      help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
+                           "2/3 of the Fourier spectrum.")
+        form.addParallelSection(threads=1, mpi=1)
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        volList = []
+        dimList = []
+        srList = []
+        volList, dimList, srList, _ = self._iterInputVolumes(volList, dimList, srList, [])
+        nVoli = 1
+        depsConvert = []
+        for _ in volList:
+            convert = self._insertFunctionStep('convertStep', volList[nVoli - 1],
+                                     dimList[nVoli - 1], srList[nVoli - 1],
+                                     min(dimList), max(srList), nVoli, prerequisites=[])
+            depsConvert.append(convert)
+            nVoli += 1
+
+        nVoli = 1
+        deps = []
+        for _ in volList:
+            nVolj = 1
+            for _ in volList:
+                if nVolj != nVoli:
+                    stepID = self._insertFunctionStep('alignStep', volList[nVoli - 1],
+                                            volList[nVolj - 1], nVoli - 1,
+                                            nVolj - 1, prerequisites=depsConvert)
+                    deps.append(stepID)
+                nVolj += 1
+            nVoli += 1
+
+        self._insertFunctionStep('correlationMatrix', volList, prerequisites=deps)
+        self._insertFunctionStep('gatherResultsStepCorr')
+
+    # --------------------------- STEPS functions ---------------------------------------------------
+    def convertStep(self, volFn, volDim, volSr, minDim, maxSr, nVoli):
+        Xdim = volDim
+        Ts = volSr
+        newTs = self.targetResolution.get() * 1.0 / 3.0
+        newTs = max(maxSr, newTs)
+        newXdim = Xdim * Ts / newTs
+        fnOut = os.path.splitext(volFn)[0]
+        fnOut = self._getExtraPath(os.path.basename(fnOut + self.OUTPUT_SUFFIX % nVoli))
+
+        ih = ImageHandler()
+        if Xdim != newXdim:
+            self.runJob("xmipp_image_resize",
+                        "-i %s -o %s --dim %d" % (volFn, fnOut, newXdim))
+
+        else:
+            ih.convert(volFn, fnOut)
+
+        if newXdim>minDim:
+            self.runJob("xmipp_transform_window", " -i %s -o %s --crop %d" %
+                        (fnOut, fnOut, (newXdim - minDim)))
+
+    def alignStep(self, inputVolFn, refVolFn, i, j):
+        inputVolFn = self._getExtraPath(os.path.basename(os.path.splitext(inputVolFn)[0]
+                                                         + self.OUTPUT_SUFFIX % (i + 1)))
+        refVolFn = self._getExtraPath(os.path.basename(os.path.splitext(refVolFn)[0]
+                                                       + self.OUTPUT_SUFFIX % (j + 1)))
+        fnOut = self._getTmpPath('vol%dAlignedTo%d.vol' % (i, j))
+
+        params = ' --i1 %s --i2 %s --apply %s --local --dontScale' % \
+                 (refVolFn, inputVolFn, fnOut)
+
+        self.runJob("xmipp_volume_align", params)
+
+        self.computeCorr(fnOut, refVolFn, i, j)
+
+        cleanPath(fnOut)
+
+    def computeCorr(self, vol1, vol2, i, j):
+        vol = xmippLib.Image(vol1)
+        defVol = xmippLib.Image(vol2)
+        corr = vol.correlation(defVol)
+        corr = 1 - corr
+        outFile = self._getExtraPath('Pair_%d_%d_correlation.txt' % (i,j))
+        with open(outFile, 'w') as f:
+            f.write('%f' % corr)
+
+    def correlationMatrix(self, volList):
+        numVol = len(volList)
+        self.corrMatrix = np.zeros((numVol, numVol))
+        for i in range(numVol):
+            for j in range(numVol):
+                if i != j:
+                    path = self._getExtraPath('Pair_%d_%d_correlation.txt' % (i,j))
+                    self.corrMatrix[i,j] = np.loadtxt(path)
+
+    def gatherResultsStepCorr(self):
+        fnRoot = self._getExtraPath("CorrMatrix.txt")
+        self.saveCorrelation(self.corrMatrix, fnRoot)
+
+    def saveCorrelation(self, matrix, fnRoot, label=''):
+        np.savetxt(fnRoot, matrix, "%f")
+        corr = np.asarray(matrix)
+        for i in range(1, 4):
+            embed = mds(corr, i)
+            embedExtended = np.pad(embed, ((0, 0), (0, i - embed.shape[1])),
+                                   "constant", constant_values=0)
+            np.savetxt(self._defineResultsName(i, label), embedExtended)
+
+    # --------------------------- UTILS functions --------------------------------------------
+    def _iterInputVolumes(self, volList, dimList, srList, idList):
+        """ Iterate over all the input volumes. """
+        count = 1
+        weight = []
+        if not isinstance(self.inputVolumes.get(), SetOfVolumes):
+            for cls in self.inputVolumes.get().iterItems():
+                vol = cls.getRepresentative()
+                volList.append(vol.getFileName())
+                dimList.append(vol.getDim()[0])
+                srList.append(vol.getSamplingRate())
+                idList.append(count)
+                weight.append(cls.getSize())
+                count += 1
+        else:
+            for vol in self.inputVolumes.get().iterItems():
+                volList.append(vol.getFileName())
+                dimList.append(vol.getDim()[0])
+                srList.append(vol.getSamplingRate())
+                idList.append(count)
+                weight.append(1)
+                count += 1
+        np.savetxt(self._getExtraPath('weights.txt'), np.asarray(weight))
+        return volList, dimList, srList, idList
+
+    def split(self, array, nrows, ncols):
+        """Split a matrix into sub-matrices."""
+        r, h = array.shape
+        return (array.reshape(h // nrows, nrows, -1, ncols)
+                .swapaxes(1, 2)
+                .reshape(-1, nrows, ncols))
+
+    def _defineResultsName(self, i, label=''):
+        return self._getExtraPath('Coordinate%sMatrixCorr%d.txt' % (label,i))
diff --git a/xmipp3/protocols/protocol_structure_map_sph.py b/xmipp3/protocols/protocol_structure_map_sph.py
deleted file mode 100644
index c58d9de2c..000000000
--- a/xmipp3/protocols/protocol_structure_map_sph.py
+++ /dev/null
@@ -1,198 +0,0 @@
-########
-
-# This protocol has been removed from Scipion as there is a newer version being developed that
-# will solve some of the problems of this current implementation. The code will be available
-# again in a newer release.
-
-########
-
-
-
-# **************************************************************************
-# *
-# * Authors:     Amaya Jimenez Moreno (ajimenez@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 pwem.protocols import ProtAnalysis3D
-# import pyworkflow.protocol.params as params
-# from pwem.emlib.image import ImageHandler
-# from pwem.objects import SetOfVolumes, Volume
-# from pyworkflow import VERSION_2_0
-# import numpy as np
-#
-#
-# def mds(d, dimensions=2):
-#     """
-#     Multidimensional Scaling - Given a matrix of interpoint distances,
-#     find a set of low dimensional points that have similar interpoint
-#     distances.
-#     """
-#
-#     (n, n) = d.shape
-#     E = (-0.5 * d ** 2)
-#
-#     # Use mat to get column and row means to act as column and row means.
-#     Er = np.mat(np.mean(E, 1))
-#     Es = np.mat(np.mean(E, 0))
-#
-#     # From Principles of Multivariate Analysis: A User's Perspective (page 107).
-#     F = np.array(E - np.transpose(Er) - Es + np.mean(E))
-#
-#     [U, S, V] = np.linalg.svd(F)
-#
-#     Y = U * np.sqrt(S)
-#
-#     return (Y[:, 0:dimensions], S)
-#
-#
-# class XmippProtStructureMapSPH(ProtAnalysis3D):
-#     """ Protocol for structure mapping based on spherical harmonics. """
-#     _label = 'sph struct map'
-#     _lastUpdateVersion = VERSION_2_0
-#
-#     # --------------------------- DEFINE param functions --------------------------------------------
-#     def _defineParams(self, form):
-#         form.addSection(label='Input')
-#         form.addParam('inputVolumes', params.MultiPointerParam,
-#                       pointerClass='SetOfVolumes,Volume',
-#                       label="Input volume(s)", important=True,
-#                       help='Select one or more volumes (Volume or SetOfVolumes)\n'
-#                            'for structure mapping.')
-#         form.addParam('targetResolution', params.FloatParam, label="Target resolution",
-#                       default=8.0,
-#                       help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
-#                            "2/3 of the Fourier spectrum.")
-#         form.addParam('depth', params.IntParam, default=3,
-#                       label='Harmonical depth',
-#                       expertLevel=params.LEVEL_ADVANCED,
-#                       help='Harmonical depth of the deformation=1,2,3,...')
-#
-#     # --------------------------- INSERT steps functions --------------------------------------------
-#     def _insertAllSteps(self):
-#
-#         volList, dimList, srList = self._iterInputVolumes()
-#         self.distanceMatrix = np.zeros((len(volList), len(volList)))
-#
-#         nVoli = 1
-#         for voli in volList:
-#             self._insertFunctionStep('convertStep', volList[nVoli - 1],
-#                                      dimList[nVoli - 1], srList[nVoli - 1],
-#                                      min(dimList), max(srList))
-#             nVolj = 1
-#             for volj in volList:
-#                 if nVolj != nVoli:
-#                     self._insertFunctionStep('deformStep', volList[nVoli - 1],
-#                                              volList[nVolj - 1], nVoli - 1,
-#                                              nVolj - 1)
-#                 else:
-#                     self.distanceMatrix[nVoli - 1][nVolj - 1] = 0.0
-#                 nVolj += 1
-#             nVoli += 1
-#
-#         self._insertFunctionStep('gatherResultsStep', volList)
-#
-#     # --------------------------- STEPS functions ---------------------------------------------------
-#     def convertStep(self, volFn, volDim, volSr, minDim, maxSr):
-#         Xdim = volDim
-#         Ts = volSr
-#         newTs = self.targetResolution.get() * 1.0 / 3.0
-#         newTs = max(maxSr, newTs)
-#         newXdim = int(Xdim * Ts / newTs)
-#
-#         ih = ImageHandler()
-#         if Xdim != newXdim:
-#             self.runJob("xmipp_image_resize",
-#                         "-i %s --dim %d" % (volFn, newXdim))
-#
-#         if newXdim>minDim:
-#             self.runJob("xmipp_transform_window", " -i %s --crop %d" %
-#                         (volFn, (newXdim - minDim)))
-#
-#     def deformStep(self, inputVolFn, refVolFn, i, j):
-#         fnOut = self._getExtraPath('vol%dDeformedTo%d.vol' % (i, j))
-#
-#         params = ' --i1 %s --i2 %s --apply %s --least_squares --local ' % \
-#                  (refVolFn, inputVolFn, fnOut)
-#         self.runJob("xmipp_volume_align", params)
-#
-#         params = ' -i %s -r %s -o %s --depth %d ' %\
-#                  (fnOut, refVolFn, fnOut, self.depth.get())
-#
-#         self.runJob("xmipp_volume_deform_sph", params)
-#         distanceValue = np.loadtxt('./deformation.txt')
-#         self.distanceMatrix[i][j] = distanceValue
-#
-#
-#     def gatherResultsStep(self, volList):
-#         fnRoot = self._getExtraPath("DistanceMatrix.txt")
-#         nVoli = 1
-#         for i in volList:
-#             nVolj = 1
-#             for j in volList:
-#                 fh = open(fnRoot, "a")
-#                 fh.write("%f\t" % self.distanceMatrix[(nVoli - 1)][(nVolj - 1)])
-#                 fh.close()
-#                 nVolj += 1
-#             fh = open(fnRoot, "a")
-#             fh.write("\n")
-#             fh.close()
-#             nVoli += 1
-#
-#         distance = np.asarray(self.distanceMatrix)
-#         for i in range(1, 4):
-#             embed, _ = mds(distance, i)
-#             embedExtended = np.pad(embed, ((0, 0), (0, i - embed.shape[1])),
-#                                    "constant", constant_values=0)
-#             np.savetxt(self._defineResultsName(i), embedExtended)
-#
-#     # --------------------------- UTILS functions --------------------------------------------
-#
-#
-#     def _iterInputVolumes(self):
-#         """ Iterate over all the input volumes. """
-#         volList = []
-#         dimList = []
-#         srList = []
-#         for pointer in self.inputVolumes:
-#             item = pointer.get()
-#             if item is None:
-#                 break
-#             itemId = item.getObjId()
-#             if isinstance(item, Volume):
-#                 volList.append(item.getFileName())
-#                 dimList.append(item.getDim()[0])
-#                 srList.append(item.getSamplingRate())
-#             elif isinstance(item, SetOfVolumes):
-#                 for vol in item:
-#                     volList.append(vol.getFileName())
-#                     dimList.append(vol.getDim()[0])
-#                     srList.append(vol.getSamplingRate())
-#         return volList, dimList, srList
-#
-#
-#     def _defineResultsName(self, i):
-#         return self._getExtraPath('CoordinateMatrix%d.txt' % i)
-
-
-
-
diff --git a/xmipp3/protocols/protocol_structure_map_zernike3d.py b/xmipp3/protocols/protocol_structure_map_zernike3d.py
new file mode 100644
index 000000000..3aa77ffa0
--- /dev/null
+++ b/xmipp3/protocols/protocol_structure_map_zernike3d.py
@@ -0,0 +1,523 @@
+
+# **************************************************************************
+# *
+# * Authors:     Amaya Jimenez Moreno (ajimenez@cnb.csic.es)
+# *              David Herreros Calero (dherreros@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 scipy.spatial.distance import pdist
+from scipy.signal import find_peaks
+from scipy.stats import entropy
+from scipy.ndimage.filters import gaussian_filter
+import os
+
+import xmippLib
+
+from pwem.protocols import ProtAnalysis3D
+import pyworkflow.protocol.params as params
+from pwem.emlib.image import ImageHandler
+from pwem.objects import SetOfVolumes, Volume
+from pyworkflow import VERSION_2_0
+from pyworkflow.utils import cleanPath, copyFile, getExt
+
+
+def mds(d, dimensions=2):
+    """
+    Multidimensional Scaling - Given a matrix of interpoint distances,
+    find a set of low dimensional points that have similar interpoint
+    distances.
+    """
+
+    # Distance matrix size
+    (n, n) = d.shape
+
+    # Centering Matrix
+    J = np.identity(n, float) - (1 / n) * np.ones([n, n], float)
+
+    # Center distance matrix
+    B = -0.5 * J @ d @ J
+
+    # Singular value decomposition
+    [U, S, V] = np.linalg.svd(B)
+    S = np.diag(S)
+
+    # Coordinates matrix from MDS
+    Y = U[:, 0:dimensions] @ np.power(S[0:dimensions, 0:dimensions], 0.5)
+
+    return Y
+
+def rigidRegistration(X, Y):
+    Xcm = np.sum(X, axis=0) / X.shape[0]
+    Ycm = np.sum(Y, axis=0) / Y.shape[0]
+    Xc = np.transpose(X - Xcm)
+    Yc = np.transpose(Y - Ycm)
+    [U, S, V] = np.linalg.svd(Xc @ Yc.T)
+    R = V @ U.T
+    t = Ycm - R @ Xcm
+    return R, t
+
+
+class XmippProtStructureMapZernike3D(ProtAnalysis3D):
+    """ Protocol for structure mapping based on Zernike3D. """
+    _label = 'struct map - Zernike3D'
+    _lastUpdateVersion = VERSION_2_0
+    OUTPUT_SUFFIX = '_%d_crop.vol'
+    ALIGNED_VOL = 'vol%dAligned.vol'
+
+    def __init__(self, **args):
+        ProtAnalysis3D.__init__(self, **args)
+        self.stepsExecutionMode = params.STEPS_PARALLEL
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addHidden(params.USE_GPU, params.BooleanParam, default=True,
+                       label="Use GPU for execution",
+                       help="This protocol has both CPU and GPU implementation.\
+                           Select the one you want to use.")
+        form.addHidden(params.GPU_LIST, params.StringParam, default='0',
+                       expertLevel=params.LEVEL_ADVANCED,
+                       label="Choose GPU IDs",
+                       help="Add a list of GPU devices that can be used")
+        form.addParam('inputVolumes', params.MultiPointerParam,
+                      pointerClass='SetOfVolumes,Volume',
+                      label="Input volume(s)", important=True,
+                      help='Select one or more volumes (Volume or SetOfVolumes)\n'
+                           'for structure mapping.')
+        form.addParam('twoSets', params.BooleanParam, label='Compare two sets?', default=False,
+                      help='Useful when two Sets are intended to be compared independently (e.g. '
+                           'comparing EMDBS and Maps coming from PDBs).')
+        form.addParam('secondSet', params.MultiPointerParam, pointerClass='SetOfVolumes,Volume',
+                      condition='twoSets==True', label='Second set of volumes', allowsNull=True,
+                      help='Select one or more volumes (Volume or SetOfVolumes)\n'
+                           'to compare to the first set.')
+        form.addParam('targetResolution', params.FloatParam, label="Target resolution",
+                      default=8.0,
+                      help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
+                           "2/3 of the Fourier spectrum.")
+        form.addParam('sigma', params.NumericListParam, label="Multiresolution", default="1 2",
+                      help="Perform the analysis comparing different filtered versions of the volumes. The values "
+                           "specified here will determine the cutoff frequency of the filter in normalized units "
+                           "(normalized to 1/2).")
+        form.addParam('Rmax', params.IntParam, default=0,
+                      label='Sphere radius',
+                      experLevel=params.LEVEL_ADVANCED,
+                      help='Radius of the sphere where the spherical harmonics will be computed (in voxels).')
+        form.addParam('l1', params.IntParam, default=3,
+                      label='Zernike Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Zernike Polynomials of the deformation=1,2,3,...')
+        form.addParam('l2', params.IntParam, default=2,
+                      label='Harmonical Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Spherical Harmonics of the deformation=1,2,3,...')
+        form.addParam('penalization', params.FloatParam, default=0.00025, label='Regularization',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Penalization to deformations (higher values penalize more the deformation).')
+        form.addParallelSection(threads=1, mpi=1)
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        volList = []
+        dimList = []
+        srList = []
+        volList, dimList, srList, _ = self._iterInputVolumes(self.inputVolumes, volList, dimList, srList, [])
+        if self.twoSets.get():
+            volList, dimList, srList, _ = self._iterInputVolumes(self.secondSet, volList, dimList, srList, [])
+
+        # Resize volumes according to new sampling rate
+        nVoli = 1
+        depsConvert = []
+        for _ in volList:
+            convert = self._insertFunctionStep('convertStep', volList[nVoli - 1],
+                                     dimList[nVoli - 1], srList[nVoli - 1],
+                                     min(dimList), max(srList), nVoli, prerequisites=[])
+            depsConvert.append(convert)
+            nVoli += 1
+
+        # Align all volumes to the first one (reference)
+        nVoli = 1
+        nVolj = 2
+        depsAlign = []
+        for _ in volList[1:]:
+            if nVolj != nVoli:
+                stepID = self._insertFunctionStep('alignStep', volList[nVolj - 1],
+                                        volList[nVoli - 1], nVoli - 1,
+                                        nVolj - 1, prerequisites=depsConvert)
+                depsAlign.append(stepID)
+            nVolj += 1
+
+        # Zernikes3D
+        nVoli = 1
+        depsZernike = []
+        count = 1
+        for _ in volList:
+            nVolj = 1
+            for _ in volList:
+                if nVolj != nVoli:
+                    stepID = self._insertFunctionStep('deformStep', volList[nVoli - 1],
+                                            volList[nVolj - 1], nVoli - 1,
+                                            nVolj - 1, count, prerequisites=depsAlign)
+                    depsZernike.append(stepID)
+                nVolj += 1
+                count += 1
+            nVoli += 1
+
+        self._insertFunctionStep('deformationMatrix', volList, prerequisites=depsZernike)
+        self._insertFunctionStep('gatherResultsStepDef')
+
+        self._insertFunctionStep('correlationMatrix', volList, prerequisites=depsZernike)
+        self._insertFunctionStep('gatherResultsStepCorr')
+
+        self._insertFunctionStep('entropyConsensus')
+
+    # --------------------------- STEPS functions ---------------------------------------------------
+    def convertStep(self, volFn, volDim, volSr, minDim, maxSr, nVoli):
+        Xdim = volDim
+        Ts = volSr
+        newTs = self.targetResolution.get() * 1.0 / 3.0
+        newTs = max(maxSr, newTs)
+        newXdim = Xdim * Ts / newTs
+        newRmax = self.Rmax.get() * Ts / newTs
+        self.newRmax = min(newRmax, self.Rmax.get())
+        fnOut = os.path.splitext(volFn)[0]
+        fnOut = self._getExtraPath(os.path.basename(fnOut + self.OUTPUT_SUFFIX % nVoli))
+
+        ih = ImageHandler()
+        volFn = volFn if getExt(volFn) == '.vol' else volFn + ':mrc'
+        if Xdim != newXdim:
+            self.runJob("xmipp_image_resize",
+                        "-i %s -o %s --dim %d" % (volFn, fnOut, newXdim))
+
+        else:
+            ih.convert(volFn, fnOut)
+
+        if newXdim>minDim:
+            self.runJob("xmipp_transform_window", " -i %s -o %s --crop %d" %
+                        (fnOut, fnOut, (newXdim - minDim)))
+
+    def alignStep(self, inputVolFn, refVolFn, i, j):
+        inputVolFn = self._getExtraPath(os.path.basename(os.path.splitext(inputVolFn)[0]
+                                                         + self.OUTPUT_SUFFIX % (j+1)))
+        refVolFn = self._getExtraPath(os.path.basename(os.path.splitext(refVolFn)[0]
+                                                       + self.OUTPUT_SUFFIX % (i+1)))
+        fnOut = self._getTmpPath(self.ALIGNED_VOL % (j + 1))
+        params = ' --i1 %s --i2 %s --apply %s --local --dontScale' % \
+                 (refVolFn, inputVolFn, fnOut)
+
+        self.runJob("xmipp_volume_align", params)
+
+    def deformStep(self, inputVolFn, refVolFn, i, j, step_id):
+        if j == 0:
+            refVolFn_aux = self._getExtraPath(os.path.basename(os.path.splitext(refVolFn)[0]
+                                                               + self.OUTPUT_SUFFIX % (j + 1)))
+        else:
+            refVolFn_aux = self._getTmpPath(self.ALIGNED_VOL % (j + 1))
+        if i == 0:
+            fnOut_aux = self._getExtraPath(os.path.basename(os.path.splitext(inputVolFn)[0]
+                                                            + self.OUTPUT_SUFFIX % (i + 1)))
+        else:
+            fnOut_aux = self._getTmpPath(self.ALIGNED_VOL % (i + 1))
+        refVolFn = self._getTmpPath("reference_%d.vol" % step_id)
+        fnOut = self._getTmpPath("input_%d.vol" % step_id)
+        copyFile(refVolFn_aux, refVolFn)
+        copyFile(fnOut_aux, fnOut)
+        fnOut2 = self._getTmpPath('vol%dDeformedTo%d.vol' % (i + 1, j + 1))
+
+        params = ' -i %s -r %s -o %s --l1 %d --l2 %d --sigma "%s" --oroot %s --regularization %f' %\
+                 (fnOut, refVolFn, fnOut2, self.l1.get(), self.l2.get(), self.sigma.get(),
+                  self._getExtraPath('Pair_%d_%d' % (i, j)), self.penalization.get())
+        if self.newRmax != 0:
+            params = params + ' --Rmax %d' % self.newRmax
+
+        if self.useGpu.get():
+            self.runJob("xmipp_cuda_volume_deform_sph", params)
+        else:
+            params = params + ' --thr 1'
+            self.runJob("xmipp_volume_deform_sph", params)
+
+        self.computeCorr(fnOut2, refVolFn, i, j)
+
+        cleanPath(fnOut)
+        cleanPath(refVolFn)
+        cleanPath(fnOut2)
+
+    def deformationMatrix(self, volList):
+        cleanPath(self._getTmpPath("*.vol"))
+        numVol = len(volList)
+        self.distanceMatrix = np.zeros((numVol, numVol))
+        for i in range(numVol):
+            for j in range(numVol):
+                if i != j:
+                    path = self._getExtraPath('Pair_%d_%d_deformation.txt' % (i, j))
+                    self.distanceMatrix[i, j] = np.loadtxt(path)
+
+    def correlationMatrix(self, volList):
+        numVol = len(volList)
+        self.corrMatrix = np.zeros((numVol, numVol))
+        for i in range(numVol):
+            for j in range(numVol):
+                if i != j:
+                    path = self._getExtraPath('Pair_%d_%d_correlation.txt' % (i,j))
+                    self.corrMatrix[i, j] = np.loadtxt(path)
+
+
+    def gatherResultsStepDef(self):
+        fnRoot = self._getExtraPath("DistanceMatrix.txt")
+        self.saveDeformation(self.distanceMatrix, fnRoot)
+        if self.twoSets.get():
+            half = int(self.distanceMatrix.shape[0] / 2)
+            subMatrixes = self.split(self.distanceMatrix, half, half)
+            for idm in range(4):
+                fnRoot = self._getExtraPath("DistanceSubMatrix_%d.txt" % (idm + 1))
+                self.saveDeformation(subMatrixes[idm], fnRoot, 'Sub_%d_' % (idm + 1))
+
+    def saveDeformation(self, matrix, fnRoot, label=''):
+        np.savetxt(fnRoot, matrix, "%f")
+        distance = np.asarray(matrix)
+        for i in range(1, 4):
+            embed = mds(distance, i)
+            embedExtended = np.pad(embed, ((0, 0), (0, i - embed.shape[1])),
+                                   "constant", constant_values=0)
+            np.savetxt(self._defineResultsName(i, label), embedExtended)
+
+    def computeCorr(self, vol1, vol2, i, j):
+        vol = xmippLib.Image(vol1)
+        defVol = xmippLib.Image(vol2)
+        corr = vol.correlation(defVol)
+        corr = 1-corr
+        outFile = self._getExtraPath('Pair_%d_%d_correlation.txt' % (i, j))
+        with open(outFile, 'w') as f:
+            f.write('%f' % corr)
+
+    def gatherResultsStepCorr(self):
+        fnRoot = self._getExtraPath("CorrMatrix.txt")
+        self.saveCorrelation(self.corrMatrix, fnRoot)
+        if self.twoSets.get():
+            half = int(self.distanceMatrix.shape[0] / 2)
+            subMatrixes = self.split(self.corrMatrix, half, half)
+            for idm in range(4):
+                fnRoot = self._getExtraPath("CorrSubMatrix_%d.txt" % (idm + 1))
+                self.saveCorrelation(subMatrixes[idm], fnRoot, 'Sub_%d_' % (idm + 1))
+
+    def saveCorrelation(self, matrix, fnRoot, label=''):
+        np.savetxt(fnRoot, matrix, "%f")
+        corr = np.asarray(matrix)
+        for i in range(1, 4):
+            embed = mds(corr, i)
+            embedExtended = np.pad(embed, ((0, 0), (0, i - embed.shape[1])),
+                                   "constant", constant_values=0)
+            np.savetxt(self._defineResultsName2(i, label), embedExtended)
+
+    def rigidAlignmentMapping(self, i):
+        X1 = np.loadtxt(self._defineResultsName(i))
+        X2 = np.loadtxt(self._defineResultsName2(i))
+
+        # Normalize the matrices
+        [_, S, _] = np.linalg.svd(X1)
+        X1 = X1 / np.amax(np.abs(S))
+        [_, S, _] = np.linalg.svd(X2)
+        X2 = X2 / np.amax(np.abs(S))
+
+        # Register the points (taking mirrors into account)
+        for idf in range(4):
+            if idf == 0:
+                R, t = rigidRegistration(X1, X2)
+                X1 = t + np.transpose(R @ X1.T)
+                cost = np.sum(pdist(X1 - X2))
+            elif idf == 1:
+                X1_attempt = np.fliplr(X1)
+                R, t = rigidRegistration(X1_attempt, X2)
+                X1_attempt = t + np.transpose(R @ X1_attempt.T)
+                cost_attempt = np.sum(pdist(X1_attempt - X2))
+                if cost_attempt < cost:
+                    cost = cost_attempt
+                    X1 = X1_attempt
+            elif idf == 2:
+                X1_attempt = np.flipud(X1)
+                R, t = rigidRegistration(X1_attempt, X2)
+                X1_attempt = t + np.transpose(R @ X1_attempt.T)
+                cost_attempt = np.sum(pdist(X1_attempt - X2))
+                if cost_attempt < cost:
+                    cost = cost_attempt
+                    X1 = X1_attempt
+            elif idf == 3:
+                X1_attempt = np.fliplr(X1)
+                X1_attempt = np.flipud(X1_attempt)
+                R, t = rigidRegistration(X1_attempt, X2)
+                X1_attempt = t + np.transpose(R @ X1_attempt.T)
+                cost_attempt = np.sum(pdist(X1_attempt - X2))
+                if cost_attempt < cost:
+                    cost = cost_attempt
+                    X1 = X1_attempt
+        return X1, X2
+
+    def gaussianKernel(self, sigma, i):
+        lbox = int(6 * sigma)
+        if lbox % 2 == 0:
+            lbox += 1
+        mid = int((lbox - 1) / 2 + 1)
+        if i == 2:
+            kernel = np.zeros((lbox, lbox))
+            kernel[mid, mid] = 1
+            kernel = gaussian_filter(kernel, sigma=sigma)
+        else:
+            kernel = np.zeros((lbox, lbox, lbox))
+            kernel[mid, mid, mid] = 1
+            kernel = gaussian_filter(kernel, sigma=sigma)
+        return kernel, lbox, mid
+
+    def convolution(self, Xr, kernel, sigma, R, C, D):
+        # Create the grid
+        lbox = int(6 * sigma)
+        if lbox % 2 == 0:
+            lbox += 1
+        mid = int((lbox - 1) / 2 + 1)
+        S_shape = [d + lbox for d in R.shape]
+        S = np.zeros(S_shape)
+
+        # Place rounded points on the grid
+        for p in range(Xr.shape[0]):
+            if D is None:
+                indy = np.argmin(np.abs(R[:, 0] - Xr[p, 0]))
+                indx = np.argmin(np.abs(C[0, :] - Xr[p, 1]))
+                S[indx:indx + 2 * mid - 1, indy:indy + 2 * mid - 1] += kernel
+            else:
+                indx = np.argmin(np.abs(R[:, 0, 0] - Xr[p, 0]))
+                indy = np.argmin(np.abs(C[0, :, 0] - Xr[p, 1]))
+                indz = np.argmin(np.abs(D[0, 0, :] - Xr[p, 2]))
+                S[indx:indx + 2 * mid - 1, indy:indy + 2 * mid - 1, indz:indz + 2 * mid - 1] += kernel
+        return S
+
+    def entropyConsensus(self):
+        for i in range(2, 4):
+            # Rigid alignment of mappings (considering mirrors as well)
+            X1, X2 = self.rigidAlignmentMapping(i)
+
+            # Round point to place them in a grid
+            Xr1 = np.round(X1, decimals=3)
+            Xr2 = np.round(X2, decimals=3)
+            size_grid = 2.5 * max((np.amax(Xr1), np.amax(Xr2)))
+
+            # Parameters needed for future convolution
+            grid_coords = np.linspace(-size_grid, size_grid, num=400)
+            if i == 2:
+                R, C = np.meshgrid(grid_coords, grid_coords, indexing='ij')
+                D = None
+            else:
+                R, C, D = np.meshgrid(grid_coords, grid_coords, grid_coords, indexing='ij')
+            sigma = R.shape[0] / 20
+
+            # Create Gaussian Kernel
+            kernel, lbox, mid = self.gaussianKernel(sigma, i)
+
+            # Consensus
+            alpha_vect = np.arange(0, 1.01, 0.01)
+            entropy_vect = []
+            X_matrices = []
+
+            for alpha in alpha_vect:
+                X = alpha * X1 + (1 - alpha) * X2
+                X_matrices.append(X)
+
+                # Round points to place them in a grid
+                Xr = np.round(X, decimals=3)
+
+                # Convolve grid and kernel
+                S = self.convolution(Xr, kernel, sigma, R, C, D)
+
+                # Compute the Shannon entropy associated to the convolved grid
+                _, counts = np.unique(S, return_counts=True)
+                entropy_vect.append(entropy(counts, base=2))
+
+            # Find optimal entropy value (minimum)
+            entropy_vect = np.asarray(entropy_vect)
+            id_peaks, _ = find_peaks(-entropy_vect)
+            if id_peaks.size > 1:
+                peaks = entropy_vect[id_peaks]
+                id_optimal = id_peaks[np.argmin(peaks)]
+            elif id_peaks.size == 1:
+                id_optimal = id_peaks[0]
+            else:
+                id_optimal = 0
+
+            X_optimal = X_matrices[id_optimal]
+
+            np.savetxt(self._defineResultsName3(i), X_optimal)
+
+    # --------------------------- UTILS functions --------------------------------------------
+    def _iterInputVolumes(self, volumes, volList, dimList, srList, idList):
+        """ Iterate over all the input volumes. """
+        count = 1
+        for pointer in volumes:
+            item = pointer.get()
+            if item is None:
+                break
+            if isinstance(item, Volume):
+                volList.append(item.getFileName())
+                dimList.append(item.getDim()[0])
+                srList.append(item.getSamplingRate())
+                idList.append(count)
+            elif isinstance(item, SetOfVolumes):
+                for vol in item:
+                    volList.append(vol.getFileName())
+                    dimList.append(vol.getDim()[0])
+                    srList.append(vol.getSamplingRate())
+                    idList.append(count)
+                    count += 1
+            count += 1
+        return volList, dimList, srList, idList
+
+    def split(self, array, nrows, ncols):
+        """Split a matrix into sub-matrices."""
+        r, h = array.shape
+        return (array.reshape(h // nrows, nrows, -1, ncols)
+                .swapaxes(1, 2)
+                .reshape(-1, nrows, ncols))
+
+    def _defineResultsName(self, i, label=''):
+        return self._getExtraPath('Coordinate%sMatrix%d.txt' % (label, i))
+
+    def _defineResultsName2(self, i, label=''):
+        return self._getExtraPath('Coordinate%sMatrixCorr%d.txt' % (label, i))
+
+    def _defineResultsName3(self, i):
+        return self._getExtraPath('ConsensusMatrix%d.txt' % i)
+
+    # ------------------------- VALIDATE functions -----------------------------
+    def validate(self):
+        """ Try to find errors on define params. """
+        errors = []
+        l1 = self.l1.get()
+        l2 = self.l2.get()
+        if (l1 - l2) < 0:
+            errors.append('Zernike degree must be higher than '
+                          'SPH degree.')
+        return errors
+
+
+
+
diff --git a/xmipp3/protocols/protocol_volume_adjust_sub.py b/xmipp3/protocols/protocol_volume_adjust_sub.py
index 9ae89cf4c..01aa2744c 100644
--- a/xmipp3/protocols/protocol_volume_adjust_sub.py
+++ b/xmipp3/protocols/protocol_volume_adjust_sub.py
@@ -35,6 +35,9 @@
 from pwem.protocols import EMProtocol
 
 
+CITE = 'Fernandez-Gimenez2021'
+
+
 class XmippProtVolAdjBase(EMProtocol):
     """ Helper class that contains some Protocol utilities methods
     used by both  XmippProtVolSubtraction and XmippProtVolAdjust."""
@@ -60,6 +63,10 @@ def _defineParams(self, form):
                       expertLevel=LEVEL_ADVANCED,
                       help='Relaxation factor for Fourier amplitude projector (POCS), it should be between 0 and 1, '
                            'being 1 no relaxation and 0 no modification of volume 2 amplitudes')
+        form.addParam('radavg', BooleanParam, label="Match the rotationally averaged Fourier amplitudes?", default=True,
+                      help='Match the rotationally averaged Fourier amplitudes when adjusting the amplitudes instead of'
+                           ' taking them directly from the reference volume. For subtraction and consensus it is '
+                           'recommended to set it to True but for sharpening it is recommended to set it to False')
         form.addParam('computeE', BooleanParam, label="Compute energy?", default=True, expertLevel=LEVEL_ADVANCED,
                       help='Compute energy difference between the different adjustment steps and iterations to see if '
                            'the method reaches convergence')
@@ -198,6 +205,8 @@ def subtractionStep(self):
         if self.saveFiles:
             args += ' --saveV1 %s --saveV2 %s' % (self._getExtraPath('vol1_filtered.mrc'),
                                                   self._getExtraPath('vol2_adjusted.mrc'))
+        if self.radavg:
+            args += ' --radavg'
         if self.computeE:
             args += ' --computeEnergy'
         self.runJob(program, args)
@@ -218,6 +227,8 @@ def _summary(self):
             summary.append("Input mask 2: %s" % self.mask2.get().getFileName())
         if self.resol.get() != 0:
             summary.append("Subtraction at resolution %f A" % self.resol.get())
+        if self.radavg:
+            summary.append("Matching the rotational averaged Fourier amplitudes")
         return summary
 
     def _methods(self):
@@ -238,6 +249,9 @@ def _validate(self):
             errors.append('Relaxation factor (lambda) must be between 0 and 1')
         return errors
 
+    def _citations(self):
+        return [CITE]
+
     # --------------------------- UTLIS functions --------------------------------------------
     def _getPdbFileName(self):
         if self.inputPdbData == self.IMPORT_OBJ:
@@ -283,6 +297,8 @@ def adjustStep(self):
             args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())
         if self.masks:
             args += ' --mask1 %s --mask2 %s' % (self.mask1.get().getFileName(), self.mask2.get().getFileName())
+        if self.radavg:
+            args += ' --radavg'
         if self.computeE:
             args += ' --computeEnergy'
         self.runJob(program, args)
@@ -295,6 +311,8 @@ def _summary(self):
             summary.append("Input mask 2: %s" % self.mask2.get().getFileName())
         if self.resol.get() != 0:
             summary.append("Filter at resolution %f A" % self.resol.get())
+        if self.radavg:
+            summary.append("Matching the rotational averaged Fourier amplitudes")
         return summary
 
     def _methods(self):
@@ -314,3 +332,6 @@ def _validate(self):
         if rfactor < 0 or rfactor > 1:
             errors.append('Relaxation factor (lambda) must be between 0 and 1')
         return errors
+
+    def _citations(self):
+        return [CITE]
diff --git a/xmipp3/protocols/protocol_volume_consensus.py b/xmipp3/protocols/protocol_volume_consensus.py
index 0fe97acf5..62d457413 100644
--- a/xmipp3/protocols/protocol_volume_consensus.py
+++ b/xmipp3/protocols/protocol_volume_consensus.py
@@ -103,3 +103,6 @@ def _methods(self):
             methods.append("We compute a consensus volume from %d input volumes at %f A/px" %
                            (self.vols.getSize(), self.vols[0].get().getSamplingRate()))
         return methods
+
+    def _citations(self):
+        return ['Fernandez-Gimenez2021']
diff --git a/xmipp3/protocols/protocol_volume_deform_sph.py b/xmipp3/protocols/protocol_volume_deform_sph.py
deleted file mode 100644
index 09a0438d7..000000000
--- a/xmipp3/protocols/protocol_volume_deform_sph.py
+++ /dev/null
@@ -1,136 +0,0 @@
-########
-
-# This protocol has been removed from Scipion as there is a newer version being developed that
-# will solve some of the problems of this current implementation. The code will be available
-# again in a newer release.
-
-########
-
-
-
-# **************************************************************************
-# *
-# * Authors:     Amaya Jimenez Moreno (ajimenez@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 pwem.protocols import ProtAnalysis3D
-# import pyworkflow.protocol.params as params
-# from pwem.emlib.image import ImageHandler
-# from pwem.objects import Volume
-# from pyworkflow import VERSION_2_0
-#
-#
-# class XmippProtVolumeDeformSPH(ProtAnalysis3D):
-#     """ Protocol for volume deformation based on spherical harmonics. """
-#     _label = 'sph volume deform'
-#     _lastUpdateVersion = VERSION_2_0
-#
-#     # --------------------------- DEFINE param functions --------------------------------------------
-#     def _defineParams(self, form):
-#         form.addSection(label='Input')
-#         form.addParam('inputVolume', params.PointerParam, label="Input volume",
-#                       pointerClass='Volume')
-#         form.addParam('refVolume', params.PointerParam, label="Reference volume",
-#                       pointerClass='Volume')
-#         form.addParam('targetResolution', params.FloatParam, label="Target resolution",
-#                       default=8.0,
-#                       help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
-#                            "2/3 of the Fourier spectrum.")
-#         form.addParam('depth', params.IntParam, default=3,
-#                       label='Harmonical depth',
-#                       expertLevel=params.LEVEL_ADVANCED,
-#                       help='Harmonical depth of the deformation=1,2,3,...')
-#
-#
-#     def _createFilenameTemplates(self):
-#         """ Centralize how files are called """
-#         myDict = {
-#             'fnRefVol': self._getExtraPath('ref_volume.vol'),
-#             'fnInputVol': self._getExtraPath('input_volume.vol'),
-#             'fnOutVol': self._getExtraPath('vol1DeformedTo2.vol')
-#                  }
-#         self._updateFilenamesDict(myDict)
-#
-#     # --------------------------- INSERT steps functions -----------------------
-#     def _insertAllSteps(self):
-#         self._createFilenameTemplates()
-#         self._insertFunctionStep("convertStep")
-#         self._insertFunctionStep("deformStep")
-#         self._insertFunctionStep("createOutputStep")
-#
-#     # --------------------------- STEPS functions ------------------------------
-#     def convertStep(self):
-#         fnInputVol = self._getFileName('fnInputVol')
-#         fnRefVol = self._getFileName('fnRefVol')
-#
-#         XdimI = self.inputVolume.get().getDim()[0]
-#         TsI = self.inputVolume.get().getSamplingRate()
-#         XdimR = self.refVolume.get().getDim()[0]
-#         TsR = self.refVolume.get().getSamplingRate()
-#         self.newTs = self.targetResolution.get() * 1.0 / 3.0
-#         self.newTs = max(TsI, TsR, self.newTs)
-#         newXdimI = int(XdimI * TsI / self.newTs)
-#         newXdimR = int(XdimR * TsR / self.newTs)
-#         self.newXdim = min(newXdimI, newXdimR)
-#
-#         ih = ImageHandler()
-#         ih.convert(self.inputVolume.get(), fnInputVol)
-#         if XdimI != newXdimI:
-#             self.runJob("xmipp_image_resize",
-#                         "-i %s --dim %d" % (fnInputVol, newXdimI))
-#         if newXdimI>self.newXdim:
-#             self.runJob("xmipp_transform_window", " -i %s --crop %d" %
-#                         (fnInputVol, (newXdimI-self.newXdim)))
-#
-#
-#         ih.convert(self.refVolume.get(), fnRefVol)
-#         if XdimR != newXdimR:
-#             self.runJob("xmipp_image_resize",
-#                         "-i %s --dim %d" % (fnRefVol, newXdimR))
-#         if newXdimR>self.newXdim:
-#             self.runJob("xmipp_transform_window", " -i %s --crop %d" %
-#                         (fnRefVol, (newXdimR-self.newXdim)))
-#
-#
-#     def deformStep(self):
-#         fnInputVol = self._getFileName('fnInputVol')
-#         fnRefVol = self._getFileName('fnRefVol')
-#         fnOutVol = self._getFileName('fnOutVol')
-#
-#         params = ' --i1 %s --i2 %s --apply %s --least_squares --local ' % \
-#                  (fnRefVol, fnInputVol, fnOutVol)
-#         self.runJob("xmipp_volume_align", params)
-#
-#
-#         params = ' -i %s -r %s -o %s --analyzeStrain --depth %d ' % \
-#                  (fnOutVol, fnRefVol, fnOutVol, self.depth.get())
-#         self.runJob("xmipp_volume_deform_sph", params)
-#
-#
-#     def createOutputStep(self):
-#         vol = Volume()
-#         vol.setLocation(self._getFileName('fnOutVol'))
-#         vol.setSamplingRate(self.newTs)
-#         self._defineOutputs(outputVolume=vol)
-#         self._defineSourceRelation(self.inputVolume, vol)
-#
diff --git a/xmipp3/protocols/protocol_volume_deform_zernike3d.py b/xmipp3/protocols/protocol_volume_deform_zernike3d.py
new file mode 100644
index 000000000..306607bcd
--- /dev/null
+++ b/xmipp3/protocols/protocol_volume_deform_zernike3d.py
@@ -0,0 +1,211 @@
+
+# **************************************************************************
+# *
+# * Authors:     Amaya Jimenez Moreno (ajimenez@cnb.csic.es)
+# *              David Herreros Calero (dherreros@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 pwem.protocols import ProtAnalysis3D
+import pyworkflow.protocol.params as params
+from pwem.emlib.image import ImageHandler
+from pwem.objects import Volume
+from pyworkflow import VERSION_2_0
+
+
+class XmippProtVolumeDeformZernike3D(ProtAnalysis3D):
+    """ Protocol for volume deformation based on Zernike3D. """
+    _label = 'volume deform - Zernike3D'
+    _lastUpdateVersion = VERSION_2_0
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addHidden(params.USE_GPU, params.BooleanParam, default=True,
+                       label="Use GPU for execution",
+                       help="This protocol has both CPU and GPU implementation.\
+                           Select the one you want to use.")
+        form.addHidden(params.GPU_LIST, params.StringParam, default='0',
+                       expertLevel=params.LEVEL_ADVANCED,
+                       label="Choose GPU IDs",
+                       help="Add a list of GPU devices that can be used")
+        form.addParam('refVolume', params.PointerParam, label="Volume 1",
+                      pointerClass='Volume')
+        form.addParam('inputVolume', params.PointerParam, label="Volume 2",
+                      pointerClass='Volume')
+        form.addParam('sigma', params.NumericListParam, label="Multiresolution", default="1 2",
+                      help="Perform the analysys comparing different filtered versions of the volumes")
+        form.addParam('targetResolution', params.FloatParam, label="Target resolution",
+                      default=8.0,
+                      help="In Angstroms, the images and the volume are rescaled so that this resolution is at "
+                           "2/3 of the Fourier spectrum.")
+        form.addParam('Rmax', params.IntParam, default=0,
+                      label='Sphere radius',
+                      experLevel=params.LEVEL_ADVANCED,
+                      help='Radius of the sphere where the spherical harmonics will be computed.')
+        form.addParam('l1', params.IntParam, default=3,
+                      label='Zernike Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Zernike Polynomials of the deformation=1,2,3,...')
+        form.addParam('l2', params.IntParam, default=2,
+                      label='Harmonical Degree',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Degree Spherical Harmonics of the deformation=1,2,3,...')
+        form.addParam('penalization', params.FloatParam, default=0.00025, label='Regularization',
+                      expertLevel=params.LEVEL_ADVANCED,
+                      help='Penalization to deformations (higher values penalize more the deformation).')
+        form.addParallelSection(threads=4, mpi=0)
+
+
+    def _createFilenameTemplates(self):
+        """ Centralize how files are called """
+        myDict = {
+            'fnRefVol': self._getExtraPath('ref_volume.vol'),
+            'fnInputVol': self._getExtraPath('input_volume.vol'),
+            'fnInputFilt': self._getExtraPath('input_volume_filt.vol'),
+            'fnRefFilt': self._getExtraPath('ref_volume_filt.vol'),
+            'fnOutVol': self._getExtraPath('vol1DeformedTo2.vol')
+                 }
+        self._updateFilenamesDict(myDict)
+
+    # --------------------------- INSERT steps functions -----------------------
+    def _insertAllSteps(self):
+        self._createFilenameTemplates()
+        self._insertFunctionStep("convertStep")
+        self._insertFunctionStep("deformStep")
+        self._insertFunctionStep("createOutputStep")
+
+    # --------------------------- STEPS functions ------------------------------
+    def convertStep(self):
+        fnInputVol = self._getFileName('fnInputVol')
+        fnRefVol = self._getFileName('fnRefVol')
+
+        XdimI = self.inputVolume.get().getDim()[0]
+        TsI = self.inputVolume.get().getSamplingRate()
+        XdimR = self.refVolume.get().getDim()[0]
+        TsR = self.refVolume.get().getSamplingRate()
+        self.newTs = self.targetResolution.get() * 1.0 / 3.0
+        self.newTs = max(TsI, TsR, self.newTs)
+        newXdimI = XdimI * TsI / self.newTs
+        newXdimR = XdimR * TsR / self.newTs
+        newRmax = self.Rmax.get() * TsI / self.newTs
+        self.newXdim = min(newXdimI, newXdimR)
+        self.newRmax = min(newRmax, self.Rmax.get())
+
+        ih = ImageHandler()
+        ih.convert(self.inputVolume.get(), fnInputVol)
+        if XdimI != newXdimI:
+            self.runJob("xmipp_image_resize",
+                        "-i %s --dim %d" % (fnInputVol, newXdimI))
+        if newXdimI > self.newXdim:
+            self.runJob("xmipp_transform_window", " -i %s --crop %d" %
+                        (fnInputVol, (newXdimI-self.newXdim)))
+
+
+        ih.convert(self.refVolume.get(), fnRefVol)
+
+        if XdimR != newXdimR:
+            self.runJob("xmipp_image_resize",
+                        "-i %s --dim %d" % (fnRefVol, newXdimR))
+        if newXdimR > self.newXdim:
+            self.runJob("xmipp_transform_window", " -i %s --crop %d" %
+                        (fnRefVol, (newXdimR-self.newXdim)))
+
+
+    def deformStep(self):
+        fnRefVol = self._getFileName('fnRefVol')
+        fnOutVol = self._getFileName('fnOutVol')
+
+        self.alignMaps()
+
+        params = ' -i %s -r %s -o %s --analyzeStrain --l1 %d --l2 %d --sigma "%s" --oroot %s --regularization %f' % \
+                 (fnOutVol, fnRefVol, fnOutVol, self.l1.get(), self.l2.get(), self.sigma.get(),
+                  self._getExtraPath('Volumes'), self.penalization.get())
+        if self.newRmax != 0:
+            params = params + ' --Rmax %d' % self.newRmax
+
+        if self.useGpu.get():
+            self.runJob("xmipp_cuda_volume_deform_sph", params)
+        else:
+            if self.numberOfThreads.get() != 0:
+                params = params + ' --thr %d' % self.numberOfThreads.get()
+            self.runJob("xmipp_volume_deform_sph", params)
+
+
+    def createOutputStep(self):
+        if self.targetResolution.get() != 3.0:
+            correctionFactor = self.targetResolution.get() / 3.0
+            with open(self._getExtraPath('Volumes_clnm.txt'), 'r') as fid:
+                lines = fid.readlines()
+                basisParams = np.fromstring(lines[0], sep=' ')
+                if self.Rmax.get() != 0:
+                    basisParams[2] = self.Rmax.get()
+                else:
+                    basisParams[2] = self.refVolume.get().getDim()[0] / 2
+                clnm = np.fromstring(lines[1], sep=' ') * correctionFactor
+            with open(self._getExtraPath('Volumes_clnm.txt'), 'w') as fid:
+                fid.write(' '.join(map(str, basisParams)) + "\n")
+                fid.write(' '.join(map(str, clnm)) + "\n")
+        vol = Volume()
+        vol.setLocation(self._getFileName('fnOutVol'))
+        vol.setSamplingRate(self.newTs)
+        self._defineOutputs(outputVolume=vol)
+        self._defineSourceRelation(self.inputVolume, vol)
+
+    def alignMaps(self):
+        fnInputVol = self._getFileName('fnInputVol')
+        fnInputFilt = self._getFileName('fnInputFilt')
+        fnRefVol = self._getFileName('fnRefVol')
+        fnRefFilt = self._getFileName('fnRefFilt')
+        fnOutVol = self._getFileName('fnOutVol')
+
+        # Filter the volumes to improve alignment quality
+        params = " -i %s -o %s --fourier real_gaussian 2" % (fnInputVol, fnInputFilt)
+        self.runJob("xmipp_transform_filter", params)
+        params = " -i %s -o %s --fourier real_gaussian 2" % (fnRefVol, fnRefFilt)
+        self.runJob("xmipp_transform_filter", params)
+
+        # Find transformation needed to align the volumes
+        params = ' --i1 %s --i2 %s --local --dontScale ' \
+                 '--copyGeo %s' % \
+                 (fnRefFilt, fnInputFilt, self._getExtraPath("geo.txt"))
+        self.runJob("xmipp_volume_align", params)
+
+        # Apply transformation of filtered volume to original volume
+        with open(self._getExtraPath("geo.txt"), 'r') as file:
+            geo_str = file.read().replace('\n', ',')
+        params = " -i %s -o %s --matrix %s" % (fnInputVol, fnOutVol, geo_str)
+        self.runJob("xmipp_transform_geometry", params)
+
+    # ------------------------- VALIDATE functions -----------------------------
+    def validate(self):
+        """ Try to find errors on define params. """
+        errors = []
+        l1 = self.l1.get()
+        l2 = self.l2.get()
+        if (l1 - l2) < 0:
+            errors.append('Zernike degree must be higher than '
+                          'SPH degree.')
+        return errors
\ No newline at end of file
diff --git a/xmipp3/viewers/__init__.py b/xmipp3/viewers/__init__.py
index a8586eff5..367bc596e 100644
--- a/xmipp3/viewers/__init__.py
+++ b/xmipp3/viewers/__init__.py
@@ -58,12 +58,12 @@
 
 from .viewer_projmatch import XmippProjMatchViewer
 
-#AJ
-# from .viewer_structure_map_sph import XmippProtStructureMapSphViewer
-# from .viewer_volume_deform_sph import XmippVolumeDeformSphViewer
+from .viewer_volume_deform_zernike3d import XmippVolumeDeformZernike3DViewer
+from .viewer_structure_map import XmippProtStructureMapViewer
+from .viewer_pdb_deform_zernike3d import XmippPDBDeformViewer
 from .viewer_metaprotocol_golden_highres import XmippMetaprotocolGoldenHighResViewer
 from .viewer_deep_align import XmippDeepAlignViewer
 
 # TODO: Import from continuousflex the modules needed to create clusters (soft dependency)
-# from .viewer_angular_alignment_sph import XmippAngularAlignmentSphViewer
+from .viewer_angular_alignment_sph import XmippAngularAlignmentSphViewer
 
diff --git a/xmipp3/viewers/viewer_angular_alignment_sph.py b/xmipp3/viewers/viewer_angular_alignment_sph.py
index 1be7b5d57..6fc605c87 100644
--- a/xmipp3/viewers/viewer_angular_alignment_sph.py
+++ b/xmipp3/viewers/viewer_angular_alignment_sph.py
@@ -32,13 +32,13 @@
 import matplotlib.pyplot as plt
 from pyworkflow.utils.path import cleanPath
 from pwem.objects import SetOfParticles
-from xmipp3.protocols.protocol_angular_alignment_sph import XmippProtAngularAlignmentSPH
+from xmipp3.protocols.protocol_angular_alignment_zernike3d import XmippProtAngularAlignmentZernike3D
 
 
 class XmippAngularAlignmentSphViewer(ProtocolViewer):
     """ Visualize the output of protocol volume strain """
     _label = 'viewer angular align sph'
-    _targets = [XmippProtAngularAlignmentSPH]
+    _targets = [XmippProtAngularAlignmentZernike3D]
     _environments = [DESKTOP_TKINTER, WEB_DJANGO]
 
     def __init__(self, **kwargs):
diff --git a/xmipp3/viewers/viewer_pdb_deform_zernike3d.py b/xmipp3/viewers/viewer_pdb_deform_zernike3d.py
new file mode 100644
index 000000000..6834be13a
--- /dev/null
+++ b/xmipp3/viewers/viewer_pdb_deform_zernike3d.py
@@ -0,0 +1,89 @@
+# **************************************************************************
+# *
+# * Authors:  David Herreros Calero (dherreros@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 os
+
+import pyworkflow.viewer as pwviewer
+import pwem.viewers.views as vi
+import pyworkflow.protocol.params as params
+from pwem.viewers import ChimeraView
+from xmipp3.protocols.protocol_pdb_deform_zernike3d import XmippProtDeformPDBZernike3D
+
+
+class XmippPDBDeformViewer(pwviewer.ProtocolViewer):
+    """ Visualize the deformation applied to the PDB file """
+    _label = 'viewer pdb deform sph'
+    _targets = [XmippProtDeformPDBZernike3D]
+    _environments = [pwviewer.DESKTOP_TKINTER, pwviewer.WEB_DJANGO]
+    OPEN_FILE = "open %s\n"
+
+    def _defineParams(self, form):
+        form.addSection(label='Show deformation')
+        form.addParam('doShowPDB', params.LabelParam,
+                      label="Display original and deformed PDB")
+        form.addParam('doShowMorph', params.LabelParam,
+                      label="Display a morphing between the original and deformed PDB")
+
+    def _getVisualizeDict(self):
+        return {'doShowPDB': self._doShowPDB,
+                'doShowMorph': self._doShowMorph}
+
+    def _doShowPDB(self, obj, **kwargs):
+        scriptFile = self.protocol._getPath('pdb_deform_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        inputFile = os.path.abspath(self.protocol.inputPDB.get().getFileName())
+        outputFile = os.path.abspath(self.protocol.outputPDB.getFileName())
+
+        fhCmd.write(self.OPEN_FILE % inputFile)
+        fhCmd.write(self.OPEN_FILE % outputFile)
+        # fhCmd.write("start Model Panel\n")
+        fhCmd.write("show cartoons\n")
+        fhCmd.write("cartoon style width 1.5 thick 1.5\n")
+        fhCmd.write("style stick\n")
+        fhCmd.write("color bymodel\n")
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
+
+    def _doShowMorph(self, obj, **kwargs):
+        scriptFile = self.protocol._getPath('pdb_deform_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        inputFile = os.path.abspath(self.protocol.inputPDB.get().getFileName())
+        outputFile = os.path.abspath(self.protocol.outputPDB.getFileName())
+
+        fhCmd.write(self.OPEN_FILE % inputFile)
+        fhCmd.write(self.OPEN_FILE % outputFile)
+        fhCmd.write("hide models\n")
+        fhCmd.write("morph #1,2 frames 50 play false\n")
+        fhCmd.write("coordset #3 1,\n")
+        fhCmd.write("wait 50\n")
+        fhCmd.write("coordset #3 50,1\n")
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
+
diff --git a/xmipp3/viewers/viewer_structure_map.py b/xmipp3/viewers/viewer_structure_map.py
new file mode 100644
index 000000000..b393be0b0
--- /dev/null
+++ b/xmipp3/viewers/viewer_structure_map.py
@@ -0,0 +1,359 @@
+# **************************************************************************
+# *
+# * Authors:     David Herreros Calero (dherreros@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 os, matplotlib, math
+from scipy import ndimage
+from scipy.spatial import KDTree
+import tkinter as tk
+matplotlib.use("TkAgg")
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from matplotlib.widgets import RadioButtons
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2Tk
+from mpl_toolkits.mplot3d.proj3d import proj_transform
+from matplotlib.text import Annotation
+
+import numpy as np
+from scipy.ndimage.filters import gaussian_filter
+
+from pyworkflow.viewer import DESKTOP_TKINTER, WEB_DJANGO, ProtocolViewer
+from pyworkflow.gui.plotter import plt
+import pyworkflow.protocol.params as params
+
+from xmipp3.protocols.protocol_structure_map import XmippProtStructureMap
+from xmipp3.protocols.protocol_structure_map_zernike3d import XmippProtStructureMapZernike3D
+
+
+class XmippProtStructureMapViewer(ProtocolViewer):
+    """ Wrapper to visualize different type of data objects
+    with the Xmipp program xmipp_showj
+    """
+    
+    _label = 'viewer structure map'
+    _environments = [DESKTOP_TKINTER, WEB_DJANGO]
+    _targets = [XmippProtStructureMap, XmippProtStructureMapZernike3D]
+
+    def _defineParams(self, form):
+        form.addSection(label='Show StructMap')
+        if isinstance(self.protocol, XmippProtStructureMapZernike3D):
+            form.addParam('map', params.EnumParam, choices=['Deformation', 'Correlation', 'Consensus'],
+                          default=0, help='Choose the type of metric to display the coordinates.')
+            form.addParam('twoSets', params.BooleanParam, default=False, label='Two set analysis',
+                          help='Activate if the analysis of the deformation included the option of analysing '
+                               'two sets of volumes independently.')
+        form.addParam('doShowPlot', params.LabelParam,
+                      label="Display the StructMap")
+    
+    def _getVisualizeDict(self):
+        return {'doShowPlot': self._visualize}
+
+    def getOutputFile(self):
+        fnOutput = ['']
+        if isinstance(self.protocol, XmippProtStructureMapZernike3D):
+            if self.map.get() == 0 and not self.twoSets:
+                fnOutput = [self.protocol._defineResultsName(3)]
+            elif self.map.get() == 1 and not self.twoSets:
+                fnOutput = [self.protocol._defineResultsName2(3)]
+            elif self.map.get() == 0 and self.twoSets:
+                fnOutput = [self.protocol._defineResultsName(3, 'Sub_1_'),
+                            self.protocol._defineResultsName(3, 'Sub_2_')]
+            elif self.map.get() == 1 and self.twoSets:
+                fnOutput = [self.protocol._defineResultsName2(3, 'Sub_1_'),
+                            self.protocol._defineResultsName2(3, 'Sub_2_')]
+            elif self.map.get() == 2 and not self.twoSets:
+                fnOutput = [self.protocol._defineResultsName3(3)]
+        else:
+            fnOutput = [self.protocol._defineResultsName(3)]
+        return fnOutput
+
+    def _visualize(self, e=None):
+        fnOutput = self.getOutputFile()
+        for file in fnOutput:
+            if not os.path.exists(file):
+                return [self.errorMessage('The necessary metadata was not produced\n'
+                                          'Execute again the protocol\n',
+                                          title='Missing result file')]
+
+        self.coordinates = (np.loadtxt(file) for file in fnOutput)
+        self.coordinates = np.vstack(self.coordinates)
+        labels = [str(idp) for idp in range(1, self.coordinates.shape[0] + 1)]
+        if os.path.isfile(self._getExtraPath('weigths.txt')):
+            weights = np.loadtxt(self._getExtraPath('weigths.txt'))
+        else:
+            weights = None
+        plot = projectionPlot(self.coordinates, labels, weights)
+        plot.initializePlot()
+        return plot
+        
+    def _validate(self):
+        errors = []
+        return errors
+
+
+class projectionPlot(object):
+
+    def __init__(self, coords, labels, weights):
+        self.coords = coords
+        self.labels = labels
+        self.weights = weights
+        try:
+            self.minimum_spanning_tree()
+        except IndexError:
+            self.T = None
+        self.proj_coords = None
+        self.radio = None
+        self.cb = None
+        self.prevlabel = 'Scatter'
+        self.root = tk.Tk()
+        self.fig = plt.Figure(figsize=(10, 4))
+        self.canvas = FigureCanvasTkAgg(self.fig, master=self.root)
+        self.ax_3d = self.fig.add_subplot(121, projection="3d")
+        self.ax_3d.set_title("3D view")
+        self.ax_3d.set_axis_off()
+        self.fig.canvas.mpl_connect('button_release_event', self.onRelease)
+        self.ax_2d = self.fig.add_subplot(122)
+        self.ax_2d.set_title("Projection from current 3D view")
+
+    def minimum_spanning_tree(self):
+        N = self.coords.shape[0]
+        tree = KDTree(self.coords)
+        distances, indices = tree.query(self.coords, k=10)
+        nn_mat = np.zeros((N, N))
+        for idn in range(N):
+            nn_mat[idn, indices[idn]] += distances[idn].reshape(-1)
+            nn_mat[idn, idn] = 0
+
+        edges = ((int(e[0]), int(e[1])) for e in zip(*np.asarray(nn_mat).nonzero()))
+        triples = list(((u, v, float(nn_mat[u, v])) for u, v in edges))
+        edges = [(triple[0], triple[1]) for triple in triples]
+        weigths_edge = [triple[2] for triple in triples]
+
+        edge_matrix = np.zeros((N, N))
+        for edge, weigth in zip(edges, weigths_edge):
+            edge_matrix[edge] = weigth
+        self.T = self.KruskalMST(triples, N)
+
+    def KruskalMST(self, triples, N):
+        result = []
+        i = 0
+        e = 0
+        graph = sorted(triples, key=lambda item: item[2])
+        parent = []
+        rank = []
+        for node in range(N):
+            parent.append(node)
+            rank.append(0)
+        while e < N - 1:
+            u, v, w = graph[i]
+            i = i + 1
+            x = self.find(parent, u)
+            y = self.find(parent, v)
+            if x != y:
+                e = e + 1
+                result.append([u, v, w])
+                self.union(parent, rank, x, y)
+        return result
+
+    def find(self, parent, i):
+        if parent[i] == i:
+            return i
+        return self.find(parent, parent[i])
+
+    def union(self, parent, rank, x, y):
+        xroot = self.find(parent, x)
+        yroot = self.find(parent, y)
+        if rank[xroot] < rank[yroot]:
+            parent[xroot] = yroot
+        elif rank[xroot] > rank[yroot]:
+            parent[yroot] = xroot
+        else:
+            parent[yroot] = xroot
+            rank[xroot] += 1
+
+    def mst_3D(self):
+        N = self.coords.shape[0]
+        degree = [0] * N
+        edges_mst = [0] * (N - 1)
+        if self.T is not None:
+            for idn in range(N - 1):
+                degree[self.T[idn][0]] += 1
+                degree[self.T[idn][1]] += 1
+                edges_mst[idn] = (self.T[idn][0], self.T[idn][1])
+            edge_max = max(degree)
+            colors = [plt.cm.plasma(val / edge_max) for val in degree]
+            for edge in edges_mst:
+                if edge != 0:
+                    x = np.array((self.coords[edge[0]][0], self.coords[edge[1]][0]))
+                    y = np.array((self.coords[edge[0]][1], self.coords[edge[1]][1]))
+                    z = np.array((self.coords[edge[0]][2], self.coords[edge[1]][2]))
+                    self.ax_3d.plot(x, y, z, c='black', alpha=0.5)
+        else:
+            colors = [plt.cm.plasma(1) for _ in self.coords]
+        for idn, row in enumerate(self.coords):
+            xi = row[0]
+            yi = row[1]
+            zi = row[2]
+            self.ax_3d.scatter(xi, yi, zi, c=[colors[idn]], s=20 + 20 * degree[idn], edgecolors='k', alpha=0.7)
+        annotate3D(self.ax_3d, s=self.labels, xyz=self.coords, fontsize=10, xytext=(-3, 3),
+                   textcoords='offset points', ha='right', va='bottom')
+
+    def projectMatrix(self, M, coords):
+        proj_coords = []
+        for coord in coords:
+            coord = np.append(coord, 1.0)
+            proj_coord = M.dot(coord)
+            proj_coords.append(proj_coord)
+        return np.asarray(proj_coords)
+
+    def onRelease(self, event):
+        if event.inaxes == self.ax_3d:
+            M = event.inaxes.get_proj()
+            self.proj_coords = self.projectMatrix(M, self.coords)
+            self.plotType(self.prevlabel)
+
+    def plotScatter(self, x, y):
+        self.ax_2d.clear()
+        self.ax_2d.scatter(x, y, color="green")
+        self.ax_2d.set_title("Projection Minimum Spanning Tree")
+        self.fig.canvas.draw()
+
+    def plotContour(self, x, y):
+        self.ax_2d.clear()
+        rangeX = np.max(x) - np.min(x)
+        rangeY = np.max(y) - np.min(y)
+        if rangeX > rangeY:
+            sigma = rangeX / 50
+        else:
+            sigma = rangeY / 50
+        xi = np.linspace(min(x) - 0.1, max(x) + 0.1, 100)
+        yi = np.linspace(min(x) - 0.1, max(x) + 0.1, 100)
+        z = np.zeros((100, 100), float)
+        zSize = z.shape
+        N = len(x)
+        for c in range(zSize[1]):
+            for r in range(zSize[0]):
+                for d in range(N):
+                    z[r, c] = z[r, c] + (1.0 / N) * (1.0 / ((2 * math.pi) * sigma ** 2)) * math.exp(
+                        -((xi[c] - x[d]) ** 2 + (yi[r] - y[d]) ** 2) / (2 * sigma ** 2))
+        zMax = np.max(z)
+        z = z / zMax
+        self.ax_2d.contour(xi, yi, z, 15, linewidths=0.5, colors='k')
+        cf = self.ax_2d.contourf(xi, yi, z, 15, cmap=plt.cm.jet)
+        self.ax_2d.set_title("Projection Scatter Plot")
+        cbaxes = self.fig.add_axes([0.92, 0.1, 0.01, 0.8])
+        self.cb = self.fig.colorbar(mappable=cf, cax=cbaxes)
+        self.cb.set_ticks([])
+        self.fig.canvas.draw()
+
+    def plotConvolution(self, x, y):
+        coordinates = np.stack((x, y), axis=1)
+        Xr = np.round(coordinates, decimals=3)
+        size_grid = 2 * np.amax(Xr)
+        grid_coords = np.arange(-size_grid, size_grid, 0.001)
+        R, C = np.meshgrid(grid_coords, grid_coords, indexing='ij')
+        S = np.zeros(R.shape)
+        sigma = R.shape[0] / (200 / 5)
+        lbox = int(6 * sigma)
+        if lbox % 2 == 0:
+            lbox += 1
+        mid = int((lbox - 1) / 2 + 1)
+        kernel = np.zeros((lbox, lbox))
+        kernel[mid, mid] = 1
+        kernel = gaussian_filter(kernel, sigma=sigma)
+        for p in range(Xr.shape[0]):
+            indx = np.argmin(np.abs(R[:, 0] - Xr[p, 0]))
+            indy = np.argmin(np.abs(C[0, :] - Xr[p, 1]))
+            if 'weights' in locals():
+                S[indx - mid:indx + mid - 1, indy - mid:indy + mid - 1] += kernel * self.weights
+            else:
+                S[indx - mid:indx + mid - 1, indy - mid:indy + mid - 1] += kernel
+        S = S[~np.all(S == 0, axis=1)]
+        S = S[:, ~np.all(S == 0, axis=0)]
+        S = ndimage.rotate(S, 90)
+        cf = self.ax_2d.imshow(S)
+        cbaxes = self.fig.add_axes([0.92, 0.1, 0.01, 0.8])
+        self.ax_2d.set_title('Projection Scatter Plot')
+        self.cb = self.fig.colorbar(mappable=cf, cax=cbaxes)
+        self.cb.set_ticks([])
+        self.fig.canvas.draw()
+
+    def plotType(self, label):
+        self.prevlabel = label
+        x = self.proj_coords[:, 0]
+        y = self.proj_coords[:, 1]
+        if self.cb != None:
+            self.cb.remove()
+            self.cb = None
+        if label == 'Scatter':
+            self.plotScatter(x, y)
+        elif label == 'Contour':
+            self.plotContour(x, y)
+        elif label == 'Convolution':
+            self.plotConvolution(x, y)
+
+    def initializePlot(self):
+        self.mst_3D()
+        M = self.ax_3d.get_proj()
+        self.proj_coords = self.projectMatrix(M, self.coords)
+        self.ax_2d.scatter(self.proj_coords[:, 0], self.proj_coords[:, 1], color="green")
+
+        # Buttons
+        axcolor = 'silver'
+        rax = self.fig.add_axes([0.01, 0.4, 0.12, 0.25], facecolor=axcolor)
+        self.radio = RadioButtons(rax, ('Scatter', 'Contour', 'Convolution'))
+        self.radio.on_clicked(self.plotType)
+        self.canvas.get_tk_widget().pack(side=tk.BOTTOM, fill=tk.BOTH, expand=1)
+        self.canvas._tkcanvas.pack(side=tk.BOTTOM, fill=tk.BOTH, expand=1)
+
+        # Toolbar
+        toolbar = NavigationToolbar2Tk(self.canvas, self.root)
+        toolbar.update()
+        self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
+
+        tk.mainloop()
+
+class Annotation3D(Annotation):
+    '''Annotate the point xyz with text s'''
+
+    def __init__(self, s, xyz, *args, **kwargs):
+        Annotation.__init__(self, "", xy=(0, 0), bbox=dict(boxstyle="round,pad=0.3", fc="whitesmoke", ec="lightgray", lw=2),
+                            *args, **kwargs)
+        self._verts3d = xyz
+        self.s = s
+
+    def draw(self, renderer):
+        for coord, text in zip(self._verts3d, self.s):
+            xs3d, ys3d, zs3d = coord[0], coord[1], coord[2]
+            xs, ys, zs = proj_transform(xs3d, ys3d, zs3d, renderer.M)
+            self.xy = (xs, ys)
+            Annotation.set_text(self, text)
+            Annotation.draw(self, renderer)
+
+def annotate3D(ax, s, *args, **kwargs):
+    '''add anotation text s to to Axes3d ax'''
+
+    tag = Annotation3D(s, *args, **kwargs)
+    ax.add_artist(tag)
diff --git a/xmipp3/viewers/viewer_structure_map_sph.py b/xmipp3/viewers/viewer_structure_map_sph.py
deleted file mode 100644
index 26df81fab..000000000
--- a/xmipp3/viewers/viewer_structure_map_sph.py
+++ /dev/null
@@ -1,160 +0,0 @@
-########
-
-# This viewer has been removed as the protocols associated to it are no longer offered in Scipion.
-# This will be solved in a future release.
-
-########
-
-
-
-# **************************************************************************
-# *
-# * Authors:     Amaya Jimenez Moreno (ajimenez@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 os
-# import numpy as np
-# from mpl_toolkits.mplot3d import proj3d
-#
-# from pyworkflow.viewer import DESKTOP_TKINTER, WEB_DJANGO, ProtocolViewer
-# from pyworkflow.gui.plotter import plt
-# import pyworkflow.protocol.params as params
-#
-# from xmipp3.protocols.protocol_structure_map_sph import XmippProtStructureMapSPH
-#
-#
-# class XmippProtStructureMapSphViewer(ProtocolViewer):
-#     """ Wrapper to visualize different type of data objects
-#     with the Xmipp program xmipp_showj
-#     """
-#
-#     _label = 'viewer structure map sph'
-#     _environments = [DESKTOP_TKINTER, WEB_DJANGO]
-#     _targets = [XmippProtStructureMapSPH]
-#
-#
-#     def _defineParams(self, form):
-#         form.addSection(label='Show StructMap')
-#         form.addParam('numberOfDimensions', params.IntParam, default=2,
-#                       label="Number of dimensions",
-#                       help='In normal practice, it should be 1, 2 or, at most, 3.')
-#         form.addParam('doShowPlot', params.LabelParam,
-#                       label="Display the StructMap")
-#
-#     def _getVisualizeDict(self):
-#         return {'doShowPlot': self._visualize}
-#
-#
-#     def _visualize(self, e=None):
-#         nDim = self.numberOfDimensions.get()
-#         fnOutput = self.protocol._defineResultsName(nDim)
-#         if not os.path.exists(fnOutput):
-#             return [self.errorMessage('The necessary metadata was not produced\n'
-#                                       'Execute again the protocol\n',
-#                                       title='Missing result file')]
-#         coordinates = np.loadtxt(fnOutput)
-#
-#         # Create labels
-#         count = 0
-#         labels = []
-#         volList, _, _ = self.protocol._iterInputVolumes()
-#         for voli in volList:
-#             count+=1
-#             labels.append("vol_%02d"%count)
-#
-#         val = 0
-#         if nDim == 1:
-#             AX = plt.subplot(111)
-#             plot = plt.plot(coordinates, np.zeros_like(coordinates) + val, 'o', c='g')
-#             plt.xlabel('Dimension 1', fontsize=11)
-#             AX.set_yticks([1])
-#             plt.title('StructMap')
-#
-#             for label, x, y in zip(labels, coordinates, np.zeros_like(coordinates)):
-#                 plt.annotate(label,
-#                              xy = (x, y), xytext = (x+val, y+val),
-#                              textcoords = 'data', ha = 'right', va = 'bottom',fontsize=9,
-#                              bbox = dict(boxstyle = 'round,pad=0.3', fc = 'yellow', alpha = 0.3))
-#             plt.grid(True)
-#             plt.show()
-#
-#         elif nDim == 2:
-#             plot = plt.scatter(coordinates[:, 0], coordinates[:, 1], marker='o', c='g')
-#             plt.xlabel('Dimension 1', fontsize=11)
-#             plt.ylabel('Dimension 2', fontsize=11)
-#             plt.title('StructMap')
-#             for label, x, y in zip(labels, coordinates[:, 0], coordinates[:, 1]):
-#                 plt.annotate(label,
-#                              xy = (x, y), xytext = (x+val, y+val),
-#                              textcoords = 'data', ha = 'right', va = 'bottom',fontsize=9,
-#                              bbox = dict(boxstyle = 'round,pad=0.3', fc = 'yellow', alpha = 0.3))
-#             plt.grid(True)
-#             plt.show()
-#
-#         else:
-#
-#             fig = plt.figure()
-#             ax = fig.add_subplot(111, projection = '3d')
-#
-#             ax.scatter(coordinates[:, 0], coordinates[:, 1], coordinates[:, 2], marker = 'o', c='g')
-#             ax.set_xlabel('Dimension 1', fontsize=11)
-#             ax.set_ylabel('Dimension 2', fontsize=11)
-#             ax.set_zlabel('Dimension 3', fontsize=11)
-#             ax.text2D(0.05, 0.95, "StructMap", transform=ax.transAxes)
-#
-#             x2, y2, _ = proj3d.proj_transform(coordinates[:, 0], coordinates[:, 1], coordinates[:, 2], ax.get_proj())
-#             Labels = []
-#             for i in range(len(coordinates[:, 0])):
-#                 text = labels[i]
-#                 label = ax.annotate(text,
-#                                     xycoords='data',
-#                                     xy = (x2[i], y2[i]), xytext = (x2[i]+val, y2[i]+val),
-#                                     textcoords = 'data', ha = 'right',
-#                                      va = 'bottom', fontsize=9,
-#                                     bbox = dict(boxstyle = 'round,pad=0.3',
-#                                                  fc = 'yellow', alpha = 0.3))
-#
-#                 Labels.append(label)
-#
-#             def update_position(e):
-#                 x2, y2, _ = proj3d.proj_transform(coordinates[:, 0], coordinates[:, 1], coordinates[:, 2], ax.get_proj())
-#                 for i in range(len(coordinates[:, 0])):
-#                     label = Labels[i]
-#                     label.xytext = (x2[i],y2[i])
-#                     label.update_positions(fig.canvas.get_renderer())
-#                 fig.canvas.draw()
-#             fig.canvas.mpl_connect('button_release_event', update_position)
-#             plt.show()
-#
-#         return plot
-#
-#     def _validate(self):
-#         errors = []
-#
-#         numberOfDimensions=self.numberOfDimensions.get()
-#         if numberOfDimensions > 3 or numberOfDimensions < 1:
-#             errors.append("The number of dimensions should be 1, 2 or, at most, 3.")
-#
-#         return errors
-#
\ No newline at end of file
diff --git a/xmipp3/viewers/viewer_volume_deform_sph.py b/xmipp3/viewers/viewer_volume_deform_sph.py
deleted file mode 100644
index 4a12481dd..000000000
--- a/xmipp3/viewers/viewer_volume_deform_sph.py
+++ /dev/null
@@ -1,180 +0,0 @@
-########
-
-# This viewer has been removed as the protocols associated to it are no longer offered in Scipion.
-# This will be solved in a future release.
-
-########
-
-
-
-# **************************************************************************
-# *
-# * Authors:  Amaya Jimenez Moreno (ajimenez@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.viewer import DESKTOP_TKINTER, WEB_DJANGO, ProtocolViewer
-# from pwem.viewers import ChimeraView
-# import pyworkflow.protocol.params as params
-# from pyworkflow.utils import getExt, removeExt
-# from os.path import abspath
-# from xmipp3.protocols.protocol_volume_deform_sph import XmippProtVolumeDeformSPH
-#
-#
-# class XmippVolumeDeformSphViewer(ProtocolViewer):
-#     """ Visualize the output of protocol volume strain """
-#     _label = 'viewer volume deform sph'
-#     _targets = [XmippProtVolumeDeformSPH]
-#     _environments = [DESKTOP_TKINTER, WEB_DJANGO]
-#
-#
-#     def _defineParams(self, form):
-#         form.addSection(label='Show deformation')
-#         form.addParam('doShowStrain', params.LabelParam,
-#                       label="Display the strain deformation")
-#         form.addParam('doShowRotation', params.LabelParam,
-#                       label="Display the rotation deformation")
-#         form.addParam('doShowMorphOrigRef', params.LabelParam,
-#                       label="Display the morphing between original and reference volumes")
-#         form.addParam('doShowMorphDeformedRef', params.LabelParam,
-#                       label="Display the morphing between deformed and reference volumes")
-#
-#     def _getVisualizeDict(self):
-#         self.protocol._createFilenameTemplates()
-#         return {'doShowStrain': self._doShowStrain,
-#                 'doShowRotation': self._doShowRotation,
-#                 'doShowMorphOrigRef': self._doShowMorphOrigRef,
-#                 'doShowMorphDeformedRef': self._doShowDeformedOrigRef
-#                 }
-#
-#
-#     def _doShowStrain(self, param=None):
-#
-#         scriptFile = self.protocol._getPath('strain_chimera.cxc')
-#         fhCmd = open(scriptFile, 'w')
-#         fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
-#         ext = getExt(self.protocol._getFileName('fnInputVol'))
-#         fninput = abspath(fnbase + ext[0:4])
-#         smprt = self.protocol.outputVolume.getSamplingRate()
-#
-#         fnbase2 = removeExt(self.protocol._getFileName('fnOutVol'))
-#         fnStrain = abspath(fnbase2)
-#
-#         fhCmd.write("open %s\n" % fninput)
-#         fhCmd.write("open %s\n" % (fnStrain+"_strain.mrc"))
-#         counter = 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         counter += 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         fhCmd.write("vol #%d hide\n" % counter)
-#         #fhCmd.write("focus\n")
-#         fhCmd.write('color sample #%d map #%d palette ^rainbow\n' % (counter -1, counter))
-#         fhCmd.close()
-#
-#         view = ChimeraView(scriptFile)
-#         return [view]
-#
-#
-#     def _doShowRotation(self, param=None):
-#
-#         scriptFile = self.protocol._getPath('rotation_chimera.cxc')
-#         fhCmd = open(scriptFile, 'w')
-#         fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
-#         ext = getExt(self.protocol._getFileName('fnInputVol'))
-#         fninput = abspath(fnbase + ext[0:4])
-#         smprt = self.protocol.outputVolume.getSamplingRate()
-#
-#         fnbase2 = removeExt(self.protocol._getFileName('fnOutVol'))
-#         fnStrain = abspath(fnbase2)
-#
-#         fhCmd.write("open %s\n" % fninput)
-#         fhCmd.write("open %s\n" % (fnStrain+"_rotation.mrc"))
-#         counter = 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         counter += 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         #fhCmd.write("focus\n")
-#         fhCmd.write("vol #%d hide\n" % (counter))
-#         fhCmd.write('color sample #%d map #%d palette ^rainbow\n' % (counter -1, counter))
-#         # fhCmd.write("scolor #0 volume #1 cmap rainbow reverseColors True\n")
-#         fhCmd.close()
-#
-#         view = ChimeraView(scriptFile)
-#         return [view]
-#
-#     def _doShowMorphOrigRef(self, param=None):
-#
-#         scriptFile = self.protocol._getPath('morph_orig_ref_chimera.cxc')
-#         fhCmd = open(scriptFile, 'w')
-#         fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
-#         ext = getExt(self.protocol._getFileName('fnInputVol'))
-#         fninput = abspath(fnbase + ext[0:4])
-#         smprt = self.protocol.outputVolume.getSamplingRate()
-#
-#         fnbase2 = removeExt(self.protocol._getFileName('fnRefVol'))
-#         ext2 = getExt(self.protocol._getFileName('fnRefVol'))
-#         fnref = abspath(fnbase2 + ext2[0:4])
-#
-#         fhCmd.write("open %s\n" % fninput)
-#         fhCmd.write("open %s\n" % fnref)
-#
-#         counter = 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         fhCmd.write("vol #%d hide\n" % (counter))
-#         counter += 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         fhCmd.write("vol #%d hide\n" % (counter))
-#         #fhCmd.write("focus\n")
-#         fhCmd.write("volume morph #%d,%d frames 500\n" % (counter - 1, counter))
-#         fhCmd.close()
-#
-#         view = ChimeraView(scriptFile)
-#         return [view]
-#
-#     def _doShowDeformedOrigRef(self, param=None):
-#
-#         scriptFile = self.protocol._getPath('morph_deformed_ref_chimera.cxc')
-#         fhCmd = open(scriptFile, 'w')
-#         fnbase = removeExt(self.protocol._getFileName('fnOutVol'))
-#         ext = getExt(self.protocol._getFileName('fnOutVol'))
-#         fninput = abspath(fnbase + ext[0:4])
-#         smprt = self.protocol.outputVolume.getSamplingRate()
-#
-#         fnbase2 = removeExt(self.protocol._getFileName('fnRefVol'))
-#         ext2 = getExt(self.protocol._getFileName('fnRefVol'))
-#         fnref = abspath(fnbase2 + ext2[0:4])
-#
-#         fhCmd.write("open %s\n" % fninput)
-#         fhCmd.write("open %s\n" % fnref)
-#         counter = 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         fhCmd.write("vol #%d hide\n" % counter)
-#         counter += 1
-#         fhCmd.write("volume #%d voxelSize %s\n" % (counter, str(smprt)))
-#         # fhCmd.write("focus\n")
-#         fhCmd.write("vol #%d hide\n" % counter)
-#         fhCmd.write("volume morph #%d,%d frames 500\n" %  (counter-1, counter))
-#         fhCmd.close()
-#
-#         view = ChimeraView(scriptFile)
-#         return [view]
diff --git a/xmipp3/viewers/viewer_volume_deform_zernike3d.py b/xmipp3/viewers/viewer_volume_deform_zernike3d.py
new file mode 100644
index 000000000..e1da88de2
--- /dev/null
+++ b/xmipp3/viewers/viewer_volume_deform_zernike3d.py
@@ -0,0 +1,176 @@
+# **************************************************************************
+# *
+# * Authors:  Amaya Jimenez Moreno  (ajimenez@cnb.csic.es)
+# *           David Herreros Calero (dherreros@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.viewer import DESKTOP_TKINTER, WEB_DJANGO, ProtocolViewer
+from pwem.viewers import ChimeraView
+import pyworkflow.protocol.params as params
+from pyworkflow.utils import getExt, removeExt
+from os.path import abspath
+from xmipp3.protocols.protocol_volume_deform_zernike3d import XmippProtVolumeDeformZernike3D
+
+
+class XmippVolumeDeformZernike3DViewer(ProtocolViewer):
+    """ Visualize the output of protocol volume strain """
+    _label = 'viewer volume deform sph'
+    _targets = [XmippProtVolumeDeformZernike3D]
+    _environments = [DESKTOP_TKINTER, WEB_DJANGO]
+    OPEN_FILE = "open %s\n"
+    VOXEL_SIZE = "volume #%d voxelSize %s\n"
+    VOL_HIDE = "vol #%d hide\n"
+    VIEW = "view\n"
+
+    def _defineParams(self, form):
+        form.addSection(label='Show deformation')
+        form.addParam('doShowStrain', params.LabelParam,
+                      label="Display the strain deformation")
+        form.addParam('doShowRotation', params.LabelParam,
+                      label="Display the rotation deformation")
+        form.addParam('doShowMorphOrigRef', params.LabelParam,
+                      label="Display the morphing between original and reference volumes")
+        form.addParam('doShowMorphDeformedRef', params.LabelParam,
+                      label="Display the morphing between deformed and reference volumes")
+
+    def _getVisualizeDict(self):
+        self.protocol._createFilenameTemplates()
+        return {'doShowStrain': self._doShowStrain,
+                'doShowRotation': self._doShowRotation,
+                'doShowMorphOrigRef': self._doShowMorphOrigRef,
+                'doShowMorphDeformedRef': self._doShowDeformedOrigRef
+                }
+
+
+    def _doShowStrain(self, param=None):
+
+        scriptFile = self.protocol._getPath('strain_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
+        ext = getExt(self.protocol._getFileName('fnInputVol'))
+        fninput = abspath(fnbase + ext[0:4])
+        smprt = self.protocol.outputVolume.getSamplingRate()
+
+        fnbase2 = removeExt(self.protocol._getFileName('fnOutVol'))
+        fnStrain = abspath(fnbase2)
+
+        fhCmd.write(self.OPEN_FILE % fninput)
+        fhCmd.write(self.OPEN_FILE % (fnStrain+"_strain.mrc"))
+        counter = 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        counter += 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        fhCmd.write(self.VOL_HIDE % counter)
+        fhCmd.write('color sample #%d map #%d palette rainbow\n' % (counter - 1, counter))
+        fhCmd.write(self.VIEW)
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
+
+
+    def _doShowRotation(self, param=None):
+
+        scriptFile = self.protocol._getPath('rotation_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
+        ext = getExt(self.protocol._getFileName('fnInputVol'))
+        fninput = abspath(fnbase + ext[0:4])
+        smprt = self.protocol.outputVolume.getSamplingRate()
+
+        fnbase2 = removeExt(self.protocol._getFileName('fnOutVol'))
+        fnStrain = abspath(fnbase2)
+
+        fhCmd.write(self.OPEN_FILE % fninput)
+        fhCmd.write(self.OPEN_FILE % (fnStrain+"_rotation.mrc"))
+        counter = 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        counter += 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        fhCmd.write(self.VOL_HIDE % (counter))
+        fhCmd.write('color sample #%d map #%d palette rainbow\n' % (counter - 1, counter))
+        fhCmd.write(self.VIEW)
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
+
+    def _doShowMorphOrigRef(self, param=None):
+
+        scriptFile = self.protocol._getPath('morph_orig_ref_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        fnbase = removeExt(self.protocol._getFileName('fnInputVol'))
+        ext = getExt(self.protocol._getFileName('fnInputVol'))
+        fninput = abspath(fnbase + ext[0:4])
+        smprt = self.protocol.outputVolume.getSamplingRate()
+
+        fnbase2 = removeExt(self.protocol._getFileName('fnRefVol'))
+        ext2 = getExt(self.protocol._getFileName('fnRefVol'))
+        fnref = abspath(fnbase2 + ext2[0:4])
+
+        fhCmd.write(self.OPEN_FILE % fninput)
+        fhCmd.write(self.OPEN_FILE % fnref)
+
+        counter = 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        fhCmd.write(self.VOL_HIDE % (counter))
+        counter += 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        fhCmd.write(self.VOL_HIDE % (counter))
+        fhCmd.write("volume morph #%d,%d frames 500\n" % (counter - 1, counter))
+        fhCmd.write(self.VIEW)
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
+
+    def _doShowDeformedOrigRef(self, param=None):
+
+        scriptFile = self.protocol._getPath('morph_deformed_ref_chimera.cxc')
+        fhCmd = open(scriptFile, 'w')
+        fnbase = removeExt(self.protocol._getFileName('fnOutVol'))
+        ext = getExt(self.protocol._getFileName('fnOutVol'))
+        fninput = abspath(fnbase + ext[0:4])
+        smprt = self.protocol.outputVolume.getSamplingRate()
+
+        fnbase2 = removeExt(self.protocol._getFileName('fnRefVol'))
+        ext2 = getExt(self.protocol._getFileName('fnRefVol'))
+        fnref = abspath(fnbase2 + ext2[0:4])
+
+        fhCmd.write(self.OPEN_FILE % fninput)
+        fhCmd.write(self.OPEN_FILE % fnref)
+        counter = 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        fhCmd.write(self.VOL_HIDE % counter)
+        counter += 1
+        fhCmd.write(self.VOXEL_SIZE % (counter, str(smprt)))
+        # fhCmd.write("focus\n")
+        fhCmd.write(self.VOL_HIDE % counter)
+        fhCmd.write("volume morph #%d,%d frames 500\n" %  (counter-1, counter))
+        fhCmd.write(self.VIEW)
+        fhCmd.close()
+
+        view = ChimeraView(scriptFile)
+        return [view]
diff --git a/xmipp3/wizards.py b/xmipp3/wizards.py
index 14a8e3fe7..10217c725 100644
--- a/xmipp3/wizards.py
+++ b/xmipp3/wizards.py
@@ -43,7 +43,7 @@
     XmippProtHelicalParameters, XmippProtConsensusPicking, XmippProtMonoRes,
     XmippProtRotSpectra, XmippProtReconstructHighRes, XmippProtExtractUnit,
     XmippProtReconstructHeterogeneous, XmippMetaProtDiscreteHeterogeneityScheduler,
-    XmippProtShiftParticles)
+    XmippProtShiftParticles, XmippProtVolumeDeformZernike3D, XmippProtStructureMapZernike3D)
 
 
 #===============================================================================
@@ -330,8 +330,7 @@ def _getParameters(self, protocol):
 
 class XmippVolumeRadiiWizard(VolumeMaskRadiiWizard):
     _targets = [(XmippProtMaskVolumes, ['innerRadius', 'outerRadius']),
-               (XmippProtExtractUnit, ['innerRadius', 'outerRadius'])
-              ]
+               (XmippProtExtractUnit, ['innerRadius', 'outerRadius'])]
 
     def _getParameters(self, protocol):
 
@@ -368,6 +367,33 @@ def _getParameters(self, protocol):
         protParams['value']= value
         return protParams
 
+class Zernike3DMaskWizard(VolumeMaskRadiusWizard):
+    _targets = [(XmippProtVolumeDeformZernike3D, ['Rmax']),
+                (XmippProtStructureMapZernike3D, ['Rmax'])]
+
+    def _getParameters(self, protocol):
+
+        label, value = self._getInputProtocol(self._targets, protocol)
+
+        protParams = {}
+        if isinstance(protocol, XmippProtVolumeDeformZernike3D):
+            protParams['input'] = protocol.inputVolume
+        else:
+            protParams['input'] = protocol.inputVolumes
+        protParams['label'] = label
+        protParams['value'] = value
+        return protParams
+
+    def _getProvider(self, protocol):
+        _objs = self._getParameters(protocol)['input']
+        return VolumeMaskRadiusWizard._getListProvider(self, _objs)
+
+    def show(self, form):
+        params = self._getParameters(form.protocol)
+        _value = params['value']
+        _label = params['label']
+        VolumeMaskRadiusWizard.show(self, form, _value, _label, UNIT_PIXEL)
+
 
 #===============================================================================
 #  FILTERS


Test patch:
diff --git a/xmipp3/tests/test_protocol_screen_deepConsensus.py b/xmipp3/tests/test_protocol_screen_deepConsensus.py
index 74cb3b38b..6cb56397d 100644
--- a/xmipp3/tests/test_protocol_screen_deepConsensus.py
+++ b/xmipp3/tests/test_protocol_screen_deepConsensus.py
@@ -33,10 +33,12 @@
                             ProtCreateStreamData)
 from pwem.protocols.protocol_create_stream_data import SET_OF_MICROGRAPHS, SET_OF_COORDINATES
 
-
 from pwem import emlib
 from xmipp3.protocols import XmippProtScreenDeepConsensus
 
+ADD_MODEL_TRAIN_NEW = 0
+ADD_MODEL_TRAIN_PRETRAIN = 1
+ADD_MODEL_TRAIN_PREVRUN = 2
 
 class TestXmippProtScreenDeepConsensus(BaseTest):
     @classmethod
@@ -96,15 +98,7 @@ def _runDeepConsensusPicking(self, inpCoords, kwargs, case = ''):
                            "There was a problem with the consensus")
       return prot
 
-    def getDeepConsensusKwargs(self, case=1):
-      ADD_MODEL_TRAIN_PRETRAIN = 1
-      ADD_MODEL_TRAIN_PREVRUN = 2
-
-      ADD_DATA_TRAIN_CUST = 2
-
-      ADD_DATA_TRAIN_CUSTOM_OPT_PARTS = 0
-      ADD_DATA_TRAIN_CUSTOM_OPT_COORS = 1
-
+    def getDeepConsensusKwargs(self, case=1, inputCase=1):
       kwargs = {
         'nEpochs' : 1.0,
         'nModels' :2,
@@ -112,58 +106,62 @@ def getDeepConsensusKwargs(self, case=1):
         'trainingBatch':3,
         'predictingBatch':3
       }
+      #inputCase controls the input model of the protocol: previous protocol model, new model, pretrained
+      inputCaseKwargs = {'numberOfThreads': 1} if inputCase<4 else {'numberOfThreads': 4}
+      if inputCase in [1, 4]:
+        inputCaseKwargs['modelInitialization'] = ADD_MODEL_TRAIN_NEW
+      elif inputCase in [2, 5]:
+        inputCaseKwargs['modelInitialization'] = ADD_MODEL_TRAIN_PRETRAIN
+      elif inputCase in [3, 6]:
+        inputCaseKwargs['modelInitialization'] = ADD_MODEL_TRAIN_PREVRUN
+        inputCaseKwargs['continueRun'] = self.lastRun
+
+      #case controls the behavior of the protocol. e.g: skip training, do preliminar predictions
       if case == 1:
-        #Simplest case
-        caseKwargs = {}
+        caseKwargs = {'doPreliminarPredictions': False, 'skipTraining': False}
       elif case == 2:
-        #Pretrained model, do testing, additional data particles, do preliminar predictions
-        caseKwargs = {'doPreliminarPredictions': True,
-                      'modelInitialization': ADD_MODEL_TRAIN_PRETRAIN,
-                      'addTrainingData': ADD_DATA_TRAIN_CUST,
-                      'trainingDataType':ADD_DATA_TRAIN_CUSTOM_OPT_PARTS,
-                      'trainTrueSetOfParticles': self.lastRun.outputParticles,
-                      'trainFalseSetOfParticles': self.lastRun.outputParticles
-                      }
+        caseKwargs = {'doPreliminarPredictions': False, 'skipTraining': True}
       elif case == 3:
-        #Previous run model, additional data coords
-        caseKwargs = {'modelInitialization':ADD_MODEL_TRAIN_PREVRUN,
-                      'continueRun':self.lastRun,
-                      'doTesting': True,
-                      'testTrueSetOfParticles': self.lastRun.outputParticles,
-                      'testFalseSetOfParticles': self.lastRun.outputParticles,
-                      'addTrainingData': ADD_DATA_TRAIN_CUST,
-                      'trainingDataType': ADD_DATA_TRAIN_CUSTOM_OPT_COORS,
-                      'trainTrueSetOfCoords': self.lastRun.outputCoordinates,
-                      'trainFalseSetOfCoords': self.lastRun.outputCoordinates
-                      }
+        caseKwargs = {'doPreliminarPredictions': True, 'skipTraining': False}
+      elif case == 4:
+        caseKwargs = {'doPreliminarPredictions': True, 'skipTraining': True}
+
       else:
         caseKwargs={}
 
       kwargs.update(caseKwargs)
+      kwargs.update(inputCaseKwargs)
       return kwargs
 
-
-    def testImportCoordinates(self):
+    def testDeepConsensusNew(self):
+      #Testing the protocol on new model
       nCoordinateSets = 3
-      protImpMics = self._runInportMicrographs()
+      inputCase, case = 4, 3
 
-      for case in range(1,3):
-        protImpCoords = self._runImportCoordinates(protImpMics, case)
+      inpCoords = self.prepareInput(case, nCoordinateSets)
+      kwargs = self.getDeepConsensusKwargs(case, inputCase)
+      self.lastRun = self._runDeepConsensusPicking(inpCoords, kwargs, case)
 
-        protsStreamCoords = []
-        for i in range(nCoordinateSets):
-          protsStreamCoords.append(self._runStreamCoordinates(protImpCoords, i))
+    def testDeepConsensusPretrain(self):
+      #Testing the protocol on pretrained model
+      nCoordinateSets = 3
+      inputCase, case = 5, 1
 
-        for i in range(nCoordinateSets):
-          self._waitOutput(protsStreamCoords[i], "outputCoordinates")
+      inpCoords = self.prepareInput(case, nCoordinateSets)
+      kwargs = self.getDeepConsensusKwargs(case, inputCase)
+      self.lastRun = self._runDeepConsensusPicking(inpCoords, kwargs, case)
 
-        inpCoords = []
-        for prot in protsStreamCoords:
-          point = Pointer(prot, extended="outputCoordinates")
-          inpCoords.append(point)
-        kwargs = self.getDeepConsensusKwargs(case)
+    def testDeepConsensusPrevRun(self):
+      #Testing the protocol on a model from a previous run
+      nCoordinateSets = 3
+      inputCase, case = 6, 4
+
+      inpCoords = self.prepareInput(case, nCoordinateSets)
+      prev_kwargs = self.getDeepConsensusKwargs(case=2, inputCase=5)
+      self.lastRun = self._runDeepConsensusPicking(inpCoords, prev_kwargs, case)
 
-        self.lastRun = self._runDeepConsensusPicking(inpCoords, kwargs, case)
+      kwargs = self.getDeepConsensusKwargs(case, inputCase)
+      self._runDeepConsensusPicking(inpCoords, kwargs, case)
 
 
     def _updateProtocol(self, prot):
@@ -181,4 +179,21 @@ def checkStreamOutput(self, prots, outputName):
           return False
       return True
 
+    def prepareInput(self, case, nCoordSets=3):
+      protImpMics = self._runInportMicrographs()
+      protImpCoords = self._runImportCoordinates(protImpMics, case)
+
+      protsStreamCoords = []
+      for i in range(nCoordSets):
+        protsStreamCoords.append(self._runStreamCoordinates(protImpCoords, i))
+
+      for i in range(nCoordSets):
+        self._waitOutput(protsStreamCoords[i], "outputCoordinates")
+
+      inpCoords = []
+      for prot in protsStreamCoords:
+        point = Pointer(prot, extended="outputCoordinates")
+        inpCoords.append(point)
+      return inpCoords
+
 
diff --git a/xmipp3/tests/test_protocols_sph.py b/xmipp3/tests/test_protocols_sph.py
deleted file mode 100644
index b10928baa..000000000
--- a/xmipp3/tests/test_protocols_sph.py
+++ /dev/null
@@ -1,145 +0,0 @@
-########
-
-# This test has been removed as the protocols associated to it are no longer offered in Scipion.
-# This will be solved in a future release.
-
-########
-
-
-
-# **************************************************************************
-# *
-# * Authors:    Carlos Oscar Sorzano (coss@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.tests import BaseTest, DataSet, setupTestProject
-#
-# from pwem.protocols import ProtImportVolumes, ProtImportParticles, ProtSubSet
-#
-# from xmipp3.protocols import (XmippProtAngularAlignmentSPH,
-#                               XmippProtVolumeDeformSPH,
-#                               XmippProtStructureMapSPH)
-#
-# class TestSPHBase(BaseTest):
-#     @classmethod
-#     def setData(cls):
-#         cls.dataset = DataSet.getDataSet('relion_tutorial')
-#         cls.particles = cls.dataset.getFile('import/case2/particles.sqlite')
-#         cls.volParticles = cls.dataset.getFile('volume')
-#         cls.volInput = cls.dataset.getFile('volumes/reference_rotated.vol')
-#         cls.volRef = cls.dataset.getFile('volumes/reference.mrc')
-#
-#     @classmethod
-#     def runImportVolume(cls, pattern, sampling):
-#         """ Run an Import volume protocol. """
-#         protImport = cls.newProtocol(ProtImportVolumes,
-#                                           filesPath=pattern,
-#                                           samplingRate=sampling)
-#         cls.launchProtocol(protImport)
-#         if protImport.outputVolume is None:
-#             raise Exception('Import of volume: %s, failed. outputVolume is None.' % pattern)
-#         return protImport
-#
-#     @classmethod
-#     def runImportParticles(cls, pattern, sampling):
-#         """ Run an Import particles protocol. """
-#         protImport = cls.newProtocol(ProtImportParticles,
-#                                  objLabel='Particles from scipion',
-#                                  importFrom=ProtImportParticles.IMPORT_FROM_SCIPION,
-#                                  sqliteFile=pattern,
-#                                  magnification=50000,
-#                                  samplingRate=sampling,
-#                                  haveDataBeenPhaseFlipped=True
-#                                  )
-#         cls.launchProtocol(protImport)
-#         if protImport.outputParticles is None:
-#             raise Exception('Import of particles: %s, failed. outputParticles is None.' % pattern)
-#         return protImport
-#
-#
-#
-# class TestAngularAlignmentSPH(TestSPHBase):
-#     @classmethod
-#     def setUpClass(cls):
-#         setupTestProject(cls)
-#         TestSPHBase.setData()
-#         cls.protImportPart = cls.runImportParticles(cls.particles, 7.08)
-#         cls.protImportVol = cls.runImportVolume(cls.volParticles, 7.08)
-#
-#     def test(self):
-#         subset = self.newProtocol(ProtSubSet,
-#                                   inputFullSet=self.protImportPart.outputParticles,
-#                                   chooseAtRandom=True,
-#                                   nElements=10)
-#         self.launchProtocol(subset)
-#
-#         alignSph = self.newProtocol(XmippProtAngularAlignmentSPH,
-#                                    inputParticles=subset.outputParticles,
-#                                    inputVolume=self.protImportVol.outputVolume,
-#                                    depth=1,
-#                                    targetResolution = 4.0,
-#                                    numberOfMpi=8)
-#         self.launchProtocol(alignSph)
-#         self.assertIsNotNone(alignSph.outputParticles,
-#                              "There was a problem with AngularAlignmentSPH")
-#
-#
-# class TestVolumeDeformSPH(TestSPHBase):
-#     @classmethod
-#     def setUpClass(cls):
-#         setupTestProject(cls)
-#         TestSPHBase.setData()
-#         cls.protImportVol = cls.runImportVolume(cls.volInput, 1.0)
-#         cls.protImportVolRef = cls.runImportVolume(cls.volRef, 1.0)
-#
-#     def test(self):
-#
-#         volDeformSph = self.newProtocol(XmippProtVolumeDeformSPH,
-#                                    inputVolume=self.protImportVol.outputVolume,
-#                                    refVolume = self.protImportVolRef.outputVolume,
-#                                    depth=2,
-#                                    targetResolution = 4.0)
-#         self.launchProtocol(volDeformSph)
-#         self.assertIsNotNone(volDeformSph.outputVolume,
-#                              "There was a problem with VolumeDeformSPH")
-#
-#
-#
-# class TestStructureMapSPH(TestSPHBase):
-#     @classmethod
-#     def setUpClass(cls):
-#         setupTestProject(cls)
-#         TestSPHBase.setData()
-#         cls.protImportVol1 = cls.runImportVolume(cls.volInput, 1.0)
-#         cls.protImportVol2 = cls.runImportVolume(cls.volRef, 1.0)
-#
-#     def test(self):
-#
-#         structureMap = self.newProtocol(XmippProtStructureMapSPH,
-#                                    depth=2,
-#                                    targetResolution = 4.0)
-#         volumeList = [self.protImportVol1.outputVolume, self.protImportVol2.outputVolume]
-#         structureMap.inputVolumes.set(volumeList)
-#         self.launchProtocol(structureMap)
-#         self.assertTrue(structureMap.isFinished(), "There was a problem with StructureMapSPH")
\ No newline at end of file
diff --git a/xmipp3/tests/test_protocols_xmipp_2d.py b/xmipp3/tests/test_protocols_xmipp_2d.py
index b76b95ce6..e9ef6756a 100644
--- a/xmipp3/tests/test_protocols_xmipp_2d.py
+++ b/xmipp3/tests/test_protocols_xmipp_2d.py
@@ -1006,7 +1006,7 @@ def setUpClass(cls):
         diameterInA = 2 * radiusInPixel * psize
         cls.protRelion2DClass = cls.newProtocol(ProtRelionClassify2D,
                                                 doCTF=False, maskDiameterA=diameterInA,
-                                                numberOfMpi=4, numberOfThreads=1)
+                                                numberOfMpi=1, numberOfThreads=1)
         cls.protRelion2DClass.numberOfClasses.set(4)
         cls.protRelion2DClass.numberOfIterations.set(3)
         cls.protRelion2DClass.inputParticles.set(cls.protNormalize.outputParticles)
diff --git a/xmipp3/tests/test_protocols_xmipp_3d.py b/xmipp3/tests/test_protocols_xmipp_3d.py
index 141b84168..b608cbcd3 100644
--- a/xmipp3/tests/test_protocols_xmipp_3d.py
+++ b/xmipp3/tests/test_protocols_xmipp_3d.py
@@ -48,6 +48,10 @@
                                                ALIGN_ALGORITHM_EXHAUSTIVE_LOCAL,
                                                ALIGN_ALGORITHM_LOCAL)
 
+MSG_WRONG_SAMPLING = "wrong sampling rate"
+MSG_WRONG_SIZE = "wrong size"
+MSG_MASK_ERROR = "There was a problem with create mask from volume"
+
 
 class TestXmippBase(BaseTest):
     """ Some utility functions to import volumes that are used in several tests."""
@@ -120,7 +124,7 @@ def testCreateMask(self):
         protMask1.setObjLabel('threshold mask')
         self.launchProtocol(protMask1)
         self.assertIsNotNone(protMask1.outputMask,
-                             "There was a problem with create mask from volume")
+                             MSG_MASK_ERROR)
 
 
         print("Run create segment mask from volume")
@@ -131,7 +135,7 @@ def testCreateMask(self):
         protMask2.setObjLabel('segmentation automatic')
         self.launchProtocol(protMask2)
         self.assertIsNotNone(protMask2.outputMask,
-                             "There was a problem with create mask from volume")
+                             MSG_MASK_ERROR)
 
     
         print("Run create mask from another mask")
@@ -789,6 +793,56 @@ def testExhaustiveLocal(self):
         protAlign.inputVolumes.append(self.protImport2.outputVolume)
         self.launchProtocol(protAlign)
 
+    def testExhaustiveCircularMask(self):
+        protAlign = self.newProtocol(XmippProtAlignVolume,
+                                     inputReference=self.protImport1.outputVolume,
+                                     alignmentAlgorithm=ALIGN_ALGORITHM_EXHAUSTIVE,
+                                     minRotationalAngle=65,
+                                     maxRotationalAngle=100,
+                                     stepRotationalAngle=10,
+                                     minTiltAngle=65,
+                                     maxTiltAngle=100,
+                                     stepTiltAngle=10,
+                                     minInplaneAngle=0,
+                                     maxInplaneAngle=0,
+                                     stepInplaneAngle=0,
+                                     applyMask=True,
+                                     maskType=0,
+                                     maskRadius=22,
+                                     numberOfMpi=1, numberOfThreads=1
+                                     )
+        protAlign.inputVolumes.append(self.protImport2.outputVolume)
+        self.launchProtocol(protAlign)
+
+    def testExhaustiveBinaryMask(self):
+        protMask = self.newProtocol(XmippProtCreateMask3D,
+                                    source=0,
+                                    volumeOperation=0,
+                                    threshold=0.0)
+        protMask.inputVolume.set(self.protImport1.outputVolume)
+        self.launchProtocol(protMask)
+        self.assertIsNotNone(protMask.outputMask,
+                             MSG_MASK_ERROR)
+        protAlign = self.newProtocol(XmippProtAlignVolume,
+                                     inputReference=self.protImport1.outputVolume,
+                                     alignmentAlgorithm=ALIGN_ALGORITHM_EXHAUSTIVE,
+                                     minRotationalAngle=65,
+                                     maxRotationalAngle=100,
+                                     stepRotationalAngle=10,
+                                     minTiltAngle=65,
+                                     maxTiltAngle=100,
+                                     stepTiltAngle=10,
+                                     minInplaneAngle=0,
+                                     maxInplaneAngle=0,
+                                     stepInplaneAngle=0,
+                                     applyMask=True,
+                                     maskType=1,
+                                     numberOfMpi=1, numberOfThreads=1
+                                     )
+        protAlign.inputVolumes.append(self.protImport2.outputVolume)
+        protAlign.maskFile.set(protMask.outputMask)
+        self.launchProtocol(protAlign)
+
 
 class TestXmippProtHelicalParameters(TestXmippBase):
     @classmethod
@@ -1032,7 +1086,8 @@ def testImportPdbFromFn(self):
         self.launchProtocol(protConvert)
         self.assertIsNotNone(protConvert.outputPdb.getFileName(), 
                              "There was a problem with the import")
-        
+
+
 class TestXmippPdbConvert(TestXmippBase):
     
     @classmethod
@@ -1048,8 +1103,9 @@ def testXmippPdbConvertFromDb(self):
         self.launchProtocol(protConvert)
         self.assertIsNotNone(protConvert.outputVolume.getFileName(), "There was a problem with the conversion")
         self.assertAlmostEqual(protConvert.outputVolume.getSamplingRate(), protConvert.sampling.get(), places=1,
-                               msg="wrong sampling rate")
-        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1, msg="wrong size")
+                               msg=MSG_WRONG_SAMPLING)
+        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1,
+                               msg=MSG_WRONG_SIZE)
         
     def testXmippPdbConvertFromObj(self):
         print("Run convert a pdb from import")
@@ -1066,8 +1122,9 @@ def testXmippPdbConvertFromObj(self):
         self.launchProtocol(protConvert)
         self.assertIsNotNone(protConvert.outputVolume.getFileName(), "There was a problem with the conversion")
         self.assertAlmostEqual(protConvert.outputVolume.getSamplingRate(), protConvert.sampling.get(), places=1,
-                               msg="wrong sampling rate")
-        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1, msg="wrong size")
+                               msg=MSG_WRONG_SAMPLING)
+        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1,
+                               msg=MSG_WRONG_SIZE)
 
     def testXmippPdbConvertFromFn(self):
         print("Run convert a pdb from file")
@@ -1075,8 +1132,8 @@ def testXmippPdbConvertFromFn(self):
         self.launchProtocol(protConvert)
         self.assertIsNotNone(protConvert.outputVolume.getFileName(), "There was a problem with the conversion")
         self.assertAlmostEqual(protConvert.outputVolume.getSamplingRate(), protConvert.sampling.get(), places=1,
-                               msg="wrong sampling rate")
-        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], 48, places=1, msg="wrong size")
+                               msg=MSG_WRONG_SAMPLING)
+        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], 48, places=1, msg=MSG_WRONG_SIZE)
 
 
 class TestXmippValidateNonTilt(TestXmippBase):
@@ -1172,20 +1229,110 @@ def testXmippVolSub(self):
         self.launchProtocol(protImportVol1)
         self.assertIsNotNone(protImportVol1.getFiles(),
                              "There was a problem with the import 1")
+
         print("Import Volume 2")
         protImportVol2 = self.newProtocol(ProtImportVolumes,
-                                         objLabel='Volume',
-                                         filesPath=self.vol2,
-                                         samplingRate=7.08)
+                                          objLabel='Volume',
+                                          filesPath=self.vol2,
+                                          samplingRate=7.08)
         self.launchProtocol(protImportVol2)
         self.assertIsNotNone(protImportVol2.getFiles(),
                              "There was a problem with the import 2")
+
+        print("Import atomic structure")
+        protImportPdb = self.newProtocol(ProtImportPdb,
+                                          pdbId='6Z9F',
+                                          inputVolume=protImportVol1.outputVolume)
+        self.launchProtocol(protImportPdb)
+        self.assertIsNotNone(protImportPdb.getFiles(),
+                             "There was a problem with the atomic structure")
+
+        print("Create mask")
+        protCreateMask = self.newProtocol(XmippProtCreateMask3D,
+                                          inputVolume=protImportVol1.outputVolume,
+                                          threshold=0.28)
+        self.launchProtocol(protCreateMask)
+        self.assertIsNotNone(protCreateMask.getFiles(),
+                             "There was a problem with the 3D mask")
+
+        print("Run volume adjust")
+        protVolAdj = self.newProtocol(XmippProtVolAdjust,
+                                      vol1=protImportVol1.outputVolume,
+                                      vol2=protImportVol2.outputVolume,
+                                      masks=False)
+        self.launchProtocol(protVolAdj)
+        self.assertIsNotNone(protVolAdj.outputVolume,
+                             "There was a problem with Volumes adjust")
+
+        protVolAdjNoE = self.newProtocol(XmippProtVolAdjust,
+                                         vol1=protImportVol1.outputVolume,
+                                         vol2=protImportVol2.outputVolume,
+                                         computeE=False,
+                                         masks=False)
+        self.launchProtocol(protVolAdjNoE)
+        self.assertIsNotNone(protVolAdjNoE.outputVolume,
+                             "There was a problem with Volumes adjust without computing energy")
+
+        protVolAdjNoRadAvg = self.newProtocol(XmippProtVolAdjust,
+                                              vol1=protImportVol1.outputVolume,
+                                              vol2=protImportVol2.outputVolume,
+                                              masks=False,
+                                              radavg=False)
+        self.launchProtocol(protVolAdjNoRadAvg)
+        self.assertIsNotNone(protVolAdjNoRadAvg.outputVolume,
+                             "There was a problem with Volumes adjust without radial average")
+
         print("Run volume subtraction")
         protVolSub = self.newProtocol(XmippProtVolSubtraction,
                                       vol1=protImportVol1.outputVolume,
                                       vol2=protImportVol2.outputVolume,
-                                      pdb=False,
-                                      masks=False)
+                                      masks=False,
+                                      radavg=False)
         self.launchProtocol(protVolSub)
         self.assertIsNotNone(protVolSub.outputVolume,
                              "There was a problem with Volumes subtraction")
+
+        protVolSubNoE = self.newProtocol(XmippProtVolSubtraction,
+                                         vol1=protImportVol1.outputVolume,
+                                         vol2=protImportVol2.outputVolume,
+                                         computeE=False,
+                                         masks=False,
+                                         radavg=False)
+        self.launchProtocol(protVolSubNoE)
+        self.assertIsNotNone(protVolSubNoE.outputVolume,
+                             "There was a problem with Volumes subtraction without computing energy")
+
+        protVolSubMask = self.newProtocol(XmippProtVolSubtraction,
+                                      vol1=protImportVol1.outputVolume,
+                                      vol2=protImportVol2.outputVolume,
+                                      radavg=False,
+                                      mask1=protCreateMask.outputMask,
+                                      mask2=protCreateMask.outputMask)
+        self.launchProtocol(protVolSubMask)
+        self.assertIsNotNone(protVolSubMask.outputVolume,
+                             "There was a problem with Volumes subtraction with masks")
+
+        protVolSubRadAvg = self.newProtocol(XmippProtVolSubtraction,
+                                      vol1=protImportVol1.outputVolume,
+                                      vol2=protImportVol2.outputVolume,
+                                      masks=False)
+        self.launchProtocol(protVolSubRadAvg)
+        self.assertIsNotNone(protVolSubRadAvg.outputVolume,
+                             "There was a problem with Volumes subtraction radial average")
+
+        protVolSubPdb = self.newProtocol(XmippProtVolSubtraction,
+                                      vol1=protImportVol1.outputVolume,
+                                      pdb=True,
+                                      pdbObj=protImportPdb.outputPdb,
+                                      masks=False)
+        self.launchProtocol(protVolSubPdb)
+        self.assertIsNotNone(protVolSubPdb.outputVolume,
+                             "There was a problem with Volumes subtraction pdb")
+
+        print("Run volume consensus")
+        protVolConsensus = self.newProtocol(XmippProtVolConsensus,
+                                            vols=[protImportVol1.outputVolume, protImportVol2.outputVolume,
+                                                  protVolAdj.outputVolume])
+        self.launchProtocol(protVolConsensus)
+        self.assertIsNotNone(protVolConsensus.outputVolume,
+                             "There was a problem with Volumes consensus")
diff --git a/xmipp3/tests/test_protocols_xmipp_movies.py b/xmipp3/tests/test_protocols_xmipp_movies.py
index 1f841e383..69ee1ebbd 100644
--- a/xmipp3/tests/test_protocols_xmipp_movies.py
+++ b/xmipp3/tests/test_protocols_xmipp_movies.py
@@ -32,6 +32,7 @@
 from xmipp3.protocols import *
 from pwem.protocols import ProtImportMovies, ProtImportCoordinates
 import pyworkflow.utils as pwutils
+from contextlib import redirect_stdout
 
 # Some utility functions to import movies that are used in several tests.
 class TestXmippBase(BaseTest):
@@ -327,6 +328,21 @@ def test_cct(self):
         self._checkAlignment(prot.outputMovies[1],
                              (1,7), [0, 0, 0, 0])
 
+    def test_controlPoints(self):
+        prot = self.newProtocol(XmippProtMovieCorr,
+                                doSaveMovie=False,
+                                doPSD=False,
+                                autoControlPoints=False,
+                                skipAutotuning=True,
+                                controlPointY=9,
+                                objLabel="TestControlPoints(ShouldFail)")
+        prot.inputMovies.set(self.protImport2.outputMovies)
+        with self.assertRaises(Exception,
+                               msg=("Protocol should fail because number of control points is higher "
+                                    "than number of patches for local alignment")):
+            with redirect_stdout(None):
+                self.launchProtocol(prot)
+
     def test_qbeta_SkipCrop(self):
         prot = self.newProtocol(XmippProtMovieCorr,
                                 alignFrame0=3, alignFrameN=5,
diff --git a/xmipp3/tests/test_protocols_zernike3d.py b/xmipp3/tests/test_protocols_zernike3d.py
new file mode 100644
index 000000000..2ffca6eba
--- /dev/null
+++ b/xmipp3/tests/test_protocols_zernike3d.py
@@ -0,0 +1,305 @@
+# **************************************************************************
+# *
+# * Authors:    Carlos Oscar Sorzano (coss@cnb.csic.es)
+# *             David Herreros Calero (dherreros@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 os
+import numpy as np
+
+from pyworkflow.tests import BaseTest, DataSet, setupTestProject
+
+from pwem.protocols import ProtImportVolumes, ProtImportParticles, ProtSubSet
+import pwem.emlib.metadata as md
+
+from xmipp3.protocols import (XmippProtAngularAlignmentZernike3D,
+                              XmippProtVolumeDeformZernike3D,
+                              XmippProtStructureMapZernike3D,
+                              XmippProtStructureMap)
+
+
+class TestZernike3DBase(BaseTest):
+    @classmethod
+    def setData(cls):
+        cls.dataset = DataSet(name='test_zernike3d',
+                              folder='test_zernike3d',
+                              files={'particles': 'particles/images_1720_norm.xmd',
+                                     'volumes': 'volumes/*.vol',
+                                     '1720': 'volumes/EMD-1720_norm.vol',
+                                     '1723': 'volumes/EMD-1723_norm.vol'})
+        cls.dataset = DataSet.getDataSet('test_zernike3d')
+        cls.particles = cls.dataset.getFile('particles')
+        cls.volumes = cls.dataset.getFile('volumes')
+        cls.volume_1720 = cls.dataset.getFile('1720')
+        cls.volume_1723 = cls.dataset.getFile('1723')
+        cls.clnm_pd_cpu_gold = cls.dataset.getFile('gold_standard_pd/CPU/Volumes_clnm.txt')
+        cls.clnm_pd_gpu_gold = cls.dataset.getFile('gold_standard_pd/GPU/Volumes_clnm.txt')
+
+    @classmethod
+    def runImportVolumes(cls, pattern):
+        """ Run an Import volumes protocol. """
+        protImport = cls.newProtocol(ProtImportVolumes,
+                                     filesPath=pattern,
+                                     samplingRate=1.0,
+                                     objLabel="Maps - Ribosomes 1718-1724")
+        cls.launchProtocol(protImport)
+        cls.assertIsNotNone(protImport.outputVolumes,
+                            'Import of volumes: %s, failed. outputVolumes is None.' % pattern)
+        return protImport
+
+    @classmethod
+    def runImportVolume(cls, pattern, label):
+        protImport = cls.newProtocol(ProtImportVolumes,
+                                     filesPath=pattern,
+                                     samplingRate=1.0,
+                                     objLabel="Map - Ribosome %s" % label)
+        cls.launchProtocol(protImport)
+        cls.assertIsNotNone(protImport.outputVolume,
+                            'Import of volume: %s, failed. outputVolumes is None.' % pattern)
+        return protImport
+
+    @classmethod
+    def runImportParticles(cls, pattern):
+        """ Run an Import particles protocol. """
+        print(pattern)
+        protImport = cls.newProtocol(ProtImportParticles,
+                                     objLabel='Particles - Ribosome 1720',
+                                     importFrom=ProtImportParticles.IMPORT_FROM_XMIPP3,
+                                     mdFile=pattern,
+                                     magnification=50000,
+                                     samplingRate=1.0,
+                                     haveDataBeenPhaseFlipped=True
+                                     )
+        cls.launchProtocol(protImport)
+        cls.assertIsNotNone(protImport.outputParticles,
+                            'Import of particles: %s, failed. outputParticles is None.' % pattern)
+        return protImport
+
+    @classmethod
+    def readZernikeParams(cls, filename):
+        with open(filename, 'r') as fid:
+            lines = fid.readlines()
+        basis_params = np.fromstring(lines[0].strip('\n'), sep=' ')
+        coeffs = np.fromstring(lines[1].strip('\n'), sep=' ')
+        size = int(coeffs.size / 3)
+        coeffs = np.asarray([coeffs[:size], coeffs[size:2 * size], coeffs[2 * size:]])
+        return [int(basis_params[0]), int(basis_params[1]), basis_params[2], coeffs]
+
+    @classmethod
+    def readMetadata(cls, file):
+        mdOut = md.MetaData(file)
+        return np.vstack(mdOut.getColumnValues(md.MDL_SPH_COEFFICIENTS))[:, :-8]
+
+
+class TestProtocolsZernike3D(TestZernike3DBase):
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+        TestZernike3DBase.setData()
+        cls.protImportPart = cls.runImportParticles(cls.particles)
+        cls.protImportVols = cls.runImportVolumes(cls.volumes)
+        cls.protImportVol_1720 = cls.runImportVolume(cls.volume_1720, '1720')
+        cls.protImportVol_1723 = cls.runImportVolume(cls.volume_1723, '1723')
+        cls.CORR_FILE = 'CoordinateMatrixCorr3.txt'
+        cls.EXPECTED_CORR_MAP = 'Unexpected correlation structure mapping'
+
+    def testPairAlignmentZernike3DCPU(self):
+        # Test for CPU implementation
+        alignZernike = self.newProtocol(XmippProtVolumeDeformZernike3D,
+                                        inputVolume=self.protImportVol_1720.outputVolume,
+                                        refVolume=self.protImportVol_1723.outputVolume,
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=6,
+                                        useGpu=False,
+                                        numberOfThreads=4,
+                                        objLabel="Pair Alignment Zernike3D - CPU")
+        self.launchProtocol(alignZernike)
+        self.assertIsNotNone(alignZernike.outputVolume,
+                             "There was a problem with XmippProtVolumeDeformZernike3D - CPU Version")
+        # Check coefficients and parameters are the ones expected
+        clnm_gold = self.readZernikeParams(self.clnm_pd_cpu_gold)
+        clnm_test = self.readZernikeParams(alignZernike._getExtraPath('Volumes_clnm.txt'))
+        self.assertEqual(clnm_gold[0], clnm_test[0], "There is a problem with Zernike degree")
+        self.assertEqual(clnm_gold[1], clnm_test[1], "There is a problem with Spherical Harmonics degree")
+        self.assertEqual(clnm_gold[2], clnm_test[2], "There is a problem with sphere radius")
+        diff_clnm = np.abs(np.sum(clnm_gold[3] - clnm_test[3]))
+        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+                               msg="Zernike coefficients do not have the expected value")
+
+    def testPairAlignmentZernike3DGPU(self):
+        # Test for GPU implementation
+        alignZernike = self.newProtocol(XmippProtVolumeDeformZernike3D,
+                                        inputVolume=self.protImportVol_1720.outputVolume,
+                                        refVolume=self.protImportVol_1723.outputVolume,
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=3,
+                                        useGpu=True,
+                                        numberOfThreads=1,
+                                        objLabel="Pair Alignment Zernike3D - GPU")
+        self.launchProtocol(alignZernike)
+        self.assertIsNotNone(alignZernike.outputVolume,
+                             "There was a problem with XmippProtVolumeDeformZernike3D - GPU Version")
+        # Check coefficients and parameters are the ones expected
+        clnm_gold = self.readZernikeParams(self.clnm_pd_gpu_gold)
+        clnm_test = self.readZernikeParams(alignZernike._getExtraPath('Volumes_clnm.txt'))
+        self.assertEqual(clnm_gold[0], clnm_test[0], "There is a problem with Zernike degree")
+        self.assertEqual(clnm_gold[1], clnm_test[1], "There is a problem with Spherical Harmonics degree")
+        self.assertEqual(clnm_gold[2], clnm_test[2], "There is a problem with sphere radius")
+        diff_clnm = np.abs(np.sum(clnm_gold[3] - clnm_test[3]))
+        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+                               msg="Zernike coefficients do not have the expected value")
+
+    def testStructMapZernike3DCPU(self):
+        # Test for CPU implementation
+        structureMap = self.newProtocol(XmippProtStructureMapZernike3D,
+                                        inputVolumes=[self.protImportVols.outputVolumes],
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=12.0,
+                                        useGpu=False,
+                                        numberOfThreads=4,
+                                        objLabel="StructMap Zernike3D - CPU")
+        self.launchProtocol(structureMap)
+        self.assertTrue(structureMap.isFinished(), "There was a problem with StructureMapZernike3D - CPU Version")
+        # Check if structure mappings (3D) are fine
+        sm_gold_def = np.loadtxt(self.dataset.getFile('gold_standard_sm/CPU/CoordinateMatrix3.txt'),
+                                     delimiter=' ')
+        sm_test_def = np.loadtxt(structureMap._getExtraPath('CoordinateMatrix3.txt'),
+                                     delimiter=' ')
+        diff_sm_def = np.abs(np.sum(sm_gold_def - sm_test_def))
+        self.assertAlmostEqual(diff_sm_def, 0, delta=0.2,
+                               msg="Unexpected deformation structure mapping")
+        sm_cpu_gold_corr = np.loadtxt(self.dataset.getFile('gold_standard_sm/CPU/CoordinateMatrixCorr3.txt'),
+                                     delimiter=' ')
+        sm_cpu_test_corr = np.loadtxt(structureMap._getExtraPath(self.CORR_FILE),
+                                     delimiter=' ')
+        diff_sm_corr = np.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
+        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+                               msg=self.EXPECTED_CORR_MAP)
+        sm_cpu_gold_cons = np.loadtxt(self.dataset.getFile('gold_standard_sm/CPU/ConsensusMatrix3.txt'),
+                                     delimiter=' ')
+        sm_cpu_test_cons = np.loadtxt(structureMap._getExtraPath('ConsensusMatrix3.txt'),
+                                     delimiter=' ')
+        diff_sm_cons = np.abs(np.sum(sm_cpu_gold_cons - sm_cpu_test_cons))
+        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.2,
+                               msg="Unexpected consensus structure mapping")
+
+    def testStructMapCPU(self):
+        # Test for CPU implementation
+        structureMap = self.newProtocol(XmippProtStructureMap,
+                                        inputVolumes=self.protImportVols.outputVolumes,
+                                        targetResolution=12.0,
+                                        numberOfThreads=4,
+                                        objLabel="StructMap - CPU")
+        self.launchProtocol(structureMap)
+        self.assertTrue(structureMap.isFinished(), "There was a problem with StructureMap - CPU Version")
+        # Check if structure mappings (3D) are fine
+        sm_cpu_gold_corr = np.loadtxt(self.dataset.getFile('gold_standard_sm/No_Zernike/CoordinateMatrixCorr3.txt'),
+                                     delimiter=' ')
+        sm_cpu_test_corr = np.loadtxt(structureMap._getExtraPath(self.CORR_FILE),
+                                     delimiter=' ')
+        diff_sm_corr = np.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
+        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+                               msg=self.EXPECTED_CORR_MAP)
+
+
+    def testStructMapZernike3DGPU(self):
+        # Test for GPU implementation
+        structureMap = self.newProtocol(XmippProtStructureMapZernike3D,
+                                        inputVolumes=[self.protImportVols.outputVolumes],
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=3.0,
+                                        useGpu=True,
+                                        numberOfThreads=1,
+                                        objLabel="StructMap Zernike3D - GPU")
+        self.launchProtocol(structureMap)
+        self.assertTrue(structureMap.isFinished(), "There was a problem with StructureMapZernike3D - GPU Version")
+        # Check if structure mappings (3D) are fine
+        sm_gold_def = np.loadtxt(self.dataset.getFile('gold_standard_sm/GPU/CoordinateMatrix3.txt'),
+                                     delimiter=' ')
+        sm_test_def = np.loadtxt(structureMap._getExtraPath('CoordinateMatrix3.txt'),
+                                     delimiter=' ')
+        diff_sm_def = np.abs(np.sum(sm_gold_def - sm_test_def))
+        self.assertAlmostEqual(diff_sm_def, 0, delta=0.2,
+                               msg="Unexpected deformation structure mapping")
+        sm_cpu_gold_corr = np.loadtxt(self.dataset.getFile('gold_standard_sm/GPU/CoordinateMatrixCorr3.txt'),
+                                     delimiter=' ')
+        sm_cpu_test_corr = np.loadtxt(structureMap._getExtraPath(self.CORR_FILE),
+                                     delimiter=' ')
+        diff_sm_corr = np.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
+        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+                               msg=self.EXPECTED_CORR_MAP)
+        sm_cpu_gold_cons = np.loadtxt(self.dataset.getFile('gold_standard_sm/GPU/ConsensusMatrix3.txt'),
+                                     delimiter=' ')
+        sm_cpu_test_cons = np.loadtxt(structureMap._getExtraPath('ConsensusMatrix3.txt'),
+                                     delimiter=' ')
+        diff_sm_cons = np.abs(np.sum(sm_cpu_gold_cons - sm_cpu_test_cons))
+        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.2,
+                               msg="Unexpected consensus structure mapping")
+
+    def testAngularAlignZernike3DCPU(self):
+        # Test for CPU implementation
+        alignZernike = self.newProtocol(XmippProtAngularAlignmentZernike3D,
+                                        inputParticles=self.protImportPart.outputParticles,
+                                        inputVolume=self.protImportVol_1723.outputVolume,
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=12,
+                                        useGpu=False,
+                                        numberOfMpi=4,
+                                        objLabel="Angular Align Zernike3D - CPU")
+        self.launchProtocol(alignZernike)
+        self.assertIsNotNone(alignZernike.outputParticles,
+                             "There was a problem with AngularAlignmentZernike3D - CPU Version")
+        # Check if coefficients are fine
+        aa_gold_clnm = self.readMetadata(self.dataset.getFile('gold_standard_aa/CPU/output_particles.xmd'))
+        aa_test_clnm = self.readMetadata(alignZernike._getExtraPath('output_particles.xmd'))
+        diff_clnm = np.abs(np.sum(aa_gold_clnm - aa_test_clnm))
+        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+                               msg="Unexpected Zernike coefficients")
+
+    def testAngularAlignZernike3DGPU(self):
+        # Test for GPU implementation
+        alignZernike = self.newProtocol(XmippProtAngularAlignmentZernike3D,
+                                        inputParticles=self.protImportPart.outputParticles,
+                                        inputVolume=self.protImportVol_1723.outputVolume,
+                                        l1=3,
+                                        l2=2,
+                                        targetResolution=3,
+                                        useGpu=True,
+                                        numberOfMpi=1,
+                                        objLabel="Angular Align Zernike3D - GPU")
+        self.launchProtocol(alignZernike)
+        self.assertIsNotNone(alignZernike.outputParticles,
+                             "There was a problem with AngularAlignmentZernike3D - GPU Version")
+        # Check if coefficients are fine
+        aa_gold_clnm = self.readMetadata(self.dataset.getFile('gold_standard_aa/GPU/output_particles.xmd'))
+        aa_test_clnm = self.readMetadata(alignZernike._getExtraPath('output_particles.xmd'))
+        diff_clnm = np.abs(np.sum(aa_gold_clnm - aa_test_clnm))
+        self.assertAlmostEqual(diff_clnm, 0, delta=4,
+                               msg="Unexpected Zernike coefficients")
\ No newline at end of file