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

Issue:
Center particle: Adding summary info 
Added number of centered particles and the average module of the distance for those displaced.
Also reviewed the general behavior


Discussion / hints:

Fix patch:
diff --git a/setup.py b/setup.py
index 028b96d3a..75a7bd2ae 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='22.04.0',  # Required
+    version='22.07.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 d9d368c5e..5aeba5ee4 100644
--- a/xmipp3/__init__.py
+++ b/xmipp3/__init__.py
@@ -39,7 +39,7 @@
 
 _logo = "xmipp_logo.png"
 _references = ['delaRosaTrevin2013', 'Sorzano2013']
-_currentVersion = '3.22.04.0'
+_currentVersion = '3.22.07.0'
 
 
 class Plugin(pwem.Plugin):
@@ -207,6 +207,10 @@ def installDeepLearningToolkit(plugin, env):
                                                env=plugin.getEnviron(),
                                                stdout=subprocess.PIPE
                                                ).stdout.read().decode('utf-8').split(".")[0]
+            try:
+                int(nvidiaDriverVer)
+            except ValueError:
+                nvidiaDriverVer = None
             if int(nvidiaDriverVer) < 390:
                 preMsgs.append("Incompatible driver %s" % nvidiaDriverVer)
                 cudaMsgs.append("Your NVIDIA drivers are too old (<390). "
diff --git a/xmipp3/protocols/__init__.py b/xmipp3/protocols/__init__.py
index 947db4f2d..c5d5815fb 100644
--- a/xmipp3/protocols/__init__.py
+++ b/xmipp3/protocols/__init__.py
@@ -96,6 +96,7 @@
 from .protocol_resolution_monogenic_signal import XmippProtMonoRes
 from .protocol_resolution_deepres import XmippProtDeepRes
 from .protocol_postProcessing_deepPostProcessing import XmippProtDeepVolPostProc
+from .protocol_rotate_volume import XmippProtRotateVolume
 from .protocol_rotational_spectra import XmippProtRotSpectra
 from .protocol_rotational_symmetry import XmippProtRotationalSymmetry
 from .protocol_screen_particles import XmippProtScreenParticles
@@ -138,3 +139,4 @@
 from .protocol_volume_adjust_sub import XmippProtVolAdjust, XmippProtVolSubtraction
 from .protocol_volume_consensus import XmippProtVolConsensus
 from .protocol_classes_2d_mapping import XmippProtCL2DMap
+from .protocol_deep_hand import XmippProtDeepHand
diff --git a/xmipp3/protocols/protocol_align_volume.py b/xmipp3/protocols/protocol_align_volume.py
index 44ad19269..24a62400c 100644
--- a/xmipp3/protocols/protocol_align_volume.py
+++ b/xmipp3/protocols/protocol_align_volume.py
@@ -82,7 +82,7 @@ def _defineParams(self, form):
                       help='Select the volume mask object')
         
         form.addSection(label='Search strategy')
-        form.addParam('alignmentAlgorithm', params.EnumParam, default=ALIGN_ALGORITHM_EXHAUSTIVE, 
+        form.addParam('alignmentAlgorithm', params.EnumParam, default=ALIGN_ALGORITHM_FAST_FOURIER,
                       choices=['exhaustive',
                                'local', 
                                'exhaustive + local', 
diff --git a/xmipp3/protocols/protocol_align_volume_and_particles.py b/xmipp3/protocols/protocol_align_volume_and_particles.py
index 188a0d4ca..57aba22b9 100644
--- a/xmipp3/protocols/protocol_align_volume_and_particles.py
+++ b/xmipp3/protocols/protocol_align_volume_and_particles.py
@@ -235,7 +235,7 @@ def _validate(self):
     
     def _summary(self):
         summary = []
-        summary.append("Alignment method: %s" % self.getEnumText('alignmentAlgorithm'))
+        summary.append("Alignment method: %s" % self.getEnumText('alignmentMode'))
         return summary
     
     def _methods(self):
diff --git a/xmipp3/protocols/protocol_analyze_local_ctf.py b/xmipp3/protocols/protocol_analyze_local_ctf.py
index 5ca07c749..928545bf6 100644
--- a/xmipp3/protocols/protocol_analyze_local_ctf.py
+++ b/xmipp3/protocols/protocol_analyze_local_ctf.py
@@ -27,14 +27,11 @@
 # **************************************************************************
 
 import numpy as np
-
 from pyworkflow import VERSION_2_0
 from pyworkflow.protocol.params import PointerParam
-
 from pwem import emlib
 from pwem.protocols import ProtAnalysis3D
 import pwem.emlib.metadata as md
-
 from xmipp3.convert import readSetOfMicrographs, writeSetOfMicrographs
 
 
@@ -47,7 +44,7 @@ class XmippProtAnalyzeLocalCTF(ProtAnalysis3D):
     def __init__(self, **args):
         ProtAnalysis3D.__init__(self, **args)
     
-    #--------------------------- DEFINE param functions --------------------------------------------
+    # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputMics', PointerParam, label="Input micrographs",
@@ -55,21 +52,21 @@ def _defineParams(self, form):
         form.addParam('inputSet', PointerParam, label="Input images",
                       pointerClass='SetOfParticles', help="Set of particles with assigned local defocus")
 
-    #--------------------------- INSERT steps functions --------------------------------------------
+    # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
         self._insertFunctionStep("analyzeDefocus")
         self._insertFunctionStep("createOutputStep")
 
-    #--------------------------- STEPS functions ---------------------------------------------------
+    # --------------------------- STEPS functions ---------------------------------------------------
     def analyzeDefocus(self):
         """compute R2 coefficient of each micrograph and prepare data to be later displayed in the viewer as a 3D
         representation of the distribution of particles in the micrograph"""
-        micIds=[]
-        particleIds=[]
-        x=[]
-        y=[]
-        defocusU=[]
-        defocusV=[]
+        micIds = []
+        particleIds = []
+        x = []
+        y = []
+        defocusU = []
+        defocusV = []
         for particle in self.inputSet.get():
             micIds.append(particle.getMicId())
             particleIds.append(particle.getObjId())
@@ -78,14 +75,12 @@ def analyzeDefocus(self):
             y.append(yi)
             defocusU.append(particle.getCTF().getDefocusU())
             defocusV.append(particle.getCTF().getDefocusV())
-
         uniqueMicIds = list(set(micIds))
-        self.R2={}
-
+        self.R2 = {}
         md = emlib.MetaData()
 
         for micId in uniqueMicIds:
-            idx = [i for i,j in enumerate(micIds) if j==micId]
+            idx = [i for i, j in enumerate(micIds) if j == micId]
             defocusUbyId = []
             defocusVbyId = []
             meanDefocusbyId = []
@@ -101,27 +96,34 @@ def analyzeDefocus(self):
                 ybyId.append(y[idxi])
                 particleIdsbyMicId.append(particleIds[idxi])
 
-            #defocus = c*y + b*x + a = A * X; A=[x(i),y(i)]
-            A = np.column_stack([np.ones(len(xbyId)),xbyId,ybyId])
-            polynomial, residuals, _, _ = np.linalg.lstsq(A,meanDefocusbyId,rcond=None)
+            # defocus = c*y + b*x + a = A * X; A=[x(i),y(i)]
+            A = np.column_stack([np.ones(len(xbyId)), xbyId, ybyId])
+            polynomial, _, _, _ = np.linalg.lstsq(A, meanDefocusbyId, rcond=None)
+            residuals = 0
+            for Ai, bi in zip(A, meanDefocusbyId):
+                residuali = bi - (Ai[0]*polynomial[0] + Ai[1]*polynomial[1] + Ai[2]*polynomial[2])
+                residuals += residuali*residuali
             meanDefocusbyIdArray = np.asarray(meanDefocusbyId)
             coefficients = np.asarray(polynomial)
-            self.R2[micId] = 1 - residuals / sum((meanDefocusbyIdArray - meanDefocusbyIdArray.mean()) ** 2)
+            den = sum((meanDefocusbyIdArray - meanDefocusbyIdArray.mean()) ** 2)
+            if den == 0:
+                self.R2[micId] = 0
+            else:
+                self.R2[micId] = 1 - residuals / den
             mdBlock = emlib.MetaData()
-            for xi, yi, deltafi, parti in zip(xbyId,ybyId,meanDefocusbyId,particleIdsbyMicId):
+            for xi, yi, deltafi, parti in zip(xbyId, ybyId, meanDefocusbyId, particleIdsbyMicId):
                 objId = mdBlock.addObject()
-                mdBlock.setValue(emlib.MDL_ITEM_ID,parti,objId)
-                mdBlock.setValue(emlib.MDL_XCOOR,xi,objId)
-                mdBlock.setValue(emlib.MDL_YCOOR,yi,objId)
-                mdBlock.setValue(emlib.MDL_CTF_DEFOCUSA,deltafi,objId)
+                mdBlock.setValue(emlib.MDL_ITEM_ID, parti, objId)
+                mdBlock.setValue(emlib.MDL_XCOOR, xi, objId)
+                mdBlock.setValue(emlib.MDL_YCOOR, yi, objId)
+                mdBlock.setValue(emlib.MDL_CTF_DEFOCUSA, deltafi, objId)
                 estimatedVal = coefficients[2]*yi + coefficients[1]*xi + coefficients[0]
                 residuali = deltafi - estimatedVal
-                mdBlock.setValue(emlib.MDL_CTF_DEFOCUS_RESIDUAL,residuali,objId)
-            mdBlock.write("mic_%d@%s"%(micId,self._getExtraPath("micrographDefoci.xmd")),emlib.MD_APPEND)
+                mdBlock.setValue(emlib.MDL_CTF_DEFOCUS_RESIDUAL, residuali, objId)
+            mdBlock.write("mic_%d@%s" % (micId, self._getExtraPath("micrographDefoci.xmd")), emlib.MD_APPEND)
             objId = md.addObject()
-            md.setValue(emlib.MDL_CTF_DEFOCUS_COEFS,coefficients.tolist(),objId)
-            md.write(self._getExtraPath("micrographCoef.xmd"),emlib.MD_APPEND)
-
+            md.setValue(emlib.MDL_CTF_DEFOCUS_COEFS, coefficients.tolist(), objId)
+            md.write(self._getExtraPath("micrographCoef.xmd"), emlib.MD_APPEND)
 
     def createOutputStep(self):
         """create as output a setOfParticles and add the columns of corresponding computed metadata"""
@@ -130,7 +132,7 @@ def createOutputStep(self):
         writeSetOfMicrographs(inputMicSet, fnMics)
         mdMics = md.MetaData(fnMics)
         for objId in mdMics:
-            micId = mdMics.getValue(emlib.MDL_ITEM_ID,objId)
+            micId = mdMics.getValue(emlib.MDL_ITEM_ID, objId)
             if micId in self.R2:
                 micR2 = float(self.R2[micId])
                 mdMics.setValue(emlib.MDL_CTF_DEFOCUS_R2, micR2, objId)
@@ -143,8 +145,7 @@ def createOutputStep(self):
         self._defineSourceRelation(self.inputSet, outputSet)
         self._defineSourceRelation(inputMicSet, outputSet)
 
-
-    #--------------------------- INFO functions --------------------------------------------
+    # --------------------------- INFO functions --------------------------------------------
     def _summary(self):
         summary = []
         summary.append("Local defocus analyzed for %i particles" % self.inputSet.get().getSize())
diff --git a/xmipp3/protocols/protocol_center_particles.py b/xmipp3/protocols/protocol_center_particles.py
index bd6a5ef22..61df8db3f 100644
--- a/xmipp3/protocols/protocol_center_particles.py
+++ b/xmipp3/protocols/protocol_center_particles.py
@@ -26,7 +26,7 @@
 # ******************************************************************************
 
 import numpy as np
-
+import os
 from pyworkflow import VERSION_2_0
 from pwem.constants import ALIGN_2D
 from pwem.objects import Class2D, Particle, Coordinate, Transform
@@ -85,7 +85,7 @@ def realignStep(self):
         centeredStack = md.MetaData(centeredStackName)
 
         listName = []
-        listTransform=[]
+        listTransform = []
         for rowStk in md.iterRows(centeredStack):
             listName.append(rowStk.getValue(md.MDL_IMAGE))
         for rowMd in md.iterRows(centeredMd):
@@ -104,19 +104,27 @@ def realignStep(self):
                            MD_APPEND)
 
         mdImages = md.MetaData()
-        i=0
+        i = 0
         mdBlocks = md.getBlocksInMetaDataFile(inputMdName)
         resultMat = Transform()
+        listDisplacament = []
+        centerSummary = self._getPath("summary.txt")
+        centerSummary = open(centerSummary, "w")
+        particlesCentered = 0
+        totalParticles = 0
+
         for block in mdBlocks:
             if block.startswith('class00'):
                 newMat = listTransform[i]
                 newMatrix = newMat.getMatrix()
                 mdClass = md.MetaData(block + "@" + inputMdName)
                 mdNewClass = md.MetaData()
-                i+=1
+                i += 1
+                j = 0
                 for rowIn in md.iterRows(mdClass):
                     #To create the transformation matrix (and its parameters)
                     #  for the realigned particles
+                    totalParticles += 1
                     if rowIn.getValue(md.MDL_ANGLE_PSI)!=0:
                         flag_psi=True
                     if rowIn.getValue(md.MDL_ANGLE_ROT)!=0:
@@ -138,20 +146,35 @@ def realignStep(self):
                     invResultMat = np.linalg.inv(resultMatrix)
                     centerPoint = np.dot(invResultMat,inPoint)
 
+                    if int(centerPoint[0]) > 1 or int(centerPoint[1]) > 1:
+                        particlesCentered += 1
+                        listDisplacament.append([int(centerPoint[0]), int(centerPoint[1])])
+
                     rowOut.setValue(md.MDL_XCOOR, rowOut.getValue(
                         md.MDL_XCOOR)+int(centerPoint[0]))
                     rowOut.setValue(md.MDL_YCOOR, rowOut.getValue(
                         md.MDL_YCOOR)+int(centerPoint[1]))
+
                     rowOut.addToMd(mdNewClass)
+                    j += 1
+
                 mdNewClass.write(block + "@" + self._getExtraPath(
                     'final_classes.xmd'), MD_APPEND)
                 mdImages.unionAll(mdNewClass)
+
+        listModule = [(np.sqrt((x[0] * x[0]) + (x[1] * x[1]))) for x in listDisplacament]
+        moduleDisp = round(sum(listModule) / len(listModule), 1)
+        centerSummary.write("Particles centered: {} \t({}%) "
+                            "\nAverage module displacement: {}px\n".
+                            format(particlesCentered,
+                            round((100 * particlesCentered/totalParticles ), 1),
+                            moduleDisp))
+        centerSummary.close()
         mdImages.write(self._getExtraPath('final_images.xmd'))
 
 
     def createOutputStep(self):
         inputParticles = self.inputClasses.get().getImages()
-
         outputClasses = self._createSetOfClasses2D(inputParticles)
         self._fillClasses(outputClasses)
         result = {'outputClasses': outputClasses}
@@ -160,6 +183,7 @@ def createOutputStep(self):
 
         outputParticles = self._createSetOfParticles()
         outputParticles.copyInfo(inputParticles)
+
         self._fillParticles(outputParticles)
         result = {'outputParticles': outputParticles}
         self._defineOutputs(**result)
@@ -202,7 +226,7 @@ def _fillClasses(self, outputClasses):
                 outputClasses.update(newClass)
 
     def _fillParticles(self, outputParticles):
-        """ Create the SetOfParticles and SetOfCoordinates"""
+        """ Create the SetOfParticles"""
         myParticles = md.MetaData(self._getExtraPath('final_images.xmd'))
         outputParticles.enableAppend()
         for row in md.iterRows(myParticles):
@@ -210,12 +234,20 @@ def _fillParticles(self, outputParticles):
             p = rowToParticle(row)
             outputParticles.append(p)
 
-
     # --------------------------- INFO functions -------------------------------
     def _summary(self):
         summary = []
         summary.append("Realignment of %s classes."
                        % self.inputClasses.get().getSize())
+
+        centerSummary = self._getPath("summary.txt")
+        if not os.path.exists(centerSummary):
+            summary.append("No summary file yet.")
+        else:
+            centerSummary = open(centerSummary, "r")
+            for line in centerSummary.readlines():
+                summary.append(line.rstrip())
+            centerSummary.close()
         return summary
 
     def _validate(self):
diff --git a/xmipp3/protocols/protocol_cl2d.py b/xmipp3/protocols/protocol_cl2d.py
index 6ff9e8ba5..bf6c32c76 100644
--- a/xmipp3/protocols/protocol_cl2d.py
+++ b/xmipp3/protocols/protocol_cl2d.py
@@ -407,19 +407,19 @@ def _summaryLevelFiles(self, summary, levelFiles, subset):
     def _summary(self):
         self._defineFileNames()
         summary = []
-        summary.append("Input Particles: *%d*\nClassified into *%d* classes\n"
-                       % (self.inputParticles.get().getSize(),
-                          self.numberOfClasses.get()))
-        #summary.append('- Used a _clustering_ algorithm to subdivide the original dataset into the given number of classes')
-
         levelFiles = self._getAllLevelMdFiles()
-        if levelFiles:
+
+        if not hasattr(self, 'outputClasses'):
+            summary.append("Output classes not ready yet.")
+        elif levelFiles:
             self._summaryLevelFiles(summary, levelFiles, CLASSES)
             self._summaryLevelFiles(summary, self._getAllLevelMdFiles(CLASSES_CORE), CLASSES_CORE)
             self._summaryLevelFiles(summary, self._getAllLevelMdFiles(CLASSES_STABLE_CORE), CLASSES_STABLE_CORE)
         else:
-            summary.append("Output classes not ready yet.")
-
+            summary.append("Input Particles: *%d*\nClassified into *%d* classes\n"
+                           % (self.inputParticles.get().getSize(),
+                              self.numberOfClasses.get()))
+            # summary.append('- Used a _clustering_ algorithm to subdivide the original dataset into the given number of classes')
         return summary
 
     def _methods(self):
diff --git a/xmipp3/protocols/protocol_compare_reprojections.py b/xmipp3/protocols/protocol_compare_reprojections.py
index cb4e4e611..75756be15 100644
--- a/xmipp3/protocols/protocol_compare_reprojections.py
+++ b/xmipp3/protocols/protocol_compare_reprojections.py
@@ -74,9 +74,6 @@ def _defineParams(self, form):
                       label='Evaluate residuals',
                       help='If this option is chosen, then the residual covariance matrix is calculated and '
                            'characterized. But this option takes time and disk space')
-        form.addParam('maskVol', PointerParam, pointerClass='VolumeMask', label='Volume mask', allowsNull=True,
-                      condition='doEvaluateResiduals', expertLevel=LEVEL_ADVANCED,
-                      help='3D mask for the input volume. This mask is not mandatory but advisable.')
         form.addParam('symmetryGroup', StringParam, default="c1", label='Symmetry group',
                       help='See http://xmipp.cnb.uam.es/twiki/bin/view/Xmipp/Symmetry for a description of the symmetry'
                            ' groups format. If no symmetry is present, give c1')
@@ -84,19 +81,12 @@ def _defineParams(self, form):
                       label='Angular sampling rate',
                       help='In degrees.'
                       ' This sampling defines how fine the projection gallery from the volume is explored.')
-        form.addParam('resol', FloatParam, label="Filter at resolution: ", default=3, allowsNull=True,
-                      expertLevel=LEVEL_ADVANCED,
+        form.addParam('resol', FloatParam, label="Filter at resolution: ", default=3, expertLevel=LEVEL_ADVANCED,
                       help='Resolution (A) at which subtraction will be performed, filtering the volume projections.'
                            'Value 0 implies no filtering.')
         form.addParam('sigma', FloatParam, label="Decay of the filter (sigma): ", default=3, condition='resol',
                       help='Decay of the filter (sigma parameter) to smooth the mask transition',
                       expertLevel=LEVEL_ADVANCED)
-        form.addParam('iter', IntParam, label="Number of iterations: ", default=5, expertLevel=LEVEL_ADVANCED,
-                      help='Number of iterations for the adjustment process of the images before the subtraction '
-                           'itself. Several iterations are recommended to improve the adjustment.')
-        form.addParam('rfactor', FloatParam, label="Relaxation factor (lambda): ",  expertLevel=LEVEL_ADVANCED,
-                      default=1, help='Relaxation factor for Fourier amplitude projector (POCS), it should be between '
-                                      '0 (no modification of the amplitudes of the volume 2) and 1 (no relaxation).')
         form.addParallelSection(threads=0, mpi=8)
     
     # --------------------------- INSERT steps functions --------------------------------------------
@@ -139,6 +129,7 @@ def convertStep(self, imgsFn):
     def produceResiduals(self, fnVol, fnAngles, Ts):
         anglesOutFn = self._getExtraPath("anglesCont.stk")
         projectionsOutFn = self._getExtraPath("projections.stk")
+        fnVol = self._getTmpPath("volume.vol")
         xdim = self._getDimensions()
         args = "-i %s -o %s --ref %s --optimizeAngles --optimizeShift --max_shift %d --oprojections %s --sampling %f" \
                % (fnAngles, anglesOutFn, fnVol, floor(xdim*0.05), projectionsOutFn, Ts)
@@ -147,7 +138,7 @@ def produceResiduals(self, fnVol, fnAngles, Ts):
             args += " --oresiduals %s" % self.fnResiduals
 
         if self.ignoreCTF:
-            args += " --ignoreCTF"
+            args += " --ignoreCTF "
         if self.optimizeGray:
             args += "--optimizeGray --max_gray_scale 0.95 "
         self.runJob("xmipp_angular_continuous_assign2", args)
@@ -271,44 +262,12 @@ def _useProjMatching(self):
             return False
 
     def _computeResiduals(self, fnVol):
-        vol = self.inputVolume.get().clone()
-        involdim = vol.getDim()
         if fnVol.endswith('.mrc'):
             fnVol += ':mrc'
         program = "xmipp_subtract_projection"
-        args = '-i %s --ref %s -o %s --iter %s --lambda %s --subAll' % \
-               (self.imgsFn, fnVol, self._getExtraPath("residuals"), self.iter.get(), self.rfactor.get())
-        args += ' --saveProj %s' % self._getExtraPath('')
-
-        if self.maskVol.get() is not None:
-            mskVol = self.maskVol.get().getFileName()
-        else:
-            fnDescr = self._getExtraPath("mask.descr")
-            fhDescr = open(fnDescr, 'w')
-            fhDescr.write("%d %d %d 0 \n sph + 1 0 0 0 %d" % (involdim[0], involdim[1], involdim[2], involdim[0] / 2))
-            fhDescr.close()
-            mskVol = self._getExtraPath("mask.mrc")
-            args_mask = "-i %s -o %s" % (fnDescr, mskVol)
-            self.runJob("xmipp_phantom_create", args_mask, numberOfMpi=1)
-            args_imageheader = "-i %s --sampling_rate %f" % (mskVol, vol.getSamplingRate())
-            self.runJob("xmipp_image_header", args_imageheader, numberOfMpi=1)
-        args += ' --maskVol %s' % mskVol
-
-        fnDescr2 = self._getExtraPath("mask.descr2")
-        fhDescr2 = open(fnDescr2, 'w')
-        fhDescr2.write("%d %d %d 0" % (involdim[0], involdim[1], involdim[2]))
-        fhDescr2.close()
-        mskKeep = self._getExtraPath("mask_keep.mrc")
-        args_mask_keep = "-i %s -o %s" % (fnDescr2, mskKeep)
-        self.runJob("xmipp_phantom_create", args_mask_keep, numberOfMpi=1)
-        args_imageheader2 = "-i %s --sampling_rate %f" % (mskKeep, vol.getSamplingRate())
-        self.runJob("xmipp_image_header", args_imageheader2, numberOfMpi=1)
-        args += ' --mask %s' % mskKeep
-
-        resol = self.resol.get()
-        if resol:
-            fc = vol.getSamplingRate() / resol
-            args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())
+        args = '-i %s --ref %s -o %s --save %s --max_resolution %f --sigma %d' % \
+               (self.imgsFn, fnVol, self._getExtraPath("residuals"), self._getExtraPath(''), self.resol.get(),
+                self.sigma.get())
         self.runJob(program, args, numberOfMpi=1)
         mrcsresiduals = self._getExtraPath("residuals.mrcs")
         args2 = " -i %s -o %s" % (mrcsresiduals, self.fnResiduals)
diff --git a/xmipp3/protocols/protocol_consensus_classes3D.py b/xmipp3/protocols/protocol_consensus_classes3D.py
index 86b64cc79..19726dd4c 100644
--- a/xmipp3/protocols/protocol_consensus_classes3D.py
+++ b/xmipp3/protocols/protocol_consensus_classes3D.py
@@ -1,6 +1,8 @@
 # **************************************************************************
-# *
 # * Authors:     David Maluenda (dmaluenda@cnb.csic.es)
+# *              Daniel del Hoyo Gomez (daniel.delhoyo.gomez@alumnos.upm.es)
+# *              Jorge Garcia Condado (jgcondado@cnb.csic.es)
+# *              Oier Lauzirika Zarrabeitia (oierlauzi@bizkaia.eu)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -24,184 +26,757 @@
 # *
 # **************************************************************************
 
-from pyworkflow import VERSION_2_0
+from unittest import result
+from pwem import emlib
 from pwem.protocols import EMProtocol
-from pyworkflow.protocol.params import MultiPointerParam
-from pwem.objects import Class3D
-
+from pwem.objects import Volume
+
+from pyworkflow.protocol.params import MultiPointerParam, EnumParam, IntParam
+from pyworkflow.protocol.params import Range, GE, GT, LE, LT
+from pyworkflow.object import Float
+from pyworkflow.constants import BETA
+from pyworkflow.protocol.constants import STEPS_PARALLEL
+from pyworkflow.utils.path import makePath
+
+from xmipp3.convert import setXmippAttribute
+
+import math
+import csv
+import collections
+import itertools
+import os.path
+import numpy as np
+import multiprocessing as mp
+from multiprocessing.pool import ThreadPool
 
 class XmippProtConsensusClasses3D(EMProtocol):
     """ Compare several SetOfClasses3D.
-        Return the intersection of the input classes.
+        Return the consensus clustering based on a objective function
+        that uses the similarity between clusters intersections and
+        the entropy of the clustering formed.
     """
-    _label = 'consensus classes 3D'
-    _lastUpdateVersion = VERSION_2_0
-    intersectsList = []
+    _label = 'consensus clustering 3D'
+    _devStatus = BETA
+
+    def __init__(self, *args, **kwargs):
+        EMProtocol.__init__(self, *args, **kwargs)
+        self._createFileNames()
+        
+    def _createFileNames(self):
+        iterFmt='i%(iter)06d'
+
+        myDict = {
+            'intersections': f'intersections.csv',
+            'clustering': f'clusterings/{iterFmt}.csv',
+            'objective_values': f'objective_values.csv',
+            'elbows': f'elbows.csv',
+            'reference_sizes': f'reference_sizes.csv',
+            'reference_relative_sizes': f'reference_relative_sizes.csv',
+        }
+        self._updateFilenamesDict(myDict)
 
     def _defineParams(self, form):
+        # Input data
         form.addSection(label='Input')
-
-        form.addParam('inputMultiClasses', MultiPointerParam, important=True,
-                      label="Input Classes", pointerClass='SetOfClasses3D',
-                      help='Select several sets of classes where '
-                           'to evaluate the intersections.')
+        form.addParam('inputClassifications', MultiPointerParam, pointerClass='SetOfClasses3D', label="Input classes", 
+                      important=True,
+                      help='Select several sets of classes where to evaluate the '
+                           'intersections.')
+
+        # Consensus
+        form.addSection(label='Consensus clustering')
+        form.addParam('clusteringManualCount', IntParam, label='Manual cluster count',
+                      validators=[GE(2)], default=3,
+                      help='Number of clusters to be obtained manually')
+        form.addParam('clusteringAngle', IntParam, label='Objective function threshold angle',
+                      validators=[Range(0, 90)], default=45,
+                      help='Angle to determine the cluster count')
+        
+        # Reference random classification
+        form.addSection(label='Reference random classification')
+        form.addParam('randomClassificationCount', IntParam, label='Number of random classifications',
+                      validators=[GE(0)], default=0,
+                      help='Number of random classifications used for computing the significance of the actual '
+                           'clusters. 0 for skipping this step')
+        form.addParallelSection(mpi=0, threads=mp.cpu_count())
 
     # --------------------------- INSERT steps functions -----------------------
     def _insertAllSteps(self):
-        """ Inserting one step for each intersections analisis
-        """
+        """ Inserting one step for each intersections analysis """
+        outputPrerequisites = []
+        self._setupThreadPool(int(self.numberOfThreads))
+
+        # Intersect and ensemble
+        self._insertFunctionStep('intersectStep', prerequisites=[])
+        self._insertFunctionStep('ensembleStep')
+        outputPrerequisites.append(self._insertFunctionStep('findElbowsStep'))
+
+        # Perform a random classification and use it as a reference
+        if self.randomClassificationCount.get() > 0:
+            outputPrerequisites.append(self._insertFunctionStep('checkSignificanceStep', prerequisites=[]))
+
+        # Lastly create the output    
+        self._insertFunctionStep('createOutputStep', prerequisites=outputPrerequisites)
+
+    def intersectStep(self):
+        classifications = self._getInputClassifications()
+
+        # Compute the intersection among all classes
+        intersections = self._calculateClassificationIntersections(classifications)
+
+        # Write the result to disk
+        self._writeIntersections(intersections)
+        
+    def ensembleStep(self, numClusters=1):
+        # Read input data
+        classifications = self._getInputClassifications()
+        intersections = self._readIntersections()
+
+        # Iteratively ensemble intersections until only 1 remains
+        allClusterings, allObValues = self._ensembleIntersections(classifications, intersections, numClusters)
+
+        # Reverse them so that their size increases with the index
+        allClusterings.reverse()
+        allObValues.reverse()
+
+        # Save the results
+        self._writeAllClusterings(allClusterings)
+        self._writeObjectiveValues(allObValues)
+
+    def findElbowsStep(self):
+        obValues = self._readObjectiveValues()
+        numClusters = list(range(1, 1+len(obValues)))
+
+        # Normalize clusters and obValues
+        normObValues = self._normalizeValues(obValues)
+        normNumClusters = self._normalizeValues(numClusters)
+
+        # Calculate the elbows
+        elbows = {
+            'origin': self._findElbowOrigin(normNumClusters, normObValues),
+            'angle': self._findElbowAngle(normNumClusters, normObValues, -np.radians(self.clusteringAngle.get())),
+            'pll': self._findElbowPll(obValues)
+        }
+
+        # Start counting at 1
+        for key in elbows.keys():
+            elbows[key] += 1
+
+        # Save the results
+        self._writeElbows(elbows)
+
+    def checkSignificanceStep(self):
+        classifications = self._getInputClassifications()
+        numRand = self.randomClassificationCount.get()
+        threadPool = self._getThreadPool()
+
+        # Compute the random consensus
+        consensusSizes, consensusRelativeSizes = self._calculateReferenceClassification(
+            classifications, numRand, threadPool
+        )
+
+        # Sort the sizes so that percentile calculation is easy
+        consensusSizes.sort()
+        consensusRelativeSizes.sort()
+
+        # Write results to disk
+        self._writeReferenceClassificationSizes(consensusSizes)
+        self._writeReferenceClassificationRelativeSizes(consensusRelativeSizes)
 
-        self._insertFunctionStep('compareFirstStep', 
-                                 self.inputMultiClasses[0].get().getObjId(),
-                                 self.inputMultiClasses[1].get().getObjId())
+    def createOutputStep(self):
+        # Read all the necessary data from disk
+        manualClusterCount = self.clusteringManualCount.get()
+        initialClusterCount = self._readClusteringCount()
+        elbows = self._readElbows()
+        consensusSizes = self._readReferenceClassificationSizes()
+        consensusRelativeSizes = self._readReferenceClassificationRelativeSizes()
+
+        # Append the manual and initial indices to the elbow
+        if manualClusterCount > 0:
+            elbows['manual'] = manualClusterCount
+        elbows['initial'] =  initialClusterCount
+
+        # Read selected clusters
+        clusterings = dict(zip(elbows.keys(), map(self._readClustering, elbows.values())))
+
+        # Create the output classes and define them
+        particles = self._getInputParticles()
+        outputClasses = dict(zip(
+            map(lambda name : 'outputClasses_'+name,
+                clusterings.keys() 
+            ),
+            map(self._createOutputClasses3D,
+                itertools.repeat(particles),
+                clusterings.values(),
+                clusterings.keys(),
+                itertools.repeat(consensusSizes),
+                itertools.repeat(consensusRelativeSizes)
+            )
+        ))
+        self._defineOutputs(**outputClasses)
+
+        # Stablish source output relationships
+        sources = list(self.inputClassifications)
+        destinations = list(outputClasses.values())
+        for src in sources:
+            for dst in destinations:
+                self._defineSourceRelation(src, dst)
 
-        if len(self.inputMultiClasses)>2:
-            for i in range(2,len(self.inputMultiClasses)):
-                self._insertFunctionStep('compareOthersStep', i,
-                                     self.inputMultiClasses[i].get().getObjId())
+    # --------------------------- INFO functions -------------------------------
+    def _validate(self):
+        errors = []
 
-        self._insertFunctionStep('createOutputStep')
+        # Ensure that all classifications are made of the same images
+        images = self._getInputParticles(0)
+        for i in range(1, len(self.inputClassifications)):
+            if self._getInputParticles(i) != images:
+                errors.append(f'Classification {i} has been done with different images')
 
-    def compareFirstStep(self, objId1, objId2):
-        """ We found the intersections for the two firsts sets of classes
-        """
-        set1Id = 0
-        set2Id = 1
-        set1 = self.inputMultiClasses[set1Id].get()
-        set2 = self.inputMultiClasses[set2Id].get()
+        return errors
 
-        print('Computing intersections between classes form set %s and set %s:'
-               % (set1.getNameId(), set2.getNameId()))
+    # --------------------------- UTILS functions ------------------------------
+    def _getThreadPool(self):
+        return getattr(self, '_threadPool', None)
+
+    def _setupThreadPool(self, maxThreads):
+        if maxThreads > 1:
+            self._threadPool = ThreadPool(processes=maxThreads)
+
+    def _getInputParticles(self, i=0):
+        return self.inputClassifications[i].get().getImages()
+
+    def _getInputClassifications(self):
+        if not hasattr(self, '_classifications'):
+            self._classifications = self._convertInputClassifications(self.inputClassifications)
+        
+        return self._classifications
+
+    def _writeTable(self, path, table, fmt='%s'):
+        np.savetxt(path, table, delimiter=',', fmt=fmt)
+
+    def _readTable(self, path, dtype=str):
+        return np.genfromtxt(path, delimiter=',', dtype=dtype)
+
+    def _writeList(self, path, lst, fmt='%s'):
+        self._writeTable(path, lst, fmt)
+
+    def _readList(self, path, dtype=str):
+        return self._readTable(path, dtype=dtype)
+        
+    def _writeClassification(self, path, classification):
+        with open(path, 'w') as file:
+            writer = csv.writer(file)
+            
+            # Write the header
+            header = ['Representative index', 'Representative path', 'Source size', 'Particle ids']
+            writer.writerow(header)
+
+            # Transform the classification
+            def f(cls):
+                result = []
+                
+                # Add the location of the representative
+                representative = cls.getRepresentative()
+                location = representative.getLocation() if representative is not None else (0, '')
+                result += list(location)
 
-        newList = []
-        for cls1 in set1:
-            cls1Id = cls1.getObjId()
-            ids1 = cls1.getIdSet()
+                # Add the source size
+                result += [cls.getSourceSize()]
 
-            for cls2 in set2:
-                cls2Id = cls2.getObjId()
-                ids2 = cls2.getIdSet()
+                # Add the particles
+                result += list(cls.getParticleIds())
 
-                interTuple = self.intersectClasses(set1Id, cls1Id, ids1,
-                                                   set2Id, cls2Id, ids2)
+                assert(len(result) == (len(cls) + 3))
+                return result
+            writer.writerows(map(f, classification))
 
-                newList.append(interTuple)
+    def _readClassification(self, path):
+        with open(path, 'r') as file:
+            reader = csv.reader(file)
 
-        self.intersectsList = newList
+            # Ignore the header
+            next(reader)
 
-    def compareOthersStep(self, set1Id, objId):
-        """ We found the intersections for the two firsts sets of classes
-        """
-        set1 = self.inputMultiClasses[set1Id].get()
+            # Transform each row to a ParticleCluster in a classification
+            def f(row):
+                ite = iter(row)
 
-        print('Computing intersections between classes form set %s and '
-              'the previous ones:' % (set1.getNameId()))
+                # Parse the representative
+                index = int(next(ite))
+                path = next(ite)
+                representative = Volume(location=(index, path))
 
-        newList = []
-        currDB = self.intersectsList
-        for cls1 in set1:
-            cls1Id = cls1.getObjId()
-            ids1 = cls1.getIdSet()
+                # Parse the source size
+                sourceSize = int(next(ite))
 
-            for currTuple in currDB:
-                ids2 = currTuple[1]
-                set2Id = currTuple[2]
-                cls2Id = currTuple[3]
-                clSize = currTuple[4]
+                # Parse the id set
+                particleIds = set(map(lambda id : int(id), ite))
 
-                interTuple = self.intersectClasses(set1Id, cls1Id, ids1,
-                                                   set2Id, cls2Id, ids2, clSize)
+                result = XmippProtConsensusClasses3D.ParticleCluster(particleIds, representative, sourceSize)
+                assert(len(row) == (len(result) + 3))
+                return result
+            return list(map(f, reader))
 
-                newList.append(interTuple)
-                
-        self.intersectsList = newList
+    def _writeDictionary(self, path, d):
+        with open(path, 'w') as file:
+            writer = csv.writer(file)
+            writer.writerow(d.keys()) # Header
+            writer.writerow(d.values()) # Values
 
-    def createOutputStep(self):
+    def _readDictionary(self, path, dtype=str):
+        with open(path, 'r') as file:
+            reader = csv.reader(file)
+            keys = list(next(reader))
+            values = list(map(lambda value : dtype(value), next(reader)))
+            return dict(zip(keys, values))
 
-        # self.intersectsList.sort(key=lambda e: e[0], reverse=True)
+    def _writeIntersections(self, intersections):
+        self._writeClassification(self._getExtraPath(self._getFileName('intersections')), intersections)
 
-        # print("   ---   S O R T E D:   ---")
-        # numberOfPart = 0
-        # for classItem in self.intersectsList:
-        #     printTuple = (classItem[0], classItem[2])
-        #     print(printTuple)
-        #     numberOfPart += classItem[0]
+    def _readIntersections(self):
+        return self._readClassification(self._getExtraPath(self._getFileName('intersections')))
 
-        # print('Number of intersections: %d' % len(self.intersectsList))
-        # print('Total of particles: %d' % numberOfPart)
+    def _writeClustering(self, iter, clustering):
+        self._writeClassification(self._getExtraPath(self._getFileName('clustering', iter=iter)), clustering)
 
-        inputParticles = self.inputMultiClasses[0].get().getImages()
-        outputClasses = self._createSetOfClasses3D(inputParticles)
+    def _readClustering(self, iter):
+        return self._readClassification(self._getExtraPath(self._getFileName('clustering', iter=iter)))
 
-        for classItem in self.intersectsList:
-            numOfPart = classItem[0]
-            partIds = classItem[1]
-            setRepId = classItem[2]
-            clsRepId = classItem[3]
+    def _readClusteringCount(self):
+        count = 0
+        while os.path.exists(self._getExtraPath(self._getFileName('clustering', iter=count+1))):
+            count += 1
+        return count
 
-            setRep = self.inputMultiClasses[setRepId].get()
-            clRep = setRep[clsRepId]
+    def _writeAllClusterings(self, clusterings):   
+        # Create the path
+        makePath(os.path.dirname(self._getExtraPath(self._getFileName('clustering', iter=0))))
 
-            newClass = Class3D()
-            # newClass.copyInfo(clRep)
-            newClass.setAcquisition(clRep.getAcquisition())
-            newClass.setRepresentative(clRep.getRepresentative())
-            newClass.setSamplingRate(clRep.getSamplingRate())
+        # Write all
+        for iter, clustering in enumerate(clusterings, start=1):
+            assert(len(clustering) == iter)
+            self._writeClustering(iter, clustering)
 
-            outputClasses.append(newClass)
+    def _readAllClusterings(self):
+        return list(map(self._readClustering, range(1, 1+self._readClusteringCount())))
 
-            enabledClass = outputClasses[newClass.getObjId()]
-            enabledClass.enableAppend()
-            for itemId in partIds:
-                enabledClass.append(inputParticles[itemId])
+    def _writeObjectiveValues(self, allObValues):
+        self._writeList(self._getExtraPath(self._getFileName('objective_values')), allObValues, '%f')
 
-            outputClasses.update(enabledClass)
+    def _readObjectiveValues(self):
+        return self._readList(self._getExtraPath(self._getFileName('objective_values')), dtype=float)
 
-        self._defineOutputs(outputClasses=outputClasses)
-        for item in self.inputMultiClasses:
-            self._defineSourceRelation(item, outputClasses)
+    def _writeElbows(self, elbows):
+        self._writeDictionary(self._getExtraPath(self._getFileName('elbows')), elbows)
 
+    def _readElbows(self):
+        return self._readDictionary(self._getExtraPath(self._getFileName('elbows')), dtype=int)
 
-    # --------------------------- INFO functions -------------------------------
-    def _summary(self):
-        summary = []
-        return summary
+    def _writeReferenceClassificationSizes(self, sizes):
+        self._writeList(self._getExtraPath(self._getFileName('reference_sizes')), sizes, '%d')
 
-    def _methods(self):
-        methods = []
-        return methods
+    def _readReferenceClassificationSizes(self):
+        path = self._getExtraPath(self._getFileName('reference_sizes'))
+        if os.path.exists(path):
+            return self._readList(path, dtype=int)
+        else:
+            return None
 
-    def _validate(self):
-        errors = [] if len(self.inputMultiClasses)>1 else \
-                 ["More than one Input Classes is needed to compute the consensus."]
-        return errors
+    def _writeReferenceClassificationRelativeSizes(self, sizes):
+        self._writeList(self._getExtraPath(self._getFileName('reference_relative_sizes')), sizes, '%f')
 
+    def _readReferenceClassificationRelativeSizes(self):
+        path = self._getExtraPath(self._getFileName('reference_relative_sizes'))
+        if os.path.exists(path):
+            return self._readList(path, dtype=float)
+        else:
+            return None
+
+    def _convertInputClassifications(self, classifications):
+        """ Returns the list of lists of sets that stores the set of ids of each class"""
+        def class3dToParticleCluster(cls):
+            ids = cls.getIdSet()
+            rep = cls.getRepresentative().clone() if cls.hasRepresentative() else None
+            return XmippProtConsensusClasses3D.ParticleCluster(ids, rep)
+
+        f = lambda classification : list(map(class3dToParticleCluster, classification.get()))
+        result = list(map(f, classifications))
+        return result
+    
+    def _calculateClassificationIntersections(self, classifications):
+        # Start with the first classification
+        result = classifications[0]
+
+        for i in range(1, len(classifications)):
+            classification = classifications[i]
+            intersections = []
+
+            # Perform the intersection with previous classes
+            for cls0 in result:
+                for cls1 in classification:
+                    intersection = cls0.intersection(cls1)
+                    if intersection:
+                        intersections.append(intersection)
+
+            # Use the new intersection for the next step
+            result = intersections
+
+        return result
+
+    def _calculateSubsetSimilarity(self, n, c):
+        """ Return the similarity of intersection subset (n) with class subset (c) """
+        inter = len(n.intersection(c))
+        return inter / len(n)
+
+    def _buildSimilarityVector(self, n, cs):
+        """ Build the similarity vector of n subset with each c subset in cs """
+        return list(map(self._calculateSubsetSimilarity, itertools.repeat(n), cs))
+
+    def _calculateEntropy(self, p):
+        """ Return the entropy of the elements in a vector """
+        result = -sum(map(lambda x : x*math.log2(x), p))
+        assert(result >= 0) # Entropy cannot be negative
+        return result
+
+    def _calculateDistribution(self, ns):
+        """ Return the vector of distribution (p) from the n subsets in ns
+        """
+        p = list(map(len, ns))
+        return np.array(p) / sum(p)
 
-    # --------------------------- UTILS functions ------------------------------
-    def intersectClasses(self, setId1, clId1, ids1,
-                               setId2, clId2, ids2, clsSize2=None):
-        size1 = len(ids1)
-        size2 = len(ids2) if clsSize2 is None else clsSize2
+    def _calculateCosineSimilarity(self, v1, v2):
+        """ Return the cosine similarity of two vectors that is used as similarity """
+        return np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
 
-        inter = ids1.intersection(ids2)
+    def _objectiveFunction(self, ns, indices, vs, eps=0.01):
+        """Objective function to minimize based on the entropy and the similarity between the clusters to merge
+        """
+        p = self._calculateDistribution(ns)
+        pi = self._calculateDistribution(self._ensembleIntersectionPair(ns, indices))
+        h = self._calculateEntropy(p)
+        hi = self._calculateEntropy(pi)
+
+        sv = self._calculateCosineSimilarity(vs[indices[0]], vs[indices[1]])
+
+        return (h - hi) / (sv + eps)
+
+    def _ensembleIntersectionPair(self, intersections, indices):
+        """ Return a list of subsets where the intersections[indices] have been merged """
+        # Divide the input intersections into the ones referred by indices and not
+        selection = [intersections[i] for i in indices]  # intersection[indices]
+        result = [intersections[i] for i in range(len(intersections)) if i not in indices]  # intersections - selection
+
+        # Merge all the selected clusters into the same one
+        merged = XmippProtConsensusClasses3D.ParticleCluster.union(*selection)
+        
+        # Add the merged items to the result
+        result.append(merged)
+        return result
+
+    def _ensembleCheapestIntersectionPair(self, allClasses, intersections):
+        """ Merges the lest costly pair of intersections """
+        # Compute the similarity vectors
+        similarityVectors = list(map(self._buildSimilarityVector, intersections, itertools.repeat(allClasses)))
+
+        # For each possible pair of intersections compute the cost of merging them
+        obValues = {}
+        for pair in itertools.combinations(range(len(intersections)), r=2):
+            obValues[pair] = self._objectiveFunction(intersections, pair, similarityVectors)
+
+        values = list(obValues.values())
+        keys = list(obValues.keys())
+
+        # Select the minimum cost pair and merge them
+        obValue = min(values)
+        mergePair = keys[values.index(obValue)]
+        mergedIntersections = self._ensembleIntersectionPair(intersections, mergePair)
+
+        return mergedIntersections, obValue
+
+    def _ensembleIntersections(self, classifications, intersections, numClusters=1):
+        """ Ensemble clustering by COSS that merges interactions of clusters based
+            on entropy minimization """
+        # Results
+        allIntersections = [intersections]
+        allObValues = [0.0]
+
+        # Obtain a list with all source classes
+        allClasses = list(itertools.chain(*classifications))
+
+        # Iteratively merge intersections
+        while len(allIntersections[-1]) > numClusters:
+            mergedIntersections, obValue = self._ensembleCheapestIntersectionPair(allClasses, allIntersections[-1])
+            allIntersections.append(mergedIntersections)
+            allObValues.append(obValue)
+
+        assert(len(allIntersections) == len(allObValues))
+        return allIntersections, allObValues
+            
+    def _normalizeValues(self, x):
+        """Normalize values to range 0 to 1"""
+        return (x-np.min(x))/(np.max(x)-np.min(x))
+
+    def _calculateGaussianMleEstimates(self, d, q):
+        """ MLE estimates of guassian distributions
+        with common variance parameter
+        """
+        p = len(d)
+        d1 = d[:q]
+        d2 = d[q:]
+        mu1 = np.mean(d1)
+        mu2 = np.mean(d2)
+        var1 = np.var(d1)
+        var2 = np.var(d2)
+        sigma = ((q-1)*var1+(p-q-1)*var2)/(p-2)
+        theta1 = [mu1, sigma]
+        theta2 = [mu2, sigma]
+        return theta1, theta2
+
+    def _calculateGaussianDistribution(self, d, theta):
+        """ Gaussian pdf """
+        m = theta[0]
+        s = theta[1]
+        gauss = 1/(np.sqrt(2*np.pi*s))*np.exp(-1/(2*s)*(d-m)**2)
+        return gauss
+
+    def _calculateProfileLogLikelihood(self, d, q, theta1, theta2):
+        """ Profile log likelihood for given parameters """
+        log_f_q = np.log(self._calculateGaussianDistribution(d[:q], theta1))
+        log_f_p = np.log(self._calculateGaussianDistribution(d[q:], theta2))
+        l_q = np.sum(log_f_q) + np.sum(log_f_p)
+        return l_q
+
+    def _calculateFullProfileLogLikelihood(self, d):
+        """ Calculate profile log likelihood for each partition
+            of the data """
+        pll = []
+        for q in range(1, len(d)):
+            theta1, theta2 = self._calculateGaussianMleEstimates(d, q)
+            pll.append(self._calculateProfileLogLikelihood(d, q, theta1, theta2))
+        return pll
+
+    def _findElbowPll(self, obValues):
+        """ Find the elbow according to the full profile likelihood
+        """
+        # Get profile log likelihood for log of objective values
+        # Remove last obValue as it is zero and log(0) is undefined
+        pll = self._calculateFullProfileLogLikelihood(np.log(obValues[:-1]))
+        return np.argmax(pll)
 
-        if size1 < size2:
-            setId = setId1
-            clsId = clId1
-            clsSize = size1
+    def _findElbowOrigin(self, x, y):
+        """ Find the point closest to origin from normalied data
+        """
+        coors = np.column_stack((x, y))
+        distances = np.linalg.norm(coors, axis=1)
+        assert(len(x) == len(distances))
+        return np.argmin(distances)
+
+    def _findElbowAngle(self, x, y, angle=-math.pi/4):
+        """ Find the angle the slope of the function makes and
+            return the point at which the slope crosses the given angle
+        """
+        dx = np.diff(x)
+        dy = np.diff(y)
+        angles = np.arctan2(dy, dx)
+
+        # Find the first point with an angle greater than the given one
+        crossing = np.argmax(angles >= angle)
+
+        return crossing
+
+    def _getClassificationLengths(self, classifications):
+        """ Returns a list of lists that stores the lengths of each classification """
+        return list(map(lambda classification : list(map(len, classification)), classifications))
+
+    def _calculateRandomClassification(self, C, N):
+        """ Randomly classifies N element indices into groups with
+        sizes defined by rows of C """
+        Cp = []
+
+        for Ci in C:
+            assert(sum(Ci) == N)  # The groups should address all the elements TODO: Maybe LEQ?
+            x = np.argsort(np.random.uniform(size=N)).tolist()  # Shuffles a [0, N) iota
+
+            # Select the number random indices requested by each group
+            Cip = []
+            first = 0
+            for s in Ci:
+                Cip.append(XmippProtConsensusClasses3D.ParticleCluster(x[first:first + s]))
+                first += s
+
+            # Add the random classification to the result
+            Cp.append(Cip)
+
+        return Cp
+
+    def _calculateRandomClassificationConsensus(self, C, N):
+        """ Obtains the intersections of a consensus
+            of a random classification of sizes defined in C of
+            N elements. C is a list of lists, where de sum of each
+            row must equal N """
+
+        # Create random partitions of same size
+        randomClassification = self._calculateRandomClassification(C, N)
+
+        # Compute the repeated classifications
+        return self._calculateClassificationIntersections(randomClassification)
+
+    def _calculateReferenceClassification(self, classifications, numExec, threadPool=None):
+        """ Create random partitions of same size to compare the quality
+         of the classification """
+        clusterLengths = self._getClassificationLengths(classifications)
+        numParticles = sum(clusterLengths[0])
+
+        # Repeatedly obtain a consensus of a random classification of same size
+        if threadPool:
+            consensus = threadPool.starmap(
+                self._calculateRandomClassificationConsensus,
+                itertools.repeat((clusterLengths, numParticles), numExec)
+            )
         else:
-            setId = setId2
-            clsId = clId2
-            clsSize = size2
-
-        # print(" ")
-        # print(" - Intersection of cl%d of set%d (%d part.) and "
-        #                          "cl%d of set%d (%d part.):"
-        #        % (clId1, setId1, len(ids1), clId2, setId2, len(ids2)))
-        # print("    Size1=%d < Size2=%d = %s" 
-        #        % (size1, size2, size1<size2))
-        # print("      -> from set %d calss %d, with %d part. in the intersection." 
-        #        % (setId, clsId, len(inter)))
-        # print(" -  -  -  -  -  -  -  -  -  -")
-
-        return (len(inter), inter, setId, clsId, clsSize)
+            consensus = list(map(
+                self._calculateRandomClassificationConsensus, 
+                itertools.repeat(clusterLengths, numExec),
+                itertools.repeat(numParticles, numExec)
+            ))
+        assert(len(consensus) == numExec)
+
+        # Concatenate all consensuses
+        consensus = list(itertools.chain(*consensus))
+
+        # Obtain the consensus sizes
+        consensusSizes = list(map(len, consensus))
+        consensusSizeRatios = list(map(XmippProtConsensusClasses3D.ParticleCluster.getRelativeSize, consensus))
+        return consensusSizes, consensusSizeRatios
+
+    def _createOutputClasses3D(self, particles, clustering, name, randomConsensusSizes=None, randomConsensusRelativeSizes=None):
+        outputClasses = self._createSetOfClasses3D(particles, suffix=name)
+
+        # Create a list with the filenames of the representatives
+        representatives = list(map(
+            XmippProtConsensusClasses3D.ParticleCluster.getRepresentative,
+            clustering
+        ))
+
+        # Fill the output
+        loader = XmippProtConsensusClasses3D.ClassesLoader(
+            clustering, 
+            representatives,
+            randomConsensusSizes,
+            randomConsensusRelativeSizes
+        )
+        loader.fillClasses(outputClasses)
+
+        return outputClasses
+    
+    class ParticleCluster:
+        """ Keeps track of the information related to successive class intersections.
+            It is instantiated with a single class and allows to perform intersections
+            and unions with it"""
+        def __init__(self, particleIds, representative=None, sourceSize=None):
+            self._particleIds = set(particleIds)  # Particle ids belonging to this intersection
+            self._representative = representative
+            self._sourceSize = sourceSize if sourceSize is not None else len(self._particleIds)  # The size of the representative class
+
+        def __len__(self):
+            return len(self.getParticleIds())
+
+        def __eq__(self, other):
+            return self.getParticleIds() == other.getParticleIds()
+
+        def intersection(self, *others):
+            return self._combine(set.intersection, min, XmippProtConsensusClasses3D.ParticleCluster.getSourceSize, *others)
+
+        def union(self, *others):
+            return self._combine(set.union, max, len, *others)
+
+        def getParticleIds(self):
+            return self._particleIds
+
+        def getRepresentative(self):
+            return self._representative
+
+        def getSourceSize(self):
+            return self._sourceSize
+
+        def getRelativeSize(self):
+            return len(self) / self.getSourceSize()
+
+        def _combine(self, operation, selector, selectionCriteria, *others):
+            allItems = (self, ) + others
+            
+            # Perform the requested operation among the id sets
+            allParticleIds = map(XmippProtConsensusClasses3D.ParticleCluster.getParticleIds, allItems)
+            particleIds = operation(*allParticleIds)
+
+            # Select the class with the appropriate criteria
+            selection = selector(allItems, key=selectionCriteria)
+            representative = selection.getRepresentative()
+            sourceSize = selection.getSourceSize()
+
+            return XmippProtConsensusClasses3D.ParticleCluster(particleIds, representative, sourceSize)
+
+    class ClassesLoader:
+        """ Helper class to produce classes
+        """
+        def __init__(self, clustering, representatives, randomConsensusSizes=None, randomConsensusRelativeSizes=None):
+            self.clustering = clustering
+            self.classification = self._createClassification(clustering)
+            self.representatives = representatives
+            self.randomConsensusSizes = randomConsensusSizes
+            self.randomConsensusRelativeSizes = randomConsensusRelativeSizes
+
+        def fillClasses(self, clsSet):
+            clsSet.classifyItems(updateItemCallback=self._updateParticle,
+                                updateClassCallback=self._updateClass,
+                                itemDataIterator=iter(self.classification.items()),
+                                doClone=False)
+
+        def _updateParticle(self, item, row):
+            particleId, classId = row
+            assert(item.getObjId() == particleId)
+            item.setClassId(classId)
+
+        def _updateClass(self, item):
+            classId = item.getObjId()
+            classIdx = classId-1
+
+            # Set the representative
+            item.setRepresentative(self.representatives[classIdx])
+
+            # Set the size p-value
+            if self.randomConsensusSizes is not None:
+                size = len(self.clustering[classIdx])
+                percentile = self._findPercentile(self.randomConsensusSizes, size)
+                pValue = 1 - percentile
+                setXmippAttribute(item, emlib.MDL_CLASS_INTERSECTION_SIZE_PVALUE, Float(pValue))
+
+
+            # Set the relative size p-value
+            if self.randomConsensusRelativeSizes is not None:
+                size = self.clustering[classIdx].getRelativeSize()
+                percentile = self._findPercentile(self.randomConsensusRelativeSizes, size)
+                pValue = 1 - percentile
+                setXmippAttribute(item, emlib.MDL_CLASS_INTERSECTION_RELATIVE_SIZE_PVALUE, Float(pValue))
+
+        def _createClassification(self, clustering):
+            # Fill the classification data (particle-class mapping)
+            result = {}
+            for cls, data in enumerate(clustering):
+                for particleId in data.getParticleIds():
+                    result[particleId] = cls + 1 
+
+            return collections.OrderedDict(sorted(result.items())) 
+
+        def _findPercentile(self, data, value):
+            """ Given an array of values (data), finds the corresponding percentile of the value.
+                Percentile is returned in range [0, 1] """
+            assert(np.array_equal(sorted(data), data))  # In order to iterate it in ascending order
+
+            # Count the number of elements that are smaller than the given value
+            i = 0
+            while i < len(data) and data[i] < value:
+                i += 1
+
+            # Convert it into a percentage
+            return float(i) / float(len(data))
diff --git a/xmipp3/protocols/protocol_ctf_consensus.py b/xmipp3/protocols/protocol_ctf_consensus.py
index e4bae3a74..6f81c5643 100644
--- a/xmipp3/protocols/protocol_ctf_consensus.py
+++ b/xmipp3/protocols/protocol_ctf_consensus.py
@@ -189,7 +189,7 @@ def _defineParams(self, form):
                            'If *No*, only the primary metadata (plus consensus '
                            'scores) will be in the resulting CTF.')
 
-        form.addParallelSection(threads=4, mpi=1)
+        form.addParallelSection(threads=1, mpi=1)
 
 # --------------------------- INSERT steps functions -------------------------
     def _insertAllSteps(self):
diff --git a/xmipp3/protocols/protocol_deep_hand.py b/xmipp3/protocols/protocol_deep_hand.py
new file mode 100644
index 000000000..b981952c8
--- /dev/null
+++ b/xmipp3/protocols/protocol_deep_hand.py
@@ -0,0 +1,169 @@
+# **************************************************************************
+# *
+# * Authors:     Jorge Garcia Condado (jorgegarciacondado@gmail.com)
+# *
+# * 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 EMProtocol
+from pwem.emlib.image import ImageHandler
+from pwem.objects import Volume
+
+from pyworkflow.protocol.params import PointerParam, FloatParam, IntParam
+from pyworkflow.object import Float
+from pyworkflow.utils.path import cleanPath
+
+from xmipp3.base import XmippProtocol
+from xmipp3.convert import getImageLocation
+
+class XmippProtDeepHand(EMProtocol, XmippProtocol):
+    """Protocol to returns handedness of structure from trained deep learning model
+    """
+
+    _label ="deep hand"
+    _conda_env = "xmipp_deepHand"
+
+    def __init__(self, *args, **kwargs):
+        EMProtocol.__init__(self, *args, **kwargs)
+        XmippProtocol.__init__(self)
+
+    def _defineParams(self, form):
+
+        form.addSection('Input')
+        form.addParam('inputVolume', PointerParam, pointerClass="Volume",
+                      label='Input Volume', allowsNull=False,
+                      important=True, help="Volume to process")
+        form.addParam('threshold', FloatParam, label='Mask Threshold',
+                      allowsNull=False, important=True,  help="Threshold for mask creation")
+        form.addParam('thresholdAlpha', FloatParam, label='Alpha Threshold',
+                      default=0.7, help="Threshold for alpha helix determination")
+        form.addParam('thresholdHand', FloatParam, label='Hand Threshold',
+                      default=0.6, help="Hand threshold to flip volume")
+
+# --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('preprocessStep')
+        self._insertFunctionStep('predictStep')
+        self._insertFunctionStep('flipStep')
+        self._insertFunctionStep('createOutputStep')
+
+    def preprocessStep(self):
+
+        # Get volume information
+        volume = self.inputVolume.get()
+        fn_vol = getImageLocation(volume)
+        T_s = volume.getSamplingRate()
+
+        # Paths to new files created
+        self.resizedVolFile = self._getPath('resizedVol.mrc')
+        self.maskFile = self._getPath('mask.mrc')
+        self.filteredVolFile = self._getPath('filteredVol.mrc')
+
+        # Resize to 1A/px
+        self.runJob("xmipp_image_resize", "-i %s -o %s --factor %f" %
+                    (fn_vol, self.resizedVolFile, T_s))
+
+        # Threshold to obtain mask
+        self.runJob("xmipp_transform_threshold",
+                    "-i %s -o %s --select below %f --substitute binarize"
+                    % (self.resizedVolFile, self.maskFile, self.threshold.get()))
+
+        # Filter to 5A
+        self.runJob("xmipp_transform_filter", "-i %s -o %s "\
+                    "--fourier low_pass %f --sampling 1"
+                    % (self.resizedVolFile, self.filteredVolFile, 5.0))
+
+    def predictStep(self):
+
+        # Get saved models DTlK
+        alpha_model= self.getModel('deepHand', '5A_SSE_experimental.pth')
+        hand_model = self.getModel('deepHand', '5A_TL_hand_alpha.pth')
+
+        # Predict hand
+        args = "--alphaModel %s --handModel %s -o %s " \
+               "--alphaThr %f --pathVf %s --pathVmask %s" % (
+                alpha_model, hand_model, self._getExtraPath(),
+                self.thresholdAlpha.get(), self._getPath('filteredVol.mrc'),
+                self._getPath('mask.mrc'))
+        self.runJob("xmipp_deep_hand", args, env=self.getCondaEnv())
+
+        # Store hand value
+        hand_file = self._getExtraPath('hand.txt')
+        f = open(hand_file)
+        self.hand = Float(float(f.read()))
+        f.close()
+
+    def flipStep(self):
+        if self.hand.get() > self.thresholdHand.get():
+            volume = self.inputVolume.get()
+            fn_vol = getImageLocation(volume)
+            pathFlipVol = self._getPath('flipVol.mrc')
+            self.runJob("xmipp_transform_mirror", "-i %s -o %s --flipX" \
+                            % (fn_vol, pathFlipVol))
+
+    def createOutputStep(self):
+        self._defineOutputs(outputHand=self.hand)
+
+        if self.hand.get() > self.thresholdHand.get():
+            vol = Volume()
+            volFile = self._getPath('flipVol.mrc')
+            vol.setFileName(volFile)
+            Ts = self.inputVolume.get().getSamplingRate()
+            vol.setSamplingRate(Ts)
+            self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(volFile,Ts))
+            self._defineOutputs(outputVol=vol)
+        else:
+            self._defineOutputs(outputVol=self.inputVolume.get())
+        self._defineSourceRelation(self.inputVolume, self.outputVol)
+
+        cleanPath(self._getPath('resizedVol.mrc'))
+        cleanPath(self._getPath('mask.mrc'))
+        cleanPath(self._getPath('filteredVol.mrc'))
+
+# --------------------------- INFO functions -------------------------------
+    def _summary(self):
+        summary = []
+
+        if hasattr(self, 'outputHand'):
+            summary.append('Hand value is: %f' %self.outputHand.get())
+            summary.append('Hand values close to 1 mean the structure is predicted to be left handed')
+            summary.append('Hand values close to 0 mean the structure is predicted to be right handed')
+            if self.outputHand.get() > self.thresholdHand.get():
+                summary.append('Volume was flipped as it was deemed to be left handed')
+            else:
+                summary.append('Volume was not flipped as it was deemed to be right handed')
+
+        return summary
+
+    def _methods(self):
+        methods = []
+        return methods
+
+    def _validate(self):
+        errors = []
+
+        if self.thresholdAlpha.get() > 1.0 or self.thresholdAlpha.get() < 0.0:
+           errors.append("Alpha threshold must be between 0.0 and 1.0")
+        if self.thresholdHand.get() > 1.0 or self.thresholdHand.get() < 0.0:
+           errors.append("Hand threshold must be between 0.0 and 1.0")
+
+        return errors
diff --git a/xmipp3/protocols/protocol_denoise_particles.py b/xmipp3/protocols/protocol_denoise_particles.py
index 3745bd120..11cafb0a7 100644
--- a/xmipp3/protocols/protocol_denoise_particles.py
+++ b/xmipp3/protocols/protocol_denoise_particles.py
@@ -29,6 +29,7 @@
 from pyworkflow.protocol.params import IntParam, PointerParam, LEVEL_ADVANCED
 from pwem.protocols import ProtProcessParticles
 from pwem.objects import SetOfVolumes
+from pwem.objects import SetOfAverages
 
 from xmipp3.convert import (writeSetOfParticles, writeSetOfClasses2D,
                             xmippToLocation)
@@ -89,10 +90,10 @@ def denoiseImages(self, inputId, inputClassesId):
         imagesMd = self._getPath('images.xmd')
         writeSetOfParticles(self.inputParticles.get(), imagesMd)
         classesMd = self._getPath('classes.xmd')
-        if isinstance(self.inputClasses.get(), SetOfVolumes):
+        if isinstance(self.inputClasses.get(), SetOfAverages):
             writeSetOfClasses2D(self.inputClasses.get(), classesMd)
         else:
-            writeSetOfParticles(self.inputParticles.get(), classesMd)
+            writeSetOfClasses2D(self.inputClasses.get(), classesMd)
 
         fnRoot = self._getExtraPath('pca')
         fnRootDenoised = self._getExtraPath('imagesDenoised')
diff --git a/xmipp3/protocols/protocol_local_ctf.py b/xmipp3/protocols/protocol_local_ctf.py
index fd0874b9f..1f72ea8cf 100644
--- a/xmipp3/protocols/protocol_local_ctf.py
+++ b/xmipp3/protocols/protocol_local_ctf.py
@@ -26,11 +26,8 @@
 # *
 # **************************************************************************
 
-from math import floor
-import os
-
 from pyworkflow import VERSION_2_0
-from pyworkflow.protocol.params import PointerParam, FloatParam
+from pyworkflow.protocol.params import PointerParam, FloatParam, BooleanParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 from pyworkflow.utils.path import cleanPattern
 from pwem.protocols import ProtAnalysis3D
@@ -51,7 +48,7 @@ class XmippProtLocalCTF(ProtAnalysis3D):
     def __init__(self, **args):
         ProtAnalysis3D.__init__(self, **args)
     
-    #--------------------------- DEFINE param functions --------------------------------------------
+    # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputSet', PointerParam, label="Input images",
@@ -63,20 +60,25 @@ def _defineParams(self, form):
         form.addParam('maxDefocusChange', FloatParam, label="Maximum defocus change (A)", default=10000,
                       expertLevel=LEVEL_ADVANCED)
         form.addParam('maxGrayScaleChange', FloatParam, label="Maximum gray scale change", default=1,
-                      expertLevel=LEVEL_ADVANCED,
-                      help="The reprojection is modified as a*P+b, a is restricted to the interval [1-maxGrayScale,1+maxGrayScale]")
+                      expertLevel=LEVEL_ADVANCED, help="The reprojection is modified as a*P+b, a is restricted to the "
+                                                       "interval [1-maxGrayScale,1+maxGrayScale]")
         form.addParam('maxGrayShiftChange', FloatParam, label="Maximum gray shift change", default=1,
+                      expertLevel=LEVEL_ADVANCED, help="The reprojection is modified as a*P+b, b is restricted to the "
+                                                       "interval [-maxGrayShift,maxGrayShift]")
+        form.addParam('sameDefocus', BooleanParam, label="Force defocusV to be equal than defocusU", default=True,
                       expertLevel=LEVEL_ADVANCED,
-                      help = "The reprojection is modified as a*P+b, b is restricted to the interval [-maxGrayShift,maxGrayShift]")
+                      help="As the CTF usually suffers from astigmatism (it is not spherical but ellipsoidal), the "
+                           "defocus vary if computed in X or Y direction, being defocus U value the defocus in X "
+                           "direction and defocus V value the defocus in Y direction.")
         form.addParallelSection(threads=0, mpi=8)
     
-    #--------------------------- INSERT steps functions --------------------------------------------
+    # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
         self._insertFunctionStep("convertStep")
         self._insertFunctionStep("refineDefocus")
         self._insertFunctionStep("createOutputStep")
 
-    #--------------------------- STEPS functions ---------------------------------------------------
+    # --------------------------- STEPS functions ---------------------------------------------------
     def convertStep(self):
         """convert input to proper format and dimensions if necessary"""
         imgSet = self.inputSet.get()
@@ -87,31 +89,35 @@ def convertStep(self):
         xDimVol=self.inputVolume.get().getDim()[0]
         xDimImg = imgSet.getDim()[0]
         if xDimVol!=xDimImg:
-            self.runJob("xmipp_image_resize","-i %s --dim %d"%(fnVol,xDimImg),numberOfMpi=1)
+            self.runJob("xmipp_image_resize", "-i %s --dim %d"%(fnVol,xDimImg),numberOfMpi=1)
 
     def refineDefocus(self):
-        """compute local defocus using Xmipp (xmipp_angular_continuous_assign2) and add to metadata columns related to defocus"""
+        """compute local defocus using Xmipp (xmipp_angular_continuous_assign2) and add to metadata columns related to
+        defocus"""
         fnVol = self._getExtraPath("volume.vol")
         fnIn = self._getExtraPath('input_imgs.xmd')
         fnOut = self._getExtraPath('output_imgs.xmd')
-        anglesOutFn=self._getExtraPath("anglesCont.stk")
-        Ts=self.inputSet.get().getSamplingRate()
-        args="-i %s -o %s --ref %s --optimizeDefocus --max_defocus_change %d --sampling %f "\
-             "--optimizeGray --max_gray_scale %f --max_gray_shift %f"%\
-                    (fnIn,anglesOutFn,fnVol,self.maxDefocusChange.get(),Ts,self.maxGrayScaleChange.get(),
-                     self.maxGrayShiftChange.get())
+        anglesOutFn = self._getExtraPath("anglesCont.stk")
+        Ts = self.inputSet.get().getSamplingRate()
+        args = "-i %s -o %s --ref %s --optimizeDefocus --max_defocus_change %d --sampling %f --optimizeGray " \
+               "--max_gray_scale %f --max_gray_shift %f " % \
+               (fnIn, anglesOutFn, fnVol, self.maxDefocusChange.get(), Ts, self.maxGrayScaleChange.get(),
+                self.maxGrayShiftChange.get())
         if self.inputSet.get().isPhaseFlipped():
             args += " --phaseFlipped"
+        if self.sameDefocus.get():
+            args += " --sameDefocus"
         self.runJob("xmipp_angular_continuous_assign2", args)
 
         fnCont = self._getExtraPath('anglesCont.xmd')
-        self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "itemId ctfDefocusU ctfDefocusV ctfDefocusChange ctfDefocusAngle"'%
-                    fnCont,numberOfMpi=1)
+        self.runJob("xmipp_metadata_utilities", '-i %s --operate keep_column "itemId ctfDefocusU ctfDefocusV '
+                                                'ctfDefocusChange ctfDefocusAngle"'%
+                    fnCont, numberOfMpi=1)
         self.runJob("xmipp_metadata_utilities",
                     '-i %s -o %s --operate drop_column "ctfDefocusU ctfDefocusV ctfDefocusChange ctfDefocusAngle"' %
-                    (fnIn,fnOut),numberOfMpi=1)
+                    (fnIn, fnOut), numberOfMpi=1)
         self.runJob("xmipp_metadata_utilities",
-                    "-i %s --set join %s itemId itemId" % (fnOut, fnCont),numberOfMpi=1)
+                    "-i %s --set join %s itemId itemId" % (fnOut, fnCont), numberOfMpi=1)
 
         cleanPattern(self._getExtraPath("anglesCont.*"))
 
@@ -120,11 +126,12 @@ def createOutputStep(self):
         outputSet = self._createSetOfParticles()
         imgSet = self.inputSet.get()
         outputSet.copyInfo(imgSet)
-        readSetOfParticles(self._getExtraPath('output_imgs.xmd'), outputSet,  extraLabels=[emlib.MDL_CTF_DEFOCUS_CHANGE])
+        readSetOfParticles(self._getExtraPath('output_imgs.xmd'), outputSet,
+                           extraLabels=[emlib.MDL_CTF_DEFOCUS_CHANGE])
         self._defineOutputs(outputParticles=outputSet)
         self._defineSourceRelation(self.inputSet, outputSet)
 
-    #--------------------------- INFO functions --------------------------------------------
+    # --------------------------- INFO functions --------------------------------------------
     def _summary(self):
         summary = []
         summary.append("Refined defocus of %i particles" % self.inputSet.get().getSize())
@@ -135,6 +142,7 @@ def _summary(self):
     def _methods(self):
         methods = []
         if hasattr(self, 'outputParticles'):
-            methods.append("We refined the defocus %i input particles %s regarding to volume %s, allowing a maximun defocus of %s"\
-                           %(self.inputSet.get().getSize(), self.getObjectTag('inputSet'), self.getObjectTag('inputVolume'), self.maxDefocusChange.get()) )
+            methods.append("We refined the defocus %i input particles %s regarding to volume %s, allowing a maximun "
+                           "defocus of %s" % (self.inputSet.get().getSize(), self.getObjectTag('inputSet'),
+                                              self.getObjectTag('inputVolume'), self.maxDefocusChange.get()))
         return methods
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_crop_resize.py b/xmipp3/protocols/protocol_preprocess/protocol_crop_resize.py
index 93d66afcb..d244469bd 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_crop_resize.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_crop_resize.py
@@ -70,6 +70,7 @@ def _defineProcessParams(cls, protocol, form):
                       'between interpolation and Fourier cropping.')
         form.addParam('resizeSamplingRate', FloatParam, default=1.0,
                       condition='doResize and resizeOption==%d' % cls.RESIZE_SAMPLINGRATE,
+                      allowsPointers=True,
                       label='Resize sampling rate (Å/px)',
                       help='Set the new output sampling rate.')
         form.addParam('doFourier', BooleanParam, default=False,
@@ -78,14 +79,17 @@ def _defineProcessParams(cls, protocol, form):
                       help='If you set to *True*, the final dimensions must be lower than the original ones.')
         form.addParam('resizeDim', IntParam, default=0,
                       condition='doResize and resizeOption==%d' % cls.RESIZE_DIMENSIONS,
+                      allowsPointers=True,
                       label='New image size (px)',
                       help='Size in pixels of the particle images <x> <y=x> <z=x>.')
         form.addParam('resizeFactor', FloatParam, default=0.5,
                       condition='doResize and resizeOption==%d' % cls.RESIZE_FACTOR,
+                      allowsPointers=True,
                       label='Resize factor',
                       help='New size is the old one x resize factor.')
         form.addParam('resizeLevel', IntParam, default=0,
                       condition='doResize and resizeOption==%d' % cls.RESIZE_PYRAMID,
+                      allowsPointers=True,
                       label='Pyramid level',
                       help='Use positive value to expand and negative to reduce.')
         form.addParam('hugeFile', BooleanParam, default=False, expertLevel=const.LEVEL_ADVANCED,
@@ -108,12 +112,14 @@ def _defineProcessParams(cls, protocol, form):
                       '_crop_: choose how many pixels to crop from each border.\n')
         form.addParam('cropSize', IntParam, default=0,
                       condition='doWindow and windowOperation == %d' % cls.WINDOW_OP_CROP,
+                      allowsPointers=True,
                       label='Crop size (px)',
                       help='Amount of pixels cropped from each border.\n'
                            'e.g: if you set 10 pixels, the dimensions of the\n'
                            'object (SetOfParticles, Volume or SetOfVolumes) will be\n'
                            'reduced in 20 pixels (2 borders * 10 pixels)')
         form.addParam('windowSize', IntParam, default=0,
+                      allowsPointers=True,
                       condition='doWindow and windowOperation == %d' % cls.WINDOW_OP_WINDOW,
                       label='Window size (px)',
                       help='Size in pixels of the output object. It will be '
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_process.py b/xmipp3/protocols/protocol_preprocess/protocol_process.py
index b0a32033d..5419a5e5a 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_process.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_process.py
@@ -105,7 +105,7 @@ def createOutputStep(self):
     def _defineFilenames(self):
         self.inputFn = self._getTmpPath('input_particles.xmd')
         self.outputMd = self._getExtraPath('output_images.xmd')
-        self.outputStk = self._getExtraPath('output_images.stk')
+        self.outputStk = self._getExtraPath('output_images.mrcs')
 
 
 class XmippProcessVolumes(ProtPreprocessVolumes):
diff --git a/xmipp3/protocols/protocol_resolution_deepres.py b/xmipp3/protocols/protocol_resolution_deepres.py
index 9f4e0a5dc..0898a6554 100644
--- a/xmipp3/protocols/protocol_resolution_deepres.py
+++ b/xmipp3/protocols/protocol_resolution_deepres.py
@@ -52,6 +52,7 @@ def updateEnviron(gpuNum):
 OPERATE_VOL = 'operateVolume.vol'
 #CHIMERA_RESOLUTION_VOL = 'deepRes_resolution.vol'
 OUTPUT_RESOLUTION_FILE = 'resolutionMap'
+OUTPUT_ORIGINAL_SIZE =  'deepRes_resolution_originalSize.vol'
 OUTPUT_RESOLUTION_FILE_CHIMERA = 'chimera_resolution.vol'
 METADATA_MASK_FILE = 'metadataresolutions'
 FN_METADATA_HISTOGRAM = 'mdhist'
@@ -116,6 +117,7 @@ def _createFilenameTemplates(self):
                  OUTPUT_RESOLUTION_FILE_CHIMERA: self._getExtraPath('chimera_resolution.vol'),
 #                 OUTPUT_RESOLUTION_FILE_CHIMERA: self._getExtraPath(CHIMERA_RESOLUTION_VOL),                                 
                  OUTPUT_RESOLUTION_FILE: self._getExtraPath('deepRes_resolution.vol'),
+                 OUTPUT_ORIGINAL_SIZE: self._getExtraPath('deepRes_resolution_originalSize.vol'),
                  FN_METADATA_HISTOGRAM: self._getExtraPath('hist.xmd')
                  }
         self._updateFilenamesDict(myDict)
@@ -128,7 +130,7 @@ def _insertAllSteps(self):
         self._createFilenameTemplates() 
         self._insertFunctionStep('convertInputStep')
         self._insertFunctionStep('transformStep')          
-        self._insertFunctionStep('resizeStep')                        
+        self._insertFunctionStep('resizeInputStep')                        
         self._insertFunctionStep('resolutionStep')
         self._insertFunctionStep('createOutputStep')
         self._insertFunctionStep("createHistrogram")
@@ -166,7 +168,7 @@ def transformStep(self):
                   
             
             
-    def resizeStep(self):
+    def resizeInputStep(self):
 
         if self.range == self.LOW_RESOL:
             sampling_new = 1.0
@@ -281,10 +283,32 @@ def createOutputStep(self):
             sampling_new = 1.0
         else:
             sampling_new = 0.5
+            
+        #convert size of output volume
+        
+        samplingFactor = sampling_new/float(self.inputVolume.get().getSamplingRate())
+        fourierValue = sampling_new/(2*float(self.inputVolume.get().getSamplingRate()))
+             
+        if sampling_new > self.inputVolume.get().getSamplingRate():
+            paramsResizeVol = ' -i %s' % self._getFileName(OUTPUT_RESOLUTION_FILE)
+            paramsResizeVol += ' -o %s' % self._getFileName(OUTPUT_ORIGINAL_SIZE)        
+            paramsResizeVol += ' --factor %s' % samplingFactor
+            self.runJob('xmipp_image_resize', paramsResizeVol )             
+        else:   
+            paramsFilterVol = ' -i %s' % self._getFileName(OUTPUT_RESOLUTION_FILE)
+            paramsFilterVol += ' -o %s' % self._getFileName(OUTPUT_ORIGINAL_SIZE)
+            paramsFilterVol += ' --fourier low_pass %s' % fourierValue            
+            paramsResizeVol = ' -i %s' % self._getFileName(OUTPUT_ORIGINAL_SIZE)  
+            paramsResizeVol += ' -o %s' % self._getFileName(OUTPUT_ORIGINAL_SIZE)           
+            paramsResizeVol += ' --factor %s' % samplingFactor            
+            self.runJob('xmipp_transform_filter', paramsFilterVol )
+            self.runJob('xmipp_image_resize', paramsResizeVol )
+            
+            
         volume=Volume()
-        volume.setFileName(self._getFileName(OUTPUT_RESOLUTION_FILE))
+        volume.setFileName(self._getFileName(OUTPUT_ORIGINAL_SIZE))
 
-        volume.setSamplingRate(sampling_new)
+        volume.setSamplingRate(self.inputVolume.get().getSamplingRate())
         self._defineOutputs(resolution_Volume=volume)
         self._defineTransformRelation(self.inputVolume, volume)
             
diff --git a/xmipp3/protocols/protocol_rotate_volume.py b/xmipp3/protocols/protocol_rotate_volume.py
new file mode 100644
index 000000000..bca41897b
--- /dev/null
+++ b/xmipp3/protocols/protocol_rotate_volume.py
@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+# **************************************************************************
+# *
+# * Authors:  Estrella Fernandez Gimenez (me.fernandez@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.protocol.params import PointerParam, EnumParam, FloatParam
+from pwem.convert.headers import setMRCSamplingRate
+from pwem.objects.data import Volume
+from pwem.protocols import EMProtocol
+
+
+class XmippProtRotateVolume(EMProtocol):
+    """ Rotate a volume around x,y,z """
+
+    _label = 'rotate volume'
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('vol', PointerParam, pointerClass='Volume', label="Volume ", help='Specify a volume.')
+        form.addParam('rotType', EnumParam, choices=['Align with Z', 'rotate'], display=EnumParam.DISPLAY_HLIST,
+                      default=1, label='Rotation mode: ', help='Align (x,y,z) with Z axis')
+        form.addParam('dirParam', EnumParam, choices=['X', 'Y', 'Z'], default=1, display=EnumParam.DISPLAY_HLIST,
+                      label='Axis: ', help='Align (x,y,z) with Z axis')
+        form.addParam('deg', FloatParam, label='Degrees: ', default=90, condition='rotType == 1',
+                      help='degrees of rotation in selected axis')
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('rotateStep')
+        self._insertFunctionStep('createOutputStep')
+
+    # --------------------------- STEPS functions --------------------------------------------
+    def rotateStep(self):
+        vol = self.vol.get()
+        fnVol = vol.getFileName()
+        if fnVol.endswith('.mrc'):
+            fnVol += ':mrc'
+        args = '-i %s -o %s --rotate_volume' % (fnVol, self._getExtraPath('rotated_vol.mrc'))
+        rotType = self.rotType.get()
+        if rotType == 0:
+            args += ' alignZ'
+        if rotType == 1:
+            args += ' axis %d' % self.deg.get()
+        dirParam = self.dirParam.get()
+        if dirParam == 0:
+            args += ' 1 0 0'
+        if dirParam == 1:
+            args += ' 0 1 0'
+        if dirParam == 2:
+            args += ' 0 0 1'
+        program = "xmipp_transform_geometry"
+        self.runJob(program, args)
+
+    def createOutputStep(self):
+        outputVol = Volume()
+        inputVol = self.vol.get()
+        outputVol.copyInfo(inputVol)
+        outputVol.setLocation(self._getExtraPath('rotated_vol.mrc'))
+        setMRCSamplingRate(outputVol.getFileName(), inputVol.getSamplingRate())
+        self._defineOutputs(outputVolume=outputVol)
+        self._defineSourceRelation(inputVol, outputVol)
+
+    # --------------------------- INFO functions --------------------------------------------
+    def _summary(self):
+        summary = 'Volume'
+        rotType = self.rotType.get()
+        if rotType == 0:
+            summary += ' aligned with'
+        if rotType == 1:
+            summary += ' rotated %d degrees around' % self.deg.get()
+        dirParam = self.dirParam.get()
+        if dirParam == 0:
+            summary += ' X axis'
+        if dirParam == 1:
+            summary += ' Y axis'
+        if dirParam == 2:
+            summary += ' Z axis'
+        return summary
diff --git a/xmipp3/protocols/protocol_screen_deeplearning.py b/xmipp3/protocols/protocol_screen_deeplearning.py
index 4ba311dff..4ebec7f1a 100644
--- a/xmipp3/protocols/protocol_screen_deeplearning.py
+++ b/xmipp3/protocols/protocol_screen_deeplearning.py
@@ -38,13 +38,14 @@
 
 N_MAX_NEG_SETS = 5
 
+
 class XmippProtScreenDeepLearning(ProtProcessParticles, XmippProtocol):
     """ Protocol for screening particles using deep learning. """
     _label = 'screen deep learning'
     _lastUpdateVersion = VERSION_2_0
     _conda_env = 'xmipp_DLTK_v0.3'
 
-    #--------------------------- DEFINE param functions --------------------------------------------
+    # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
         # GPU settings
         form.addHidden(params.USE_GPU, params.BooleanParam, default=True,
@@ -70,9 +71,9 @@ def _defineParams(self, form):
                       help='Select a previous run to continue from.')
         form.addParam('keepTraining', params.BooleanParam,
                       label='Continue training on previously trainedModel?',
-                      default=True,condition='doContinue',
+                      default=True, condition='doContinue',
                       help='If you set to *Yes*, you should provide training set')
-            
+
         form.addParam('inTrueSetOfParticles', params.PointerParam,
                       label="True particles", pointerClass='SetOfParticles',
                       allowsNull=True, condition="not doContinue or keepTraining",
@@ -91,14 +92,14 @@ def _defineParams(self, form):
                                     'and (not doContinue or keepTraining)' % num,
                           pointerClass='SetOfParticles', allowsNull=True,
                           help='Select the set of negative particles for training.')
-                          
-            form.addParam('inNegWeight_%d'%num, params.IntParam,
+
+            form.addParam('inNegWeight_%d' % num, params.IntParam,
                           label="Weight of negative train particles %d" % num,
                           expertLevel=params.LEVEL_ADVANCED,
                           default='1', allowsNull=True,
                           condition='(numberOfNegativeSets<=0 or numberOfNegativeSets >=%d) and '
-                                    '(not doContinue or keepTraining)'%num,
-                          help='Select the weigth for the negative set of particles. '
+                                    '(not doContinue or keepTraining)' % num,
+                          help='Select the weigh for the negative set of particles. '
                                'The weight value indicates the number of times '
                                'each image may be included at most per epoch. '
                                'Positive particles are weighted with 1. '
@@ -110,10 +111,10 @@ def _defineParams(self, form):
                       label="Set of putative particles to score",
                       pointerClass='SetOfParticles',
                       help='Select the set of putative particles to classify as good (score close '
-                            'to 1.0) or bad (score close to 0.0).')
+                           'to 1.0) or bad (score close to 0.0).')
 
         form.addSection(label='Training')
-        
+
         form.addParam('nEpochs', params.FloatParam,
                       label="Number of epochs", default=5.0,
                       condition="not doContinue or keepTraining",
@@ -122,15 +123,15 @@ def _defineParams(self, form):
                       label="Learning rate", default=1e-4,
                       condition="not doContinue or keepTraining",
                       help='Learning rate for neural network training')
-        form.addParam('auto_stopping',params.BooleanParam,
-                      label='Auto stop training when convergency is detected?',
+        form.addParam('auto_stopping', params.BooleanParam,
+                      label='Auto stop training when convergence is detected?',
                       default=True, condition="not doContinue or keepTraining",
                       help='If you set to *Yes*, the program will automatically '
                            'stop training if there is no improvement for '
                            'consecutive 2 epochs, learning rate will be '
                            'decreased by a factor 10. '
                            'If learningRate_t < 0.01*learningrate_0 training will stop. '
-                           'Warning: Sometimes convergency seems to be reached, '
+                           'Warning: Sometimes convergence seems to be reached, '
                            'but after time, improvement can still happen. '
                            'Not recommended for very small data sets (<100 true particles)')
         form.addParam('l2RegStrength', params.FloatParam,
@@ -144,18 +145,18 @@ def _defineParams(self, form):
                       label="Number of models for ensemble",
                       default=2, expertLevel=params.LEVEL_ADVANCED,
                       condition="not doContinue",
-                      help='Number of models to fit in order to build an ensamble. '
+                      help='Number of models to fit in order to build an ensemble. '
                            'Tipical values are 1 to 5. The more the better '
                            'until a point where no gain is obtained. '
                            'Each model increases running time linearly')
-                           
+
         form.addParam('doTesting', params.BooleanParam, default=False,
                       label='Perform testing after training?', expertLevel=params.LEVEL_ADVANCED,
                       help='If you set to *Yes*, you should select a testing '
                            'positive set and a testing negative set')
         form.addParam('testPosSetOfParticles', params.PointerParam,
                       label="Set of positive test particles", expertLevel=params.LEVEL_ADVANCED,
-                      pointerClass='SetOfParticles',condition='doTesting',
+                      pointerClass='SetOfParticles', condition='doTesting',
                       help='Select the set of ground true positive particles.')
         form.addParam('testNegSetOfParticles', params.PointerParam,
                       label="Set of negative test particles", expertLevel=params.LEVEL_ADVANCED,
@@ -165,10 +166,11 @@ def _defineParams(self, form):
     def _validate(self):
         return self.validateDLtoolkit()
 
-    #--------------------------- INSERT steps functions --------------------------------------------
+    # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
         """            
         """
+
         def _getFname2WeightDict(fnameToSetAndWeight):
             """
             arg: fnameToSetAndWeight= { fname: [(SetOfParticles, weight:int)]}
@@ -178,8 +180,8 @@ def _getFname2WeightDict(fnameToSetAndWeight):
             if fnameToSetAndWeight is None:
                 return None
             dictONameToWeight = {fname: fnameToSetAndWeight[fname][-1]
-                                   for fname in fnameToSetAndWeight
-                                     if not fnameToSetAndWeight[fname][0] is None}
+                                 for fname in fnameToSetAndWeight
+                                 if not fnameToSetAndWeight[fname][0] is None}
 
             if len(dictONameToWeight) == 0:
                 return None
@@ -187,68 +189,68 @@ def _getFname2WeightDict(fnameToSetAndWeight):
                 return dictONameToWeight
 
         posSetTrainDict = {self._getExtraPath("inputTrueParticlesSet.xmd"): 1}
-        
+
         negSetTrainDict = {}
         for num in range(1, N_MAX_NEG_SETS):
             if self.numberOfNegativeSets <= 0 or self.numberOfNegativeSets >= num:
-                negSetTrainDict[self._getExtraPath("negativeSet_%d.xmd"%num)] = \
-                                          self.__dict__["inNegWeight_%d"%num].get()
-                    
+                negSetTrainDict[self._getExtraPath("negativeSet_%d.xmd" % num)] = \
+                    self.__dict__["inNegWeight_%d" % num].get()
+
         setPredictDict = {self._getExtraPath("predictSetOfParticles.xmd"): 1}
 
         if self.doTesting.get() and self.testPosSetOfParticles.get() and self.testNegSetOfParticles.get():
-          setTestPosDict = {self._getExtraPath("testTrueParticlesSet.xmd"): 1}
-          setTestNegDict = {self._getExtraPath("testFalseParticlesSet.xmd"): 1}
+            setTestPosDict = {self._getExtraPath("testTrueParticlesSet.xmd"): 1}
+            setTestNegDict = {self._getExtraPath("testFalseParticlesSet.xmd"): 1}
         else:
-          setTestPosDict = None
-          setTestNegDict = None
-          
+            setTestPosDict = None
+            setTestNegDict = None
+
         self._insertFunctionStep('convertInputStep', posSetTrainDict,
                                  negSetTrainDict, setPredictDict,
                                  setTestPosDict, setTestNegDict)
         if not self.doContinue.get() or self.keepTraining.get():
-          self._insertFunctionStep('train', posSetTrainDict, negSetTrainDict)
-        self._insertFunctionStep('predict', setTestPosDict, setTestNegDict,setPredictDict)
+            self._insertFunctionStep('train', posSetTrainDict, negSetTrainDict)
+        self._insertFunctionStep('predict', setTestPosDict, setTestNegDict, setPredictDict)
         self._insertFunctionStep('createOutputStep')
 
-    #--------------------------- STEPS functions -------------------------------
+    # --------------------------- STEPS functions -------------------------------
     def convertInputStep(self, *dataDicts):
         def __getSetOfParticlesFromFname(fname):
-          if fname== self._getExtraPath("inputTrueParticlesSet.xmd"):
-            return self.inTrueSetOfParticles.get()
-          elif fname== self._getExtraPath("predictSetOfParticles.xmd"):
-            return self.predictSetOfParticles.get()
-          elif fname== self._getExtraPath("testTrueParticlesSet.xmd"):
-            return self.testPosSetOfParticles.get()
-          elif fname== self._getExtraPath("testFalseParticlesSet.xmd"):
-            return self.testNegSetOfParticles.get()
-          else:
-            matchOjb= re.match( self._getExtraPath(r"negativeSet_(\d+).xmd"), fname)
-            if matchOjb:
-              num= matchOjb.group(1)
-              return self.__dict__["negativeSet_%s"%num].get()
+            if fname == self._getExtraPath("inputTrueParticlesSet.xmd"):
+                return self.inTrueSetOfParticles.get()
+            elif fname == self._getExtraPath("predictSetOfParticles.xmd"):
+                return self.predictSetOfParticles.get()
+            elif fname == self._getExtraPath("testTrueParticlesSet.xmd"):
+                return self.testPosSetOfParticles.get()
+            elif fname == self._getExtraPath("testFalseParticlesSet.xmd"):
+                return self.testNegSetOfParticles.get()
             else:
-              raise ValueError("Error, unexpected fname")
-                     
+                matchOjb = re.match(self._getExtraPath(r"negativeSet_(\d+).xmd"), fname)
+                if matchOjb:
+                    num = matchOjb.group(1)
+                    return self.__dict__["negativeSet_%s" % num].get()
+                else:
+                    raise ValueError("Error, unexpected fname")
+
         if ((not self.doContinue.get() or self.keepTraining.get())
                 and self.nEpochs.get() > 0):
             assert not self.inTrueSetOfParticles.get() is None, \
-                    "Positive particles must be provided for training if nEpochs!=0"
+                "Positive particles must be provided for training if nEpochs!=0"
 
         for dataDict in dataDicts:
             if not dataDict is None:
                 for fnameParticles in sorted(dataDict):
-                    setOfParticles= __getSetOfParticlesFromFname(fnameParticles)
+                    setOfParticles = __getSetOfParticlesFromFname(fnameParticles)
                     writeSetOfParticles(setOfParticles, fnameParticles)
-                    
+
     def __dataDict_toStrs(self, dataDict):
-        fnamesStr=[]
-        weightsStr=[]
+        fnamesStr = []
+        weightsStr = []
         for fname in dataDict:
-          fnamesStr.append(fname)
-          weightsStr.append(str(dataDict[fname]) )
+            fnamesStr.append(fname)
+            weightsStr.append(str(dataDict[fname]))
         return ":".join(fnamesStr), ":".join(weightsStr)
-        
+
     def train(self, posTrainDict, negTrainDict):
         """
             posTrainDict, negTrainDict: { fnameToMetadata: weight (int) }
@@ -260,7 +262,7 @@ def train(self, posTrainDict, negTrainDict):
             copyTree(prevRunPath, netDataPath)
             if not self.keepTraining.get():
                 nEpochs = 0
-                
+
         if self.usesGpu():
             numberOfThreads = None
             gpuToUse = self.getGpuList()[0]
@@ -268,20 +270,19 @@ def train(self, posTrainDict, negTrainDict):
             numberOfThreads = self.numberOfThreads.get()
             gpuToUse = None
 
-          
-        fnamesPos, weightsPos= self.__dataDict_toStrs(posTrainDict)
-        fnamesNeg, weightsNeg= self.__dataDict_toStrs(negTrainDict)
-        args= " -n %s --mode train -p %s -f %s --trueW %s --falseW %s"%(netDataPath, 
-                      fnamesPos, fnamesNeg, weightsPos, weightsNeg)
-        args+= " -e %s -l %s -r %s -m %s "%(nEpochs, self.learningRate.get(), self.l2RegStrength.get(),
-                                          self.nModels.get())
+        fnamesPos, weightsPos = self.__dataDict_toStrs(posTrainDict)
+        fnamesNeg, weightsNeg = self.__dataDict_toStrs(negTrainDict)
+        args = " -n %s --mode train -p %s -f %s --trueW %s --falseW %s" % (netDataPath,
+                                                                           fnamesPos, fnamesNeg, weightsPos, weightsNeg)
+        args += " -e %s -l %s -r %s -m %s " % (nEpochs, self.learningRate.get(), self.l2RegStrength.get(),
+                                               self.nModels.get())
         if not self.auto_stopping.get():
-          args+=" -s"
-          
+            args += " -s"
+
         if not gpuToUse is None:
-          args+= " -g %s"%(gpuToUse)
+            args += " -g %s" % (gpuToUse)
         if not numberOfThreads is None:
-          args+= " -t %s"%(numberOfThreads)
+            args += " -t %s" % (numberOfThreads)
         self.runJob('xmipp_deep_consensus', args, numberOfMpi=1,
                     env=self.getCondaEnv())
 
@@ -302,19 +303,19 @@ def predict(self, posTestDict, negTestDict, predictDict):
             gpuToUse = None
 
         outParticlesPath = self._getPath("particles.xmd")
-        fnamesPred, weightsPred= self.__dataDict_toStrs(predictDict)
+        fnamesPred, weightsPred = self.__dataDict_toStrs(predictDict)
 
-        args= " -n %s --mode score -i %s -o %s "%(netDataPath, fnamesPred, outParticlesPath)
+        args = " -n %s --mode score -i %s -o %s " % (netDataPath, fnamesPred, outParticlesPath)
 
         if posTestDict and negTestDict:
-          fnamesPosTest, weightsPosTest= self.__dataDict_toStrs(posTestDict)
-          fnamesNegTest, weightsNegTest= self.__dataDict_toStrs(negTestDict)
-          args+= " --testingTrue %s --testingFalse %s "%(fnamesPosTest, fnamesNegTest)
+            fnamesPosTest, weightsPosTest = self.__dataDict_toStrs(posTestDict)
+            fnamesNegTest, weightsNegTest = self.__dataDict_toStrs(negTestDict)
+            args += " --testingTrue %s --testingFalse %s " % (fnamesPosTest, fnamesNegTest)
 
         if not gpuToUse is None:
-          args+= " -g %s"%(gpuToUse)
+            args += " -g %s" % (gpuToUse)
         if not numberOfThreads is None:
-          args+= " -t %s"%(numberOfThreads)
+            args += " -t %s" % (numberOfThreads)
         self.runJob('xmipp_deep_consensus', args, numberOfMpi=1,
                     env=self.getCondaEnv())
 
@@ -328,8 +329,7 @@ def createOutputStep(self):
         self._defineOutputs(outputParticles=partSet)
         self._defineSourceRelation(imgSet, partSet)
 
-
-    #--------------------------- INFO functions --------------------------------------------
+    # --------------------------- INFO functions --------------------------------------------
     def _summary(self):
         summary = []
         return summary
@@ -337,11 +337,10 @@ def _summary(self):
     def _methods(self):
         pass
 
-    #--------------------------- UTILS functions --------------------------------------------
+    # --------------------------- UTILS functions --------------------------------------------
     def _updateParticle(self, item, row):
         setXmippAttributes(item, row, md.MDL_ZSCORE_DEEPLEARNING1)
         if row.getValue(md.MDL_ENABLED) <= 0:
             item._appendItem = False
         else:
             item._appendItem = True
-
diff --git a/xmipp3/protocols/protocol_screen_particles.py b/xmipp3/protocols/protocol_screen_particles.py
index 58cc8de09..9a86ca5a5 100644
--- a/xmipp3/protocols/protocol_screen_particles.py
+++ b/xmipp3/protocols/protocol_screen_particles.py
@@ -42,19 +42,16 @@
 from pwem.objects import SetOfParticles
 from pwem.protocols import ProtProcessParticles
 
-
 from pwem import emlib
 from xmipp3.convert import readSetOfParticles, writeSetOfParticles
 
 
 class XmippProtScreenParticles(ProtProcessParticles):
-    """ Attach different merit values to every particle in order to
-        prune the set.
-        zScore evaluates the similarity of a particles with an average
-        (lower zScore -> high similarity).
-        SSNR evaluates the signal/noise ration in the Fourier space.
-        Variance evaluates the varaince on the micrographs context where
-        the particle was picked.
+    """Protocol to attach different merit values to every particle metadata for subsequent pruning the set.
+There are different merit values to be calculated:
+    - zScore evaluates the similarity of a particles with an average (lower zScore -> higher similarity).
+    - SSNR evaluates the signal/noise ration in the Fourier space.
+    - Variance evaluates the varaince on the micrographs context where the particle was picked.
     """
 
     _label = 'screen particles'
@@ -136,7 +133,7 @@ def _getDefaultParallel(self):
         """This protocol doesn't have mpi version"""
         return (0, 0)
 
-    #--------------------------- INSERT steps functions ----------------------
+    # --------------------------- INSERT steps functions ----------------------
     def _insertAllSteps(self):
         self._initializeZscores()
         self.outputSize = 0
@@ -275,7 +272,7 @@ def _loadOutputSet(self, SetClass, baseName):
 
         return outputSet
 
-    #--------------------------- STEP functions -----------------------------
+    # --------------------------- STEP functions -----------------------------
     def sortImagesStep(self):
         args = "-i Particles@%s -o %s --addToInput " % (self.fnInputMd,
                                                         self.fnOutputMd)
@@ -312,7 +309,7 @@ def sortImagesStep(self):
                     valuesList = varList
                     self.mdLabels = [emlib.MDL_SCORE_BY_VAR]
                 else:  # not working pretty well
-                    valuesList = [var*(1-gini) for var, gini in zip(varList, giniList)]
+                    valuesList = [var * (1 - gini) for var, gini in zip(varList, giniList)]
                     self.mdLabels = [emlib.MDL_SCORE_BY_VAR, emlib.MDL_SCORE_BY_GINI]
 
                 self.varThreshold.set(histThresholding(valuesList))
@@ -360,7 +357,7 @@ def _summary(self):
 
         if self.autoParRejectionVar is not None:
             sumRejMet['Var'] = ("Variance rejection method: " +
-                               self.VAR_CHOICES[self.autoParRejectionVar.get()])
+                                self.VAR_CHOICES[self.autoParRejectionVar.get()])
 
         # If no output yet, just the form choices are shown plus a no-ready text
         if not hasattr(self, 'outputParticles'):
@@ -384,7 +381,7 @@ def _summary(self):
                 else:
                     summary.append(" - Variance threshold not calculed yet.")
         return summary
-    
+
     def _validate(self):
         validateMsgs = []
         if self.autoParRejectionVar != self.REJ_NONE:
@@ -395,18 +392,18 @@ def _validate(self):
                                     'scoreByVariance attribute. Use Xmipp to '
                                     'extract the particles.')
         return validateMsgs
-        
+
     def _citations(self):
         return ['Vargas2013b']
-    
+
     def _methods(self):
         methods = []
         if hasattr(self, 'outputParticles'):
             outParticles = (len(self.outputParticles) if self.outputParticles
                                                          is not None else None)
-            particlesRejected = (len(self.inputParticles.get())-outParticles
+            particlesRejected = (len(self.inputParticles.get()) - outParticles
                                  if outParticles is not None else None)
-            particlesRejectedText = (' ('+str(particlesRejected)+')' if
+            particlesRejectedText = (' (' + str(particlesRejected) + ')' if
                                      particlesRejected is not None else '')
             rejectionText = ['',  # REJ_NONE
                              ' and removing those not reaching %s%s'
@@ -438,8 +435,8 @@ def histThresholding(valuesList, nBins=256, portion=4, takeNegatives=False):
         valuesList = [x for x in valuesList if not x < 0]
 
     import numpy as np
-    while len(valuesList)*1.0/nBins < 5:
-        nBins = int(nBins/2)
+    while len(valuesList) * 1.0 / nBins < 5:
+        nBins = int(nBins / 2)
 
     print('Thresholding with %d bins for the histogram.' % nBins)
 
@@ -448,9 +445,10 @@ def histThresholding(valuesList, nBins=256, portion=4, takeNegatives=False):
     histRight = hist
     histRight[0:hist.argmax()] = 0
 
-    idx = (np.abs(histRight-hist.max()/portion)).argmin()
+    idx = (np.abs(histRight - hist.max() / portion)).argmin()
     return bin_edges[idx]
 
+
 def rejectByVariance(inputMdFn, outputMdFn, threshold, mode):
     """ Sets MDL_ENABLED to -1 to those items with a higher value
         than the threshold
@@ -462,7 +460,7 @@ def rejectByVariance(inputMdFn, outputMdFn, threshold, mode):
                 mdata.setValue(emlib.MDL_ENABLED, -1, objId)
         elif mode == XmippProtScreenParticles.REJ_VARGINI:
             if (mdata.getValue(emlib.MDL_SCORE_BY_VAR, objId) *
-                (1 - mdata.getValue(emlib.MDL_SCORE_BY_GINI, objId)) > threshold):
+                    (1 - mdata.getValue(emlib.MDL_SCORE_BY_GINI, objId)) > threshold):
                 mdata.setValue(emlib.MDL_ENABLED, -1, objId)
 
     mdata.write(outputMdFn)
diff --git a/xmipp3/protocols/protocol_subtract_projection.py b/xmipp3/protocols/protocol_subtract_projection.py
index ce183d8af..9eecea84a 100644
--- a/xmipp3/protocols/protocol_subtract_projection.py
+++ b/xmipp3/protocols/protocol_subtract_projection.py
@@ -29,15 +29,13 @@
 from os.path import basename
 from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
-import pwem.emlib.metadata as md
+from pwem import emlib
 from pwem.protocols import EMProtocol
-from xmipp3.convert import writeSetOfParticles
+from xmipp3.convert import writeSetOfParticles, readSetOfParticles
 
 
 class XmippProtSubtractProjection(EMProtocol):
-    """ This protocol computes the subtraction between particles and a initial volume, by computing its projections
-    with the same angles that input particles have. Then, each particle and the correspondent projection of the initial
-    volume are numerically adjusted and subtracted using a mask which denotes the region to keep. """
+    """ This protocol computes the subtraction between particles and a reference volume, by computing its projections with the same angles that input particles have. Then, each particle and the correspondent projection of the reference volume are numerically adjusted and subtracted using a mask which denotes the region to keep. """
 
     _label = 'subtract projection'
     INPUT_PARTICLES = "input_particles.xmd"
@@ -46,27 +44,22 @@ class XmippProtSubtractProjection(EMProtocol):
     def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('particles', PointerParam, pointerClass='SetOfParticles', label="Particles: ",
-                      help='Specify a SetOfParticles.')
-        form.addParam('vol', PointerParam, pointerClass='Volume', label="Initial volume ", help='Specify a volume.')
-        form.addParam('maskVol', PointerParam, pointerClass='VolumeMask', label='Volume mask', allowsNull=True,
-                      help='3D mask for the input volume. This mask is not mandatory but advisable.')
-        form.addParam('mask', PointerParam, pointerClass='VolumeMask', label="Mask for region to keep", allowsNull=True,
-                      help='Specify a 3D mask for the region of the input volume that you want to keep. If no mask is '
-                           'selected, the whole image will be subtracted')
-        form.addParam('resol', FloatParam, label="Filter at resolution: ", default=3, allowsNull=True,
-                      expertLevel=LEVEL_ADVANCED,
-                      help='Resolution (A) at which subtraction will be performed, filtering the volume projections.'
-                           'Value 0 implies no filtering.')
-        form.addParam('sigma', FloatParam, label="Decay of the filter (sigma): ", default=3, condition='resol',
-                      help='Decay of the filter (sigma parameter) to smooth the mask transition',
-                      expertLevel=LEVEL_ADVANCED)
-        form.addParam('iter', IntParam, label="Number of iterations: ", default=5, expertLevel=LEVEL_ADVANCED,
-                      help='Number of iterations for the adjustment process of the images before the subtraction itself'
-                           'several iterations are recommended to improve the adjustment.')
-        form.addParam('rfactor', FloatParam, label="Relaxation factor (lambda): ", default=1,
-                      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')
+                      help='Specify a SetOfParticles')
+        form.addParam('vol', PointerParam, pointerClass='Volume', label="Reference volume ", help='Specify a volume.')
+        form.addParam('mask', PointerParam, pointerClass='VolumeMask', label='Mask for region to keep', allowsNull=True,
+                      help='Specify a 3D mask for the region of the input volume that you want to keep. '
+                           'If no mask is given, the subtraction is performed in whole images.')
+        form.addParam('mwidth', FloatParam, label="Extra width in final mask: ", default=40, expertLevel=LEVEL_ADVANCED,
+                      help='Length (in A) to add for each side to the final result mask. -1 means no mask.')
+        form.addParam('resol', FloatParam, label="Maximum resolution: ", default=3,  expertLevel=LEVEL_ADVANCED,
+                      help='Maximum resolution (in A) of the data ')
+        form.addParam('limit_freq', BooleanParam, label="Limit frequency?: ", default=False, expertLevel=LEVEL_ADVANCED,
+                      help='Limit frequency in the adjustment process to the frequency correspondent to the resolution'
+                           ' indicated in "Maximum resolution" field above')
+        form.addParam('sigma', FloatParam, label="Decay of the filter (sigma): ", default=3, expertLevel=LEVEL_ADVANCED,
+                      help='Decay of the filter (sigma) to smooth the mask transition')
+        form.addParam('pad', IntParam, label="Fourier padding factor: ", default=2, expertLevel=LEVEL_ADVANCED,
+                      help='The volume is zero padded by this factor to produce projections')
 
     # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
@@ -83,52 +76,47 @@ def subtractionStep(self):
         fnVol = vol.getFileName()
         if fnVol.endswith('.mrc'):
             fnVol += ':mrc'
-        resol = self.resol.get()
-        iters = self.iter.get()
-        program = "xmipp_subtract_projection"
-        args = '-i %s --ref %s -o %s --iter %s --lambda %s' % (self._getExtraPath(self.INPUT_PARTICLES), fnVol,
-                                                               self._getExtraPath("output_particles"), iters,
-                                                               self.rfactor.get())
-        args += ' --saveProj %s' % self._getExtraPath('')
-        if resol:
-            fc = vol.getSamplingRate()/resol
-            args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())
-        if self.maskVol.get() is not None:
-            args += ' --maskVol %s' % self.maskVol.get().getFileName()
-        if self.mask.get() is not None:
-            args += ' --mask %s' % self.mask.get().getFileName()
-        else:
-            involdim = vol.getDim()
-            fnDescr = self._getExtraPath("mask.descr")
-            fhDescr = open(fnDescr, 'w')
-            fhDescr.write("%d %d %d 0" % (involdim[0], involdim[1], involdim[2]))
-            fhDescr.close()
-            mskKeep = self._getExtraPath("mask_keep.mrc")
-            args_mask_keep = "-i %s -o %s" % (fnDescr, mskKeep)
-            self.runJob("xmipp_phantom_create", args_mask_keep, numberOfMpi=1)
-            args_imageheader2 = "-i %s --sampling_rate %f" % (mskKeep, vol.getSamplingRate())
-            self.runJob("xmipp_image_header", args_imageheader2, numberOfMpi=1)
-            args += ' --mask %s --subAll' % mskKeep
-        self.runJob(program, args)
+        args = '-i %s --ref %s -o %s --sampling %f --max_resolution %f --fmask_width %f --padding %f ' \
+               '--sigma %d --limit_freq %d' % \
+               (self._getExtraPath(self.INPUT_PARTICLES), fnVol, self._getExtraPath("output_particles"),
+                vol.getSamplingRate(), self.resol.get(), self.mwidth.get(), self.pad.get(), self.sigma.get(),
+                int(self.limit_freq.get()))
+        mask = self.mask.get()
+        if mask is not None:
+            args += ' --mask %s' % mask.getFileName()
+        args += ' --save %s' % self._getExtraPath()
+        self.runJob("xmipp_subtract_projection", args)
 
     def createOutputStep(self):
-        self.ix = 0  # initiate counter for particle file name
         inputSet = self.particles.get()
         outputSet = self._createSetOfParticles()
         outputSet.copyInfo(inputSet)
-        outputSet.copyItems(inputSet, updateItemCallback=self._updateItem,
-                            itemDataIterator=md.iterRows(self._getExtraPath(self.INPUT_PARTICLES)))
+        readSetOfParticles(self._getExtraPath('input_particles.xmd'), outputSet, extraLabels=[emlib.MDL_SUBTRACTION_R2])
         self._defineOutputs(outputParticles=outputSet)
         self._defineSourceRelation(inputSet, outputSet)
 
     # --------------------------- INFO functions --------------------------------------------
+    def _validate(self):
+        errors = []
+        part = self.particles.get().getFirstItem()
+        vol = self.vol.get()
+        mask = self.mask.get()
+        if part.getDim()[0] != vol.getDim()[0]:
+            errors.append("Input particles and volume should have same X and Y dimensions")
+        if part.getSamplingRate() != vol.getSamplingRate():
+            errors.append("Input particles and volume should have same sampling rate")
+        if mask:
+            if vol.getSamplingRate() != mask.getSamplingRate():
+                errors.append("Input volume and mask should have same sampling rate")
+            if vol.getDim() != mask.getDim():
+                errors.append("Input volume and mask should have same dimensions")
+        if self.resol.get() == 0:
+            errors.append("Resolution (angstroms) should be bigger than 0")
+        return errors
+
     def _summary(self):
-        summary = ["Volume: %s" % self.vol.get().getFileName()]
-        summary.append("Set of particles: %s" % self.particles.get())
-        if self.mask:
-            summary.append("Mask: %s" % self.mask.get().getFileName())
-        if self.resol.get() != 0:
-            summary.append("Subtraction at resolution %f A" % self.resol.get())
+        summary = ["Volume: %s\nSet of particles: %s\nMask: %s" %
+                   (self.vol.get().getFileName(), self.particles.get(), self.mask.get().getFileName())]
         return summary
 
     def _methods(self):
@@ -136,24 +124,6 @@ def _methods(self):
         if not hasattr(self, 'outputParticles'):
             methods.append("Output particles not ready yet.")
         else:
-            methods.append("Volume projections from %s subtracted from particles" %
-                           basename(self.vol.get().getFileName()))
-            if self.mask:
-                methods.append("with mask %s" % basename(self.mask.get().getFileName()))
-            if self.resol.get() != 0:
-                methods.append(" at resolution %f A" % self.resol.get())
+            methods.append("Volume projections subtracted to particles keeping the region in %s"
+                           % basename(self.mask.get().getFileName()))
         return methods
-
-    def _validate(self):
-        errors = []
-        rfactor = self.rfactor.get()
-        if rfactor < 0 or rfactor > 1:
-            errors.append('Relaxation factor (lambda) must be between 0 and 1')
-        return errors
-
-    # --------------------------- UTLIS functions --------------------------------------------
-    def _updateItem(self, item, row):
-        newFn = row.getValue(md.MDL_IMAGE)
-        self.ix = self.ix + 1
-        newFn = newFn.split('@')[1]
-        item.setLocation(self.ix, newFn)
diff --git a/xmipp3/protocols/protocol_trigger_data.py b/xmipp3/protocols/protocol_trigger_data.py
index b729f81f0..8eff545c0 100644
--- a/xmipp3/protocols/protocol_trigger_data.py
+++ b/xmipp3/protocols/protocol_trigger_data.py
@@ -2,6 +2,7 @@
 # *
 # * Authors:     Tomas Majtner (tmajtner@cnb.csic.es)
 # *              David Maluenda (dmaluenda@cnb.csic.es)
+# *              Daniel Marchán (da.marchan@cnb.csic.es)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -35,22 +36,23 @@
 from pyworkflow.object import Set
 from pyworkflow.protocol.params import BooleanParam, IntParam, PointerParam, GT
 
+SIGNAL_FILENAME = "STOP_STREAM.TXT"
 
 class XmippProtTriggerData(EMProtocol):
     """
     Waits until certain number of images is prepared and then
     send them to output.
     It can be done in 3 ways:
-        - If *Send all particles to output?*' is _No_:
-            Once the number of images is reached, a setOfImages is returned and
-            the protocols finished (ending the streaming from this point).
-        - If *Send all particles to output?*' is _Yes_ and:
-            - If *Split particles to multiple sets?* is _Yes_:
+        - If "Send all items to output?" is _No_:
+            Once the number of items is reached, a setOfImages is returned and
+            the protocol finishes (ending the streaming from this point).
+        - If "Send all items to output?" is _Yes_ and:
+            - If "Split items to multiple sets?" is _Yes_:
                 Multiple closed outputs will be returned as soon as
-                the number of images is reached.
-            - If *Split particles to multiple sets?* is _No_:
+                the number of items is reached.
+            - If "Split items to multiple sets?" is _No_:
                 Only one output is returned and it is growing up in batches of
-                a certain number of images (completely in streaming).
+                a certain number of items (completely in streaming).
     """
     _label = 'trigger data'
     _lastUpdateVersion = VERSION_2_0
@@ -63,6 +65,12 @@ def _defineParams(self, form):
         form.addParam('inputImages', PointerParam,
                       pointerClass='SetOfImages',
                       label='Input images', important=True)
+        form.addParam('triggerWait', BooleanParam, default=False,
+                      label='Wait for signal to stop the stream?',
+                      help='If NO is selected, normal functionality.\n'
+                           'If YES is selected it will wait for a signal to stop the stream.'
+                           '\n For this option, select send all items to output with a '
+                           'minimum size of 1')
         form.addParam('outputSize', IntParam, default=10000,
                       label='Minimum output size',
                       help='How many particles need to be on input to '
@@ -79,6 +87,16 @@ def _defineParams(self, form):
                            '"Output size" are returned.\n'
                            'If NO is selected, only one open and growing output '
                            'is returned')
+        form.addParam('triggerSignal', BooleanParam, default=False,
+                      label='Send signal to stop a stream?',
+                      help='If NO is selected, normal functionality.\n'
+                           'If YES is selected it will send a signal to a connected Trigger data protocol.'
+                           '\n For this option, select the option send all items to output.')
+        form.addParam('triggerProt', PointerParam,
+                      pointerClass=self.getClassName(),
+                      condition='triggerSignal', allowsNull=True,
+                      label='Trigger data protocol',
+                      help='Select the trigger data protocol that you will send a signal to stop the stream.')
         form.addParam('delay', IntParam, default=10, label="Delay (sec)",
                       validators=[GT(3, "must be larger than 3sec.")],
                       help="Delay in seconds before checking new output")
@@ -149,11 +167,24 @@ def _checkNewOutput(self):
 
         if self.streamClosed:
             self.finished = True
-        elif not self.allImages.get():
+        elif not self.allImages.get() and not self.triggerSignal.get():
             self.finished = len(self.images) >= self.outputSize
         else:
             self.finished = False
 
+        # Send the signal to the connected protocol
+        if self.triggerSignal.get():
+            if len(self.images) >= self.outputSize:
+                print('Sending signal to stop the input trigger data protocol')
+                self.stopWait()
+
+        # Wait for trigger data signal
+        if self.triggerWait.get():
+            print('Waiting for signal to stop the stream')
+            if self.waitingHasFinished():
+                print('Stopped by received signal from a trigger data protocol')
+                self.finished = True
+
         outputStep = self._getFirstJoinStep()
         deps = []
         if self.finished:  # Unlock createOutputStep if finished all jobs
@@ -319,4 +350,23 @@ def getImagesType(self, letters='default'):
             return typeStr
 
     def getOututName(self):
-        return 'output%s' % self.getImagesType()
\ No newline at end of file
+        return 'output%s' % self.getImagesType()
+
+    def stopWait(self):
+        f = open(self._getStopConnectingFilename(), 'w')
+        f.close()
+
+    def _getStopConnectingFilename(self):
+        triggerProtocol = self.triggerProt.get()
+        if triggerProtocol is not None:
+            fileName = triggerProtocol._getExtraPath(SIGNAL_FILENAME)
+        else:
+            fileName = None
+
+        return fileName
+
+    def waitingHasFinished(self):
+        return os.path.exists(self._getStopWaitingFilename())
+
+    def _getStopWaitingFilename(self):
+        return self._getExtraPath(SIGNAL_FILENAME)
diff --git a/xmipp3/viewers/__init__.py b/xmipp3/viewers/__init__.py
index f425ae317..1c39faa8c 100644
--- a/xmipp3/viewers/__init__.py
+++ b/xmipp3/viewers/__init__.py
@@ -48,6 +48,8 @@
 from .viewer_volume_strain import XmippVolumeStrainViewer
 from .viewer_reconstruct_highres import XmippReconstructHighResViewer
 from .viewer_solid_angles import SolidAnglesViewer
+from .viewer_consensus_classes3D import XmippConsensusClasses3DViewer
+from .viewer_classes3D import XmippClasses3DViewer
 from .viewer_extract_asymmetric_unit import viewerXmippProtExtractUnit
 from .viewer_deepEMHancer import viewerXmippProtDeepVolPostProc
 from .viewer_local_sharpening import viewerXmippProtLocSharp
diff --git a/xmipp3/viewers/viewer_analyze_local_ctf.py b/xmipp3/viewers/viewer_analyze_local_ctf.py
index 28b7ddae6..c8cf3841a 100644
--- a/xmipp3/viewers/viewer_analyze_local_ctf.py
+++ b/xmipp3/viewers/viewer_analyze_local_ctf.py
@@ -28,11 +28,9 @@
 import numpy as np
 from pwem import emlib
 from os.path import exists
-
 from pyworkflow.viewer import (ProtocolViewer, DESKTOP_TKINTER, WEB_DJANGO)
 from pyworkflow.protocol.params import IntParam, LabelParam
 from pwem.viewers import ObjectView, showj
-
 from xmipp3.protocols.protocol_analyze_local_ctf import XmippProtAnalyzeLocalCTF
 from .plotter import XmippPlotter
 
@@ -56,37 +54,36 @@ def getData(self):
     def _defineParams(self, form):
         form.addSection(label='Visualization')
         form.addParam('displayR2', LabelParam, default=False, label='Display micrograph R2')
-        form.addParam('displayLocalDefocus', IntParam, default=1,
-                      label='Display local defocus of micrograph:',
+        form.addParam('displayLocalDefocus', IntParam, label='Display local defocus of micrograph:', allowsNull=True,
                       help="""Type the ID of the micrograph to see particle local defocus of the selected micrograph. 
-                      It is possible that not all the micrographs are available""")
+                      It is possible that not all the micrographs are available.
+                      Please check the ID of the micrographs in the output set of micrographs of this protocol.""")
 
     def _getVisualizeDict(self):
-        return {'displayR2' : self._viewR2,
+        return {'displayR2': self._viewR2,
                 'displayLocalDefocus': self._viewLocalDefocus,
                 }
 
     def _viewLocalDefocus(self, paramName=None):
         """display a 3D view of where the particles are placed in the micrograph taking as height the value estimated
-        for local defocues"""
-        views=[]
-        fnDefoci="%s"%(self.protocol._getExtraPath("micrographDefoci.xmd"))
+        for local defocus"""
+        views = []
+        fnDefoci = "%s" % self.protocol._getExtraPath("micrographDefoci.xmd")
         if exists(fnDefoci):
             try:
-                mdPoints = emlib.MetaData("mic_%d@%s"%(self.displayLocalDefocus.get(),fnDefoci))
+                mdPoints = emlib.MetaData("mic_%d@%s" % (self.displayLocalDefocus.get(), fnDefoci))
 
                 x = mdPoints.getColumnValues(emlib.MDL_XCOOR)
                 y = mdPoints.getColumnValues(emlib.MDL_YCOOR)
                 defocusA = mdPoints.getColumnValues(emlib.MDL_CTF_DEFOCUSA)
                 residuals = mdPoints.getColumnValues(emlib.MDL_CTF_DEFOCUS_RESIDUAL)
-
-                title="Micrograph %d defocus"%self.displayLocalDefocus.get()
+                title = "Micrograph %d defocus" % self.displayLocalDefocus.get()
                 xplotter = XmippPlotter(windowTitle=title)
                 a = xplotter.createSubPlot(title, 'x', 'y', projection='3d')
                 a.set_zlabel('Defocus')
                 a.scatter(x, y, defocusA, c='r', marker='o')
                 a.scatter(x, y, np.asarray(defocusA)-np.asarray(residuals), c='b', marker='^')
-                legends = ['Avg. defocus','Adjusted defocus']
+                legends = ['Avg. defocus', 'Adjusted defocus']
                 xplotter.showLegend(legends, loc=1)
                 views.append(xplotter)
             except Exception as e:
@@ -99,12 +96,8 @@ def _viewR2(self, paramName=None):
             obj = self.protocol.outputMicrographs
             fn = obj.getFileName()
             labels = 'id _filename _xmipp_ctfDefocusR2'
-            views.append(ObjectView(self._project, obj.strId(), fn,
-                                          viewParams={showj.ORDER: labels,
-                                                      showj.VISIBLE: labels,
-                                                      showj.RENDER: None,
-                                                      showj.MODE: showj.MODE_MD}))
+            views.append(ObjectView(self._project, obj.strId(), fn, viewParams={showj.ORDER: labels,
+                                                                                showj.VISIBLE: labels,
+                                                                                showj.RENDER: None,
+                                                                                showj.MODE: showj.MODE_MD}))
         return views
-
-
-
diff --git a/xmipp3/viewers/viewer_classes3D.py b/xmipp3/viewers/viewer_classes3D.py
new file mode 100644
index 000000000..bc58eba13
--- /dev/null
+++ b/xmipp3/viewers/viewer_classes3D.py
@@ -0,0 +1,161 @@
+# **************************************************************************
+# *
+# * Authors:     J.M. De la Rosa Trevin (jmdelarosa@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'
+# *
+# **************************************************************************
+
+"""
+This module implement the wrappers aroung Xmipp CL2D protocol
+visualization program.
+"""
+from pyworkflow.gui import IconButton, Icon
+from pyworkflow.gui.form import FormWindow
+from pyworkflow.viewer import (ProtocolViewer, DESKTOP_TKINTER, WEB_DJANGO)
+from pyworkflow.protocol.params import StringParam, LabelParam, IntParam
+
+from xmipp3.protocols.protocol_consensus_classes3D import XmippProtConsensusClasses3D
+from pwem.objects import Class3D
+
+import pickle
+import matplotlib.pyplot as plt
+
+FILE_OBJECTIVE_FDATA = 'ObjectiveFData.pkl'
+FILE_CLUSTERINGS = 'clusterings.pkl'
+FILE_ELBOWCLUSTERS = 'elbowclusters.pkl'
+
+
+class XmippClasses3DViewer(ProtocolViewer):
+    """ Visualization of results from the consensus classes 3D protocol
+    """
+    _label = 'viewer classes 3D'
+    _targets = [XmippProtConsensusClasses3D]
+    _environments = [DESKTOP_TKINTER, WEB_DJANGO]
+
+    def __init__(self, **kwargs):
+        ProtocolViewer.__init__(self, **kwargs)
+        self.nclusters, self.objFValues = self.getObjectiveFData()
+        self._elbowIdx = self.getElbowIndex()
+        self._allClusterings = self.getClusterings()
+
+    def getObjectiveFData(self):
+        loadPath = self.protocol._getExtraPath(FILE_OBJECTIVE_FDATA)
+        with open(loadPath, 'rb') as f:
+            self._ObjFData = pickle.load(f)
+        return self._ObjFData
+
+    def getClusterings(self):
+        loadPath = self.protocol._getExtraPath(FILE_CLUSTERINGS)
+        with open(loadPath, 'rb') as f:
+            self._allClusterings = pickle.load(f)
+        return self._allClusterings
+
+    def getElbowIndex(self):
+        loadPath = self.protocol._getExtraPath(FILE_ELBOWCLUSTERS)
+        with open(loadPath, 'rb') as f:
+            self.elbowIndex = pickle.load(f)
+        return self.elbowIndex
+
+    def _defineParams(self, form):
+        form.addSection(label='Visualization')
+        form.addParam('displayObjectiveFunction', LabelParam,
+                      label='Visualize Objective Function',
+                      help='Open a GUI to visualize the objective function for each number of clusters ')
+        form.addParam('chooseNumberOfClusters', IntParam, default=-1,
+                      label='Number of Clusters',
+                      help='Choose the final number of clusters in the consensus clustering. Press the eye')
+
+        #TODO: set an apply button to change the number of clusters instead of using the eye
+        '''command = self._chooseNumberOfClusters()
+        btn = IconButton(self.getWindow(), "Apply", Icon.ACTION_CONTINUE, command=command)
+        btn.config(relief="flat", activebackground=None, compound='left',
+                   fg='black', overrelief="raised")
+        btn.bind('<Button-1>')
+        btn.grid(row=2, column=2, sticky='nw')'''
+
+
+    def _getVisualizeDict(self):
+        return {'chooseNumberOfClusters': self._chooseNumberOfClusters,
+                'displayObjectiveFunction': self._visualizeObjectiveFunction}
+
+    def _visualizeObjectiveFunction(self, e=None):
+        nclusters, ob_values = self.nclusters, self.objFValues
+        elbowIdx = self._elbowIdx
+
+        plt.plot(nclusters, ob_values)
+        plt.scatter([nclusters[elbowIdx]], [ob_values[elbowIdx]], color='green')
+        plt.xlabel('Number of clusters')
+        plt.ylabel('Objective values')
+        plt.title('Objective values for each number of clusters')
+        plt.show()
+
+    def _chooseNumberOfClusters(self, e=None):
+        views=[]
+        nclusters = self.nclusters
+        elbow_nclust = nclusters[self._elbowIdx]
+
+        nclust = self.chooseNumberOfClusters.get()
+        if nclust == -1:
+            nclust = elbow_nclust
+        elif nclust > max(nclusters):
+            nclust = max(nclusters)
+
+        scipion_clustering = self.buildSetOfClasses(nclust)
+        self.protocol._defineOutputs(outputClasses=scipion_clustering)
+        for item in self.protocol.inputMultiClasses:
+            self.protocol._defineSourceRelation(item, scipion_clustering)
+
+        #TODO: esto es un objeto SetOfClasses3D pero no sabemos como lanzar un viewer con el sin tener el protocolo
+        #views.append(ObjectView(self._project, self.protocol.strId()))
+        return views
+
+    def buildSetOfClasses(self, nclust):
+        '''From a list of clusterings, nclust index it and that clustering
+        is converted into a scipion setOfClasses'''
+        clustIdx = self.nclusters.index(nclust)
+        clustering = self._allClusterings[clustIdx]
+
+        inputParticles = self.protocol.inputMultiClasses[0].get().getImages()
+        outputClasses = self.protocol._createSetOfClasses3D(inputParticles)
+        for classItem in clustering:
+            numOfPart = classItem[0]
+            partIds = classItem[1]
+            setRepId = classItem[2]
+            clsRepId = classItem[3]
+
+            setRep = self.protocol.inputMultiClasses[setRepId].get()
+            clRep = setRep[clsRepId]
+
+            newClass = Class3D()
+            # newClass.copyInfo(clRep)
+            newClass.setAcquisition(clRep.getAcquisition())
+            newClass.setRepresentative(clRep.getRepresentative())
+
+            outputClasses.append(newClass)
+
+            enabledClass = outputClasses[newClass.getObjId()]
+            enabledClass.enableAppend()
+            for itemId in partIds:
+                enabledClass.append(inputParticles[itemId])
+            outputClasses.update(enabledClass)
+
+        return outputClasses
diff --git a/xmipp3/viewers/viewer_consensus_classes3D.py b/xmipp3/viewers/viewer_consensus_classes3D.py
new file mode 100644
index 000000000..56f5d49f8
--- /dev/null
+++ b/xmipp3/viewers/viewer_consensus_classes3D.py
@@ -0,0 +1,291 @@
+# **************************************************************************
+# *
+# * Authors:     J.M. De la Rosa Trevin (jmdelarosa@cnb.csic.es)
+# *              Oier Lauzirika Zarrabeitia (oierlauzi@bizkaia.eu)
+# *
+# * 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'
+# *
+# **************************************************************************
+
+"""
+This module implement the wrappers aroung Xmipp CL2D protocol
+visualization program.
+"""
+from pyworkflow.viewer import (ProtocolViewer, DESKTOP_TKINTER, WEB_DJANGO)
+from pyworkflow.protocol.params import StringParam, LabelParam, IntParam
+from pyworkflow.protocol.params import GE, GT, LE, LT
+
+from pwem.viewers import TableView, ObjectView
+from pwem.viewers.showj import *
+
+from xmipp3.protocols.protocol_consensus_classes3D import XmippProtConsensusClasses3D
+
+from scipy.cluster import hierarchy
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+class XmippConsensusClasses3DViewer(ProtocolViewer):
+    """ Visualization of results from the consensus classes 3D protocol
+    """
+    _label = 'viewer consensus classes 3D'
+    _targets = [XmippProtConsensusClasses3D]
+    _environments = [DESKTOP_TKINTER, WEB_DJANGO]
+
+    def __init__(self, **kwargs):
+        ProtocolViewer.__init__(self, **kwargs)
+
+    def _defineParams(self, form):
+        # Particles section
+        form.addSection(label='Classes')
+        form.addParam('displayClasses', IntParam, label='Custom number of classes',
+                      validators=[GE(1)], default=1,
+                      help='Open a GUI to visualize the class of the given iteration. Warning: It takes time to open')
+
+        form.addParam('displayInitialClasses', LabelParam, label='Initial classes')
+        form.addParam('displayManualClasses', LabelParam, label='Manual number of classes')
+        form.addParam('displayOriginClasses', LabelParam, label='Number of classes by proximity to origin')
+        form.addParam('displayAngleClasses', LabelParam, label='Number of classes by angle')
+        form.addParam('displayPllClasses', LabelParam, label='Number of classes by profile likelihood')
+
+        # Graph section
+        form.addSection(label='Graphs')
+        form.addParam('displayDendrogram', LabelParam, label='Dendrogram',
+                        help='Open a GUI to visualize the dendogram each number of clusters')
+        form.addParam('displayDendrogramLog', LabelParam, label='Dendrogram (log)',
+                        help='Open a GUI to visualize the dendogram each number'
+                        ' of clusters with a logarithmic "y" axis')
+
+        form.addParam('displayObjectiveFunction', LabelParam, label='Objective Function',
+                        help='Open a GUI to visualize the objective function for each number of clusters')
+        form.addParam('displayObjectiveFunctionLog', LabelParam, label='Objective Function (log)',
+                        help='Open a GUI to visualize the objective function'
+                        ' for each number of clusters with a logarithmic "y" axis')
+
+        # Reference classification section (only if exists)
+        if self._readReferenceClassificationSizes() is not None:
+            form.addSection(label='Reference classification consensus')
+            form.addParam('displayReferenceClassificationSizePercentiles', LabelParam, label='Reference classification consensus size percentiles',
+                            help='Open a GUI to visualize a table with the most'
+                            ' common percentiles of the sizes of a random classification consensus')
+
+    def _getVisualizeDict(self):
+        return {
+            'displayClasses': self._visualizeClasses,
+            'displayInitialClasses': self._visualizeInitialClasses,
+            'displayManualClasses': self._visualizeManualClasses,
+            'displayOriginClasses': self._visualizeOriginClasses,
+            'displayAngleClasses': self._visualizeAngleClasses,
+            'displayPllClasses': self._visualizePllClasses,
+            'displayDendrogram': self._visualizeDendrogram,
+            'displayDendrogramLog': self._visualizeDendrogramLog,
+            'displayObjectiveFunction': self._visualizeObjectiveFunction,
+            'displayObjectiveFunctionLog': self._visualizeObjectiveFunctionLog,
+            'displayReferenceClassificationSizePercentiles': self._visualizeReferenceClassificationSizes,
+        }
+
+    # --------------------------- UTILS functions ------------------------------
+    def _getInputParticles(self):
+        return self.protocol._getInputParticles()
+
+    def _getInputClassifications(self):
+        return self.protocol._getInputClassifications()
+
+    def _readIntersections(self):
+        return self.protocol._readIntersections()
+
+    def _readClustering(self, iter):
+        return self.protocol._readClustering(iter)
+
+    def _readClusteringCount(self):
+        return self.protocol._readClusteringCount()
+
+    def _readAllClusterings(self):
+        return self.protocol._readAllClusterings()
+
+    def _readObjectiveValues(self):
+        return self.protocol._readObjectiveValues()
+    
+    def _readElbows(self):
+        return self.protocol._readElbows()
+
+    def _readReferenceClassificationSizes(self):
+        return self.protocol._readReferenceClassificationSizes()
+
+    def _readReferenceClassificationRelativeSizes(self):
+        return self.protocol._readReferenceClassificationRelativeSizes()
+
+    def _visualizeClasses(self, e):
+        # Read data
+        count = self.displayClasses.get()
+        particles = self._getInputParticles()
+        clustering = self._readClustering(count)
+        randomConsensusSizes = self._readReferenceClassificationSizes()
+        randomConsensusRelativeSizes = self._readReferenceClassificationRelativeSizes()
+
+        # Create the output classes and show them
+        outputClasses = self.protocol._createOutputClasses3D(
+            particles,
+            clustering, 
+            'viewer_'+str(count),
+            randomConsensusSizes, 
+            randomConsensusRelativeSizes
+        )
+
+        return self._showSetOfClasses3D(outputClasses)
+    
+    def _visualizeInitialClasses(self, e):
+        return self._showSetOfClasses3D(self.protocol.outputClasses_initial)
+
+    def _visualizeManualClasses(self, e):
+        return self._showSetOfClasses3D(self.protocol.outputClasses_manual)
+
+    def _visualizeOriginClasses(self, e):
+        return self._showSetOfClasses3D(self.protocol.outputClasses_origin)
+
+    def _visualizeAngleClasses(self, e):
+        return self._showSetOfClasses3D(self.protocol.outputClasses_angle)
+
+    def _visualizePllClasses(self, e):
+        return self._showSetOfClasses3D(self.protocol.outputClasses_pll)
+
+    def _visualizeDendrogram(self, e):
+        clusterings = self._readAllClusterings()
+        objectiveFunction = self._readObjectiveValues()
+
+        fig, ax = plt.subplots()
+        self._plotDendrogram(fig, ax, list(reversed(clusterings)), list(reversed(objectiveFunction)))
+        return [fig]
+        
+    def _visualizeDendrogramLog(self, e):
+        clusterings = self._readAllClusterings()
+        objectiveFunction = self._readObjectiveValues()
+
+        fig, ax = plt.subplots()
+        self._plotDendrogram(fig, ax, list(reversed(clusterings)), list(reversed(objectiveFunction)))
+        ax.set_yscale('log')
+        return [fig]
+        
+    def _visualizeObjectiveFunction(self, e):
+        objectiveFunction = self._readObjectiveValues()
+        elbows = self._readElbows()
+
+        fig, ax = plt.subplots()
+        self._plotObjectiveFunctionAndElbows(fig, ax, objectiveFunction, elbows)
+        return [fig]
+    
+    def _visualizeObjectiveFunctionLog(self, e):
+        objectiveFunction = self._readObjectiveValues()
+        elbows = self._readElbows()
+
+        fig, ax = plt.subplots()
+        self._plotObjectiveFunctionAndElbows(fig, ax, objectiveFunction[:-1], elbows)
+        ax.set_yscale('log')
+        return [fig]
+
+    def _visualizeReferenceClassificationSizes(self, e):
+        randomConsensusSizes = self._readReferenceClassificationSizes()
+        randomConsensusRelativeSizes = self._readReferenceClassificationRelativeSizes()
+
+        # Calculate the percentiles
+        percentiles = [1, 5, 10, 25, 50, 75, 90, 95, 99]
+        randomConsensusSizePercentiles = np.percentile(randomConsensusSizes, percentiles)
+        randomConsensusRelativeSizePercentiles = np.percentile(randomConsensusRelativeSizes, percentiles)
+        
+        # Create a table with the data
+        data = list(zip(percentiles, randomConsensusSizePercentiles, randomConsensusRelativeSizePercentiles))
+        header = ['Percentile (%)', 'Size percentile', 'Relative size percentile']
+        title = 'Reference classification consensus size percentiles'
+        return [TableView(header, data, title=title)]
+
+    def _showSetOfClasses3D(self, classes):
+        labels = 'enabled id _size _representative._filename _xmipp_classIntersectionSizePValue _xmipp_classIntersectionRelativeSizePValue'
+        labelRender = '_representative._filename'
+        return [ObjectView( self._project, classes.strId(), classes.getFileName(),
+                            viewParams={ORDER: labels,
+                                        VISIBLE: labels,
+                                        RENDER: labelRender,
+                                        SORT_BY: '_size desc',
+                                        MODE: MODE_MD})]
+
+    def _plotDendrogram(self, fig, ax, clusterings, obValues):
+        """ Plot the dendrogram from the objective functions of the merge
+            between the groups of images """
+
+        # Initialize required values
+        linkageMatrix = np.zeros((len(clusterings)-1, 4))
+        clsIds = np.arange(len(clusterings))
+        nIntersections = len(clusterings[0])
+
+        # Loop over each iteration of clustering
+        for i in range(len(clusterings)-1):
+            # Find the sets that were merged
+            clsMerged = []
+            assert(len(clusterings[i]) == len(clsIds))
+            for cls, clsId in zip(clusterings[i], clsIds):
+                if cls not in clusterings[i+1]:
+                    clsMerged.append(clsId)
+            assert(len(clsMerged) == 2)
+            clsMerged = np.array(clsMerged)
+            
+            # Find original number of sets within new set
+            nCls = list(map(lambda x : linkageMatrix[x,3] if (x >= 0) else 1, clsMerged - nIntersections))
+
+            # Create linkage matrix
+            linkageMatrix[i, 0:2] = clsMerged  # ids of merged sets
+            linkageMatrix[i, 2] = obValues[i+1]   # objective function as distance
+            linkageMatrix[i, 3] = sum(nCls)  # total number of original sets
+
+            # Change set ids to reflect new set of clusters
+            for id in clsMerged:
+                clsIds = np.delete(clsIds, np.argwhere(clsIds == id))
+            clsIds = np.append(clsIds, len(clusterings)+i)
+
+        # Plot resulting dendrogram
+        hierarchy.dendrogram(linkageMatrix, ax=ax)
+        ax.set_title('Dendrogram')
+        ax.set_xlabel('sets ids')
+        ax.set_ylabel('objective function')
+        ax.set_ylim([np.min(linkageMatrix[:, 2]), np.max(linkageMatrix[:, 2])])
+
+    def _plotObjectiveFunctionAndElbows(self, fig, ax, obValues, elbows):
+        """ Plots the objective function and elbows """
+
+        # Shorthands for some variables
+        numClusters = np.arange(1, 1+len(obValues))
+
+        # Plot obValues vs numClusters
+        ax.plot(numClusters, obValues)
+
+        # Show elbows as scatter points
+        for key, value in elbows.items():
+            value -= 1 # Value uses 1-based indexing
+            label = key + ': ' + str(numClusters[value])
+            ax.scatter([numClusters[value]], [obValues[value]], label=label, color='green')
+
+        # Configure the figure
+        ax.legend()
+        ax.set_xlabel('Number of clusters')
+        ax.set_ylabel('Objective values')
+        ax.set_title('Objective values for each number of clusters')
+        ax.set_ylim([np.min(obValues), np.max(obValues)])
+        
diff --git a/xmipp3/xmipp_logo_devel.png b/xmipp3/xmipp_logo_devel.png
new file mode 100644
index 000000000..0845cac12
Binary files /dev/null and b/xmipp3/xmipp_logo_devel.png differ
diff --git a/xmipp3/xmipp_logo_original.png b/xmipp3/xmipp_logo_original.png
new file mode 100644
index 000000000..1ada98399
Binary files /dev/null and b/xmipp3/xmipp_logo_original.png differ


Test patch:
diff --git a/xmipp3/tests/test_protocol_consensus_classes3D.py b/xmipp3/tests/test_protocol_consensus_classes3D.py
index d3682b004..9d110ca0f 100644
--- a/xmipp3/tests/test_protocol_consensus_classes3D.py
+++ b/xmipp3/tests/test_protocol_consensus_classes3D.py
@@ -31,13 +31,13 @@
 from pyworkflow.tests import BaseTest, setupTestProject, DataSet
 from pyworkflow.utils import greenStr
 from pwem.objects import Class3D
+from pwem.protocols import ProtClassify3D
 from pwem.protocols import ProtImportParticles
 
 from xmipp3.protocols import XmippProtConsensusClasses3D
 
 
 class TestConsensusClasses3D(BaseTest):
-
     @classmethod
     def setUpClass(cls):
         """Prepare the data that we will use later on."""
@@ -49,26 +49,19 @@ def setUpClass(cls):
         cls.xmippDataTest = DataSet.getDataSet('xmipp_tutorial')
 
 
-        def _importClasses(label, numClasses, randomSeed, numPart=None):
+        def _createClasses(label, particles, numClasses, randomSeed, numPart=None):
             """ We import a set of classes with an alterate import particles
                 because there is not a import classes protocol
             """
-            pImpClasses = cls.proj.newProtocol(ProtImportParticles,
-                               filesPath=cls.xmippDataTest.getFile('particles'),
-                               samplingRate=3.5)
-
-            pImpClasses.setObjLabel('Import %s' % label)
-
-            # we launch the protocol to obtain the particles
+            # Create a dummy class with the particles
+            pImpClasses = cls.proj.newProtocol(ProtClassify3D)
+            pImpClasses.setObjLabel(label)
             cls.proj.launchProtocol(pImpClasses, wait=True)
 
-            # fake importing classes handmade
-            partSet = pImpClasses.outputParticles
+            setOfClasses = pImpClasses._createSetOfClasses3D(particles)
 
-            setOfClasses = pImpClasses._createSetOfClasses3D(partSet)
-
-            numOfPart = partSet.getSize() if numPart is None else numPart
-            partIds = list(partSet.getIdSet())
+            numOfPart = particles.getSize() if numPart is None else numPart
+            partIds = list(particles.getIdSet())
             m = int(numOfPart/numClasses)
             
             # random shuffle with a certain seed to get always the same classes
@@ -78,8 +71,8 @@ def _importClasses(label, numClasses, randomSeed, numPart=None):
                 currIds = partIds[clInx*m:(clInx+1)*m]
 
                 newClass = Class3D()
-                newClass.copyInfo(partSet)
-                newClass.setAcquisition(partSet.getAcquisition())
+                newClass.copyInfo(particles)
+                newClass.setAcquisition(particles.getAcquisition())
                 # newClass.setRepresentative(clRep.getRepresentative())
 
                 setOfClasses.append(newClass)
@@ -87,7 +80,7 @@ def _importClasses(label, numClasses, randomSeed, numPart=None):
                 enabledClass = setOfClasses[newClass.getObjId()]
                 enabledClass.enableAppend()
                 for itemId in currIds:
-                    item = partSet[itemId]
+                    item = particles[itemId]
                     enabledClass.append(item)
 
                 setOfClasses.update(enabledClass)
@@ -100,12 +93,23 @@ def _importClasses(label, numClasses, randomSeed, numPart=None):
 
             return pImpClasses.outputClasses
 
-        cls.set1 = _importClasses('3 Classes I', 3, 65)
-        cls.set2 = _importClasses('3 Classes II', 3, 123)
-        cls.set3 = _importClasses('3 Classes III', 3, 256)
-        cls.set4 = _importClasses('3 Classes IV (not all part)', 3, 568, 70)
-        cls.set5 = _importClasses('5 Classes I', 5, 745)
-        cls.set6 = _importClasses('5 Classes II', 5, 1025)
+        # Import particles
+        pImpParticles = cls.proj.newProtocol(ProtImportParticles,
+                            filesPath=cls.xmippDataTest.getFile('particles'),
+                            samplingRate=3.5)
+        pImpParticles.setObjLabel('Import particles')
+
+        # Launch the protocol to obtain the particles
+        cls.proj.launchProtocol(pImpParticles, wait=True)
+
+        # fake importing classes handmade
+        particles = pImpParticles.outputParticles
+
+        cls.set1 = _createClasses('3 Classes I', particles, 3, 65)
+        cls.set2 = _createClasses('3 Classes II', particles, 3, 123)
+        cls.set3 = _createClasses('3 Classes III', particles, 3, 256)
+        cls.set4 = _createClasses('5 Classes I', particles, 5, 745)
+        cls.set5 = _createClasses('5 Classes II', particles, 5, 1025)
 
 
     def checkIntersections(self, setOfIntersections, classId, partIds):
@@ -133,16 +137,16 @@ def testConsensus(self):
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set1, 
-                                                              self.set2, 
-                                                              self.set3])
+                                           inputClassifications=[self.set1, 
+                                                                 self.set2, 
+                                                                 self.set3])
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
+        setOfIntersections = pConsClass.outputClasses_initial
 
         # some general assertions
         self.assertIsNotNone(setOfIntersections,
                              "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 27,
+        self.assertEqual(setOfIntersections.getSize(), 26,
                          "The number of the outputClasses is wrong")
         self.checkPopulation(setOfIntersections, 75)
 
@@ -152,18 +156,18 @@ def testConsensus2(self):
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set1, 
-                                                              self.set2, 
-                                                              self.set4])
+                                           inputClassifications=[self.set1, 
+                                                                 self.set2, 
+                                                                 self.set4])
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
+        setOfIntersections = pConsClass.outputClasses_initial
 
         # some general assertions
         self.assertIsNotNone(setOfIntersections,
                              "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 27,
+        self.assertEqual(setOfIntersections.getSize(), 38,
                          "The number of the outputClasses is wrong")
-        self.checkPopulation(setOfIntersections, 69)
+        self.checkPopulation(setOfIntersections, 75)
 
 
     def testConsensus3(self):
@@ -171,18 +175,18 @@ def testConsensus3(self):
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set4, 
-                                                              self.set2, 
-                                                              self.set1])
+                                           inputClassifications=[self.set4, 
+                                                                 self.set2, 
+                                                                 self.set1])
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
+        setOfIntersections = pConsClass.outputClasses_initial
 
         # some general assertions
         self.assertIsNotNone(setOfIntersections,
                              "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 27,
+        self.assertEqual(setOfIntersections.getSize(), 38,
                          "The number of the outputClasses is wrong")
-        self.checkPopulation(setOfIntersections, 69)
+        self.checkPopulation(setOfIntersections, 75)
 
 
     def testConsensus4(self):
@@ -190,54 +194,52 @@ def testConsensus4(self):
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set5, 
-                                                              self.set2, 
-                                                              self.set1])
+                                           inputClassifications=[self.set5, 
+                                                                 self.set2, 
+                                                                 self.set1])
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
+        setOfIntersections = pConsClass.outputClasses_initial
 
         # some general assertions
         self.assertIsNotNone(setOfIntersections,
                              "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 45,
+        self.assertEqual(setOfIntersections.getSize(), 36,
                          "The number of the outputClasses is wrong")
         self.checkPopulation(setOfIntersections, 75)
 
-
-    def testConsensus5(self):
-        print("\n", greenStr(" Test Consensus with two sets".center(75, '-')))
+    def testConsensus6(self):
+        print("\n", greenStr(" Test Consensus with four sets".center(75, '-')))
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set5, 
-                                                              self.set6])
+                                           inputClassifications=[self.set1,
+                                                                 self.set2,
+                                                                 self.set3,
+                                                                 self.set5 ])
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
+        setOfIntersections = pConsClass.outputClasses_initial
 
         # some general assertions
         self.assertIsNotNone(setOfIntersections,
                              "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 25,
+        self.assertEqual(setOfIntersections.getSize(), 61,
                          "The number of the outputClasses is wrong")
         self.checkPopulation(setOfIntersections, 75)
 
-
-    def testConsensus6(self):
-        print("\n", greenStr(" Test Consensus with four sets".center(75, '-')))
+    def testConsensus7(self):
+        print("\n", greenStr("Test Consensus with four sets and a random reference classification".center(75, '-')))
 
         # preparing and launching the protocol
         pConsClass = self.proj.newProtocol(XmippProtConsensusClasses3D,
-                                           inputMultiClasses=[self.set1,
-                                                              self.set2,
-                                                              self.set5, 
-                                                              self.set6])
+                                           inputClassifications=[self.set1,
+                                                                 self.set2,
+                                                                 self.set3,
+                                                                 self.set5 ],
+                                           randomClassificationCount=1000 )
         self.proj.launchProtocol(pConsClass, wait=True)
-        setOfIntersections = pConsClass.outputClasses
-
-        # some general assertions
-        self.assertIsNotNone(setOfIntersections,
-                             "There was some problem with the output")
-        self.assertEqual(setOfIntersections.getSize(), 225,
-                         "The number of the outputClasses is wrong")
-        self.checkPopulation(setOfIntersections, 75)
+        setOfIntersections = pConsClass.outputClasses_initial
 
+        # Ensure that the intersection size has been written
+        for cls in setOfIntersections:
+            self.assertIsNotNone(cls._xmipp_classIntersectionSizePValue)
+            self.assertIsNotNone(cls._xmipp_classIntersectionRelativeSizePValue)
\ No newline at end of file
diff --git a/xmipp3/tests/test_protocol_reconstruct_fourier.py b/xmipp3/tests/test_protocol_reconstruct_fourier.py
new file mode 100644
index 000000000..7fcd6a879
--- /dev/null
+++ b/xmipp3/tests/test_protocol_reconstruct_fourier.py
@@ -0,0 +1,83 @@
+# ***************************************************************************
+# * Authors:     David Herreros Calero (dherreros@cnb.csic.es)
+# *
+# *
+# * 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, setupTestProject
+
+from pwem import emlib
+
+from xmipp3.protocols import XmippProtReconstructFourier, \
+                             XmippProtCreateMask3D, XmippProtCreateGallery
+
+
+class TestReconstructFourier(BaseTest):
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+
+    def checkOutput(self, prot, outName, conditions=[]):
+        """ Check that an ouput was generated and
+        the condition is valid.
+        """
+        o = getattr(prot, outName, None)
+        locals()[outName] = o
+        self.assertIsNotNone(o, "Output: %s is None" % outName)
+        for c in conditions:
+            self.assertTrue(eval(c), 'Condition failed: ' + c)
+
+    def createPhantomCube(self):
+        spherePhantom = self.newProtocol(XmippProtCreateMask3D,
+                                         objLabel='Target volume',
+                                         source=1,
+                                         geo=1,
+                                         size=128)
+        self.launchProtocol(spherePhantom)
+        self.assertIsNotNone(spherePhantom.outputMask)
+        return spherePhantom
+
+    def createProjections(self, volume):
+        projections = self.newProtocol(XmippProtCreateGallery,
+                                       objLabel='Projections',
+                                       inputVolume=volume,
+                                       rotStep=10,
+                                       tiltStep=10)
+        self.launchProtocol(projections)
+        self.assertIsNotNone(projections.outputReprojections)
+        return projections
+
+    def test_reconstruct_fourier(self):
+        cubePhantom = self.createPhantomCube()
+        projections = self.createProjections(cubePhantom.outputMask)
+
+        validateProt = self.newProtocol(XmippProtReconstructFourier,
+                                objLabel='Fourier reconstruction',
+                                inputParticles=projections.outputReprojections)
+
+        self.launchProtocol(validateProt)
+        self.checkOutput(validateProt, 'outputVolume')
+
+        # Check if cubes are similar
+        sphere_phantom = emlib.Image(cubePhantom.outputMask.getFileName())
+        sphere_reconstructed = emlib.Image(validateProt.outputVolume.getFileName())
+        corr = sphere_phantom.correlation(sphere_reconstructed)
+        self.assertAlmostEqual(corr, 0.98, delta=0.01, msg="There was an error during the reconstruction process")
diff --git a/xmipp3/tests/test_protocols_deepVolPostprocessing.py b/xmipp3/tests/test_protocols_deepVolPostprocessing.py
index 6d96e61af..9eb723076 100644
--- a/xmipp3/tests/test_protocols_deepVolPostprocessing.py
+++ b/xmipp3/tests/test_protocols_deepVolPostprocessing.py
@@ -124,7 +124,7 @@ def testDeepVolPostpro2(self):
                                         useHalfMapsInsteadVol=True,
                                         halfMapsAttached=False,
                                         inputHalf1=self.protImportHalf1.outputVolume,
-                                        inputHalf2=self.protImportHalf1.outputVolume,
+                                        inputHalf2=self.protImportHalf2.outputVolume,
                                         normalization= XmippProtDeepVolPostProc.NORMALIZATION_AUTO
                                    )
         self.launchProtocol(deepPostProc)
diff --git a/xmipp3/tests/test_protocols_local_defocus.py b/xmipp3/tests/test_protocols_local_defocus.py
new file mode 100644
index 000000000..7643f868a
--- /dev/null
+++ b/xmipp3/tests/test_protocols_local_defocus.py
@@ -0,0 +1,144 @@
+# ***************************************************************************
+# *
+# * Authors:     Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *
+# * 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
+import pyworkflow.tests as tests
+
+from xmipp3.protocols import XmippProtPhantom, XmippProtCreateGallery, XmippProtSimulateCTF, \
+    XmippProtLocalCTF, XmippProtConsensusLocalCTF, XmippProtAnalyzeLocalCTF
+from pwem.protocols import ProtImportMicrographs, ProtImportParticles, ProtImportVolumes
+
+
+class TestXmippLocalDefocusEstimation(BaseTest):
+    """ Testing local defocus estimation """
+    @classmethod
+    def setUpClass(cls):
+        tests.setupTestProject(cls)
+        cls.protCreatePhantom = cls.newProtocol(XmippProtPhantom)
+        cls.launchProtocol(cls.protCreatePhantom)
+        cls.assertIsNotNone(cls.protCreatePhantom.getFiles(), "There was a problem with phantom creation")
+        cls.protCreateGallery = cls.newProtocol(XmippProtCreateGallery,
+                                                inputVolume=cls.protCreatePhantom.outputVolume,
+                                                rotStep=15.0,
+                                                tiltStep=90.0)
+        cls.launchProtocol(cls.protCreateGallery)
+        cls.assertIsNotNone(cls.protCreateGallery.getFiles(), "There was a problem with create gallery")
+
+    def testXmippLocalDefocus(self):
+        protSimulateCTF = self.newProtocol(XmippProtSimulateCTF)
+        protSimulateCTF.inputParticles.set(self.protCreateGallery.outputReprojections)
+        self.launchProtocol(protSimulateCTF)
+        self.assertIsNotNone(protSimulateCTF.outputParticles, "There was a problem with CTF simulation")
+
+        protEstimateLocalDefocus = self.newProtocol(XmippProtLocalCTF,
+                                                    inputSet=protSimulateCTF.outputParticles,
+                                                    inputVolume=self.protCreatePhantom.outputVolume)
+        self.launchProtocol(protEstimateLocalDefocus)
+        self.assertIsNotNone(protEstimateLocalDefocus.outputParticles.getFiles(),
+                             "There was a problem with local CTF estimation")
+        self.assertEqual(protEstimateLocalDefocus.outputParticles.getDim(), (40, 40, 181),
+                         "There was a problem with size of set of particles")
+        self.assertEqual(protEstimateLocalDefocus.outputParticles.getFirstItem().getSamplingRate(), 4,
+                         "There was a problem with the sampling rate value of output particles")
+        self.assertEqual(protEstimateLocalDefocus.outputParticles.getFirstItem().getCTF().getDefocusU(),
+                         protEstimateLocalDefocus.outputParticles.getFirstItem().getCTF().getDefocusV(),
+                         "Defocus U and defocus V should be equal if sameDefocus option is set to True")
+
+
+class TestXmippLocalDefocusEstimationConsensus(TestXmippLocalDefocusEstimation):
+    """ Testing local defocus estimation consensus """
+
+    def testXmippLocalDefocusConsensus(self):
+        protSimulateCTF2 = self.newProtocol(XmippProtSimulateCTF,
+                                            inputParticles=self.protCreateGallery.outputReprojections)
+        self.launchProtocol(protSimulateCTF2)
+        self.assertIsNotNone(protSimulateCTF2.outputParticles, "There was a problem with second CTF simulation")
+
+        protEstimateLocalDefocus2 = self.newProtocol(XmippProtLocalCTF,
+                                                     inputSet=protSimulateCTF2.outputParticles,
+                                                     inputVolume=self.protCreatePhantom.outputVolume)
+        self.launchProtocol(protEstimateLocalDefocus2)
+        self.assertIsNotNone(protEstimateLocalDefocus2.outputParticles.getFiles(),
+                             "There was a problem with second CTF estimation")
+
+        protSimulateCTF3 = self.newProtocol(XmippProtSimulateCTF,
+                                            inputParticles=self.protCreateGallery.outputReprojections)
+        self.launchProtocol(protSimulateCTF3)
+        self.assertIsNotNone(protSimulateCTF3.outputParticles, "There was a problem with third CTF simulation")
+
+        protEstimateLocalDefocus3 = self.newProtocol(XmippProtLocalCTF,
+                                                     inputSet=protSimulateCTF3.outputParticles,
+                                                     inputVolume=self.protCreatePhantom.outputVolume)
+        self.launchProtocol(protEstimateLocalDefocus3)
+        self.assertIsNotNone(protEstimateLocalDefocus3.outputParticles.getFiles(),
+                             "There was a problem with third CTF estimation")
+
+        protConsensusLocalDefocus = self.newProtocol(XmippProtConsensusLocalCTF,
+                                                     inputSet=protSimulateCTF2.outputParticles,
+                                                     inputSets=[protEstimateLocalDefocus2.outputParticles,
+                                                                protEstimateLocalDefocus3.outputParticles])
+        self.launchProtocol(protConsensusLocalDefocus)
+        self.assertIsNotNone(protConsensusLocalDefocus.outputParticles.getFiles(),
+                             "There was a problem with consensus of local defocus estimation")
+        self.assertEqual(protConsensusLocalDefocus.outputParticles.getDim(), (40, 40, 181),
+                         "There was a problem with size of set of particles after consensus")
+        self.assertEqual(protConsensusLocalDefocus.outputParticles.getFirstItem().getSamplingRate(), 4,
+                         "There was a problem with the sampling rate value of output consensus particles")
+
+
+class TestXmippAnalyzeLocalDefocus(BaseTest):
+    """ Testing protocol analyze local defocus """
+
+    @classmethod
+    def setUpClass(cls):
+        tests.setupTestProject(cls)
+        cls.dataset = DataSet.getDataSet('xmipp_tutorial')
+        cls.micsFn = cls.dataset.getFile('micrographs/BPV_1386.mrc')
+        cls.partsFn = cls.dataset.getFile('particles/BPV_particles.sqlite')
+        cls.volFn = cls.dataset.getFile('volumes/BPV_scale_filtered_windowed_64.vol')
+
+    def testXmippAnalyzeLocalDefocus(self):
+        protImportMics = self.newProtocol(ProtImportMicrographs,
+                                          filesPath=self.micsFn,
+                                          samplingRate=1.237,
+                                          voltage=300)
+        self.launchProtocol(protImportMics)
+        self.assertIsNotNone(protImportMics.outputMicrographs, "There was a problem with the import of micrographs")
+
+        protImportParticles = self.newProtocol(ProtImportParticles,
+                                               importFrom=4,
+                                               sqliteFile=self.partsFn,
+                                               samplingRate=1.237)
+        self.launchProtocol(protImportParticles)
+        self.assertIsNotNone(protImportParticles.getFiles(), "There was a problem with the import of particles")
+
+        protAnalyzeLocalDefocus = self.newProtocol(XmippProtAnalyzeLocalCTF,
+                                                   inputMics=protImportMics.outputMicrographs,
+                                                   inputSet=protImportParticles.outputParticles)
+        self.launchProtocol(protAnalyzeLocalDefocus)
+        self.assertIsNotNone(protAnalyzeLocalDefocus.outputMicrographs.getFiles(),
+                             "There was a problem with analysis of local CTF estimation")
+        self.assertEqual(protAnalyzeLocalDefocus.outputMicrographs.getDim(), (9216, 9441, 1),
+                         "There was a problem with size of set of particles after analyze local defocus")
+        self.assertEqual(protAnalyzeLocalDefocus.outputMicrographs.getFirstItem().getSamplingRate(), 1.237,
+                         "There was a problem with the sampling rate value of output analyze local defocus particles")
diff --git a/xmipp3/tests/test_protocols_metaprotocol_heterogeneity.py b/xmipp3/tests/test_protocols_metaprotocol_heterogeneity.py
index 3162f224e..794f7c48a 100644
--- a/xmipp3/tests/test_protocols_metaprotocol_heterogeneity.py
+++ b/xmipp3/tests/test_protocols_metaprotocol_heterogeneity.py
@@ -24,6 +24,7 @@
 # *
 # **************************************************************************
 
+'''
 from pyworkflow.tests import BaseTest, setupTestProject, DataSet
 from pwem.protocols import ProtImportParticles, ProtImportVolumes, ProtSubSet
 
@@ -74,4 +75,4 @@ def test(self):
         self.launchProtocol(protHeterogeneity)
         self.assertFalse(protHeterogeneity.isFailed(), 'Metaprotocol Heterogeneity has failed.')
 
-
+'''
diff --git a/xmipp3/tests/test_protocols_xmipp_2d.py b/xmipp3/tests/test_protocols_xmipp_2d.py
index 1afab6fc6..340d36c2c 100644
--- a/xmipp3/tests/test_protocols_xmipp_2d.py
+++ b/xmipp3/tests/test_protocols_xmipp_2d.py
@@ -1217,7 +1217,7 @@ def setUpClass(cls):
     
     def test_particles1(self):
         print("Run Compare Reprojections from classes")
-        prot = self.newProtocol(XmippProtCompareReprojections, 
+        prot = self.newProtocol(XmippProtCompareReprojections,
                                         symmetryGroup="d6", numberOfMpi=5)
         prot.inputSet.set(self.protClassify.outputClasses)
         prot.inputVolume.set(self.protImportVol.outputVolume)
@@ -1226,7 +1226,7 @@ def test_particles1(self):
 
     def test_particles2(self):
         print("Run Compare Reprojections from averages")
-        prot = self.newProtocol(XmippProtCompareReprojections, 
+        prot = self.newProtocol(XmippProtCompareReprojections,
                                         symmetryGroup="d6", numberOfMpi=5)
         prot.inputSet.set(self.protImportAvgs.outputAverages)
         prot.inputVolume.set(self.protImportVol.outputVolume)
@@ -1235,7 +1235,7 @@ def test_particles2(self):
 
     def test_particles3(self):
         print("Run Compare Reprojections from projections with angles")
-        prot = self.newProtocol(XmippProtCompareReprojections, 
+        prot = self.newProtocol(XmippProtCompareReprojections,
                                         symmetryGroup="d6", numberOfMpi=5)
         prot.inputSet.set(self.protProjMatch.outputParticles)
         prot.inputVolume.set(self.protImportVol.outputVolume)
@@ -1249,7 +1249,6 @@ def test_particles4(self):
                                 symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
         prot.inputSet.set(self.protClassify.outputClasses)
         prot.inputVolume.set(self.protImportVol.outputVolume)
-        prot.maskVol.set(self.protCreateMask.outputMask)
         self.launchProtocol(prot)
         self.assertIsNotNone(prot.reprojections,
                              "There was a problem with Compare Reprojections from classes evaluating residuals")
@@ -1260,7 +1259,6 @@ def test_particles5(self):
                                 symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
         prot.inputSet.set(self.protImportAvgs.outputAverages)
         prot.inputVolume.set(self.protImportVol.outputVolume)
-        prot.maskVol.set(self.protCreateMask.outputMask)
         self.launchProtocol(prot)
         self.assertIsNotNone(prot.reprojections,
                              "There was a problem with Compare Reprojections from averages evaluating residuals")
@@ -1271,7 +1269,6 @@ def test_particles6(self):
                                 symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
         prot.inputSet.set(self.protProjMatch.outputParticles)
         prot.inputVolume.set(self.protImportVol.outputVolume)
-        prot.maskVol.set(self.protCreateMask.outputMask)
         self.launchProtocol(prot)
         self.assertIsNotNone(prot.reprojections, "There was a problem with Compare Reprojections from projections"
                                                  " with angles evaluating residuals")
@@ -1384,7 +1381,6 @@ def test_CorrectWiener(self):
         self.launchProtocol(protCorrect)
         self.assertIsNotNone(protCorrect.outputParticles, "There was a problem with Wiener Correction")
 
-
 if __name__ == "__main__":
     if len(sys.argv) > 1:
         className = sys.argv[1]
diff --git a/xmipp3/tests/test_protocols_xmipp_3d.py b/xmipp3/tests/test_protocols_xmipp_3d.py
index 04b515a77..78223d91a 100644
--- a/xmipp3/tests/test_protocols_xmipp_3d.py
+++ b/xmipp3/tests/test_protocols_xmipp_3d.py
@@ -48,9 +48,17 @@
                                                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"
+MSG_WRONG_SAMPLING = "There was a problem with the sampling rate value of the output "
+MSG_WRONG_SIZE = "There was a problem with the size of the output "
+MSG_WRONG_DIM = "There was a problem with the dimensions of the output "
+MSG_WRONG_MASK = "There was a problem with create mask from volume"
+MSG_WRONG_ALIGNMENT = "There was a problem with the alignment of the output "
+MSG_WRONG_SHIFT = "There was a problem with output shift "
+MSG_WRONG_GALLERY = "There was a problem with the gallery creation"
+MSG_WRONG_ROTATION = "There was a problem with the rotation"
+
+
+
 
 
 class TestXmippBase(BaseTest):
@@ -124,7 +132,7 @@ def testCreateMask(self):
         protMask1.setObjLabel('threshold mask')
         self.launchProtocol(protMask1)
         self.assertIsNotNone(protMask1.outputMask,
-                             MSG_MASK_ERROR)
+                             MSG_WRONG_MASK)
 
 
         print("Run create segment mask from volume")
@@ -135,7 +143,7 @@ def testCreateMask(self):
         protMask2.setObjLabel('segmentation automatic')
         self.launchProtocol(protMask2)
         self.assertIsNotNone(protMask2.outputMask,
-                             MSG_MASK_ERROR)
+                             MSG_WRONG_MASK)
 
     
         print("Run create mask from another mask")
@@ -822,7 +830,7 @@ def testExhaustiveBinaryMask(self):
         protMask.inputVolume.set(self.protImport1.outputVolume)
         self.launchProtocol(protMask)
         self.assertIsNotNone(protMask.outputMask,
-                             MSG_MASK_ERROR)
+                             MSG_WRONG_MASK)
         protAlign = self.newProtocol(XmippProtAlignVolume,
                                      inputReference=self.protImport1.outputVolume,
                                      alignmentAlgorithm=ALIGN_ALGORITHM_EXHAUSTIVE,
@@ -1103,9 +1111,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=MSG_WRONG_SAMPLING)
+                               msg=(MSG_WRONG_SAMPLING, "volume"))
         self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1,
-                               msg=MSG_WRONG_SIZE)
+                               msg=(MSG_WRONG_SIZE, "volume"))
         
     def testXmippPdbConvertFromObj(self):
         print("Run convert a pdb from import")
@@ -1122,9 +1130,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=MSG_WRONG_SAMPLING)
+                               msg=(MSG_WRONG_SAMPLING, "volume"))
         self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], protConvert.size_z.get(), places=1,
-                               msg=MSG_WRONG_SIZE)
+                               msg=(MSG_WRONG_SIZE, "volume"))
 
     def testXmippPdbConvertFromFn(self):
         print("Run convert a pdb from file")
@@ -1132,8 +1140,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=MSG_WRONG_SAMPLING)
-        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], 48, places=1, msg=MSG_WRONG_SIZE)
+                               msg=(MSG_WRONG_SAMPLING, "volume"))
+        self.assertAlmostEqual(protConvert.outputVolume.getDim()[0], 48, places=1, msg=(MSG_WRONG_SIZE, "volume"))
 
 
 class TestXmippValidateNonTilt(TestXmippBase):
@@ -1248,9 +1256,9 @@ def testXmippVolSub(self):
                                       vol2=protImportVol2.outputVolume,
                                       masks=False)
         self.launchProtocol(protVolAdj)
-        self.assertIsNotNone(protVolAdj.outputVolume,
-                             "There was a problem with Volumes adjust")
-
+        self.assertIsNotNone(protVolAdj.outputVolume, "There was a problem with Volumes adjust")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
         protVolAdjNoE = self.newProtocol(XmippProtVolAdjust,
                                          vol1=protImportVol1.outputVolume,
                                          vol2=protImportVol2.outputVolume,
@@ -1259,6 +1267,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolAdjNoE)
         self.assertIsNotNone(protVolAdjNoE.outputVolume,
                              "There was a problem with Volumes adjust without computing energy")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         protVolAdjNoRadAvg = self.newProtocol(XmippProtVolAdjust,
                                               vol1=protImportVol1.outputVolume,
@@ -1268,6 +1278,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolAdjNoRadAvg)
         self.assertIsNotNone(protVolAdjNoRadAvg.outputVolume,
                              "There was a problem with Volumes adjust without radial average")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         print("Run volume subtraction")
         protVolSub = self.newProtocol(XmippProtVolSubtraction,
@@ -1278,6 +1290,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolSub)
         self.assertIsNotNone(protVolSub.outputVolume,
                              "There was a problem with Volumes subtraction")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         protVolSubNoE = self.newProtocol(XmippProtVolSubtraction,
                                          vol1=protImportVol1.outputVolume,
@@ -1288,6 +1302,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolSubNoE)
         self.assertIsNotNone(protVolSubNoE.outputVolume,
                              "There was a problem with Volumes subtraction without computing energy")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         protVolSubMask = self.newProtocol(XmippProtVolSubtraction,
                                       vol1=protImportVol1.outputVolume,
@@ -1298,6 +1314,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolSubMask)
         self.assertIsNotNone(protVolSubMask.outputVolume,
                              "There was a problem with Volumes subtraction with masks")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         protVolSubRadAvg = self.newProtocol(XmippProtVolSubtraction,
                                       vol1=protImportVol1.outputVolume,
@@ -1306,6 +1324,8 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolSubRadAvg)
         self.assertIsNotNone(protVolSubRadAvg.outputVolume,
                              "There was a problem with Volumes subtraction radial average")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         protVolSubPdb = self.newProtocol(XmippProtVolSubtraction,
                                       vol1=protImportVol1.outputVolume,
@@ -1313,16 +1333,18 @@ def testXmippVolSub(self):
                                       pdbObj=protImportPdb.outputPdb,
                                       masks=False)
         self.launchProtocol(protVolSubPdb)
-        self.assertIsNotNone(protVolSubPdb.outputVolume,
-                             "There was a problem with Volumes subtraction pdb")
+        self.assertIsNotNone(protVolSubPdb.outputVolume, "There was a problem with Volumes subtraction pdb")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
         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")
+        self.assertIsNotNone(protVolConsensus.outputVolume, "There was a problem with Volumes consensus")
+        self.assertEqual(protVolAdj.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protVolAdj.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
 
 
 class TestXmippVolPhantom(TestXmippBase):
@@ -1350,12 +1372,6 @@ def setUpClass(cls):
         cls.vol1 = cls.dataset.getFile('volumes/volume_1_iter_002.mrc')
 
     def testXmippShiftParticlesAndVolume(self):
-
-        ERROR_SIZE_PART = "There was a problem with the size of output set of particles"
-        ERROR_DIM = "There was a problem with the dimensions of output "
-        ERROR_SR = "There was a problem with the sampling rate value of output "
-        ERROR_SHIFT = "There was a problem with output shift "
-
         protImportVol = self.newProtocol(ProtImportVolumes,
                                           objLabel='Volume',
                                           filesPath=self.vol1,
@@ -1370,7 +1386,7 @@ def testXmippShiftParticlesAndVolume(self):
                                              tiltStep=90.0)
         self.launchProtocol(protCreateGallery)
         self.assertIsNotNone(protCreateGallery.getFiles(),
-                             "There was a problem with create gallery")
+                             MSG_WRONG_GALLERY)
 
         protShiftParticles = self.newProtocol(XmippProtShiftParticles,
                                               inputParticles=protCreateGallery.outputReprojections,
@@ -1378,13 +1394,13 @@ def testXmippShiftParticlesAndVolume(self):
         self.launchProtocol(protShiftParticles)
         self.assertIsNotNone(protShiftParticles.getFiles(),
                              "There was a problem with shift particles")
-        self.assertEqual(protShiftParticles.outputParticles.getSamplingRate(), 7.08, (ERROR_SR, "particles"))
+        self.assertEqual(protShiftParticles.outputParticles.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protShiftParticles.outputParticles.getFirstItem().getDim(), (64, 64, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protShiftParticles.outputParticles.getSize(), 181, ERROR_SIZE_PART)
-        self.assertEqual(protShiftParticles.shiftX.get(), 2.0, (ERROR_SHIFT, "x"))
-        self.assertEqual(protShiftParticles.shiftY.get(), 3.0, (ERROR_SHIFT, "y"))
-        self.assertEqual(protShiftParticles.shiftZ.get(), 4.0, (ERROR_SHIFT, "z"))
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protShiftParticles.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+        self.assertEqual(protShiftParticles.shiftX.get(), 2.0, (MSG_WRONG_SHIFT, "x"))
+        self.assertEqual(protShiftParticles.shiftY.get(), 3.0, (MSG_WRONG_SHIFT, "y"))
+        self.assertEqual(protShiftParticles.shiftZ.get(), 4.0, (MSG_WRONG_SHIFT, "z"))
 
         protCreateMask = self.newProtocol(XmippProtCreateMask3D,
                                           inputVolume=protImportVol.outputVolume,
@@ -1400,13 +1416,13 @@ def testXmippShiftParticlesAndVolume(self):
         self.launchProtocol(protShiftParticlesCenterOfMass)
         self.assertIsNotNone(protShiftParticlesCenterOfMass.getFiles(),
                              "There was a problem with shift particles to center of mass")
-        self.assertEqual(protShiftParticlesCenterOfMass.outputParticles.getSamplingRate(), 7.08, (ERROR_SR, "particles"))
+        self.assertEqual(protShiftParticlesCenterOfMass.outputParticles.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protShiftParticlesCenterOfMass.outputParticles.getFirstItem().getDim(), (64, 64, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protShiftParticlesCenterOfMass.outputParticles.getSize(), 181, ERROR_SIZE_PART)
-        self.assertEqual(protShiftParticlesCenterOfMass.shiftX.get(), 32.0, (ERROR_SHIFT, "x"))
-        self.assertEqual(protShiftParticlesCenterOfMass.shiftY.get(), 32.0, (ERROR_SHIFT, "y"))
-        self.assertEqual(protShiftParticlesCenterOfMass.shiftZ.get(), 32.0, (ERROR_SHIFT, "z"))
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protShiftParticlesCenterOfMass.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+        self.assertEqual(protShiftParticlesCenterOfMass.shiftX.get(), 32.0, (MSG_WRONG_SHIFT, "x"))
+        self.assertEqual(protShiftParticlesCenterOfMass.shiftY.get(), 32.0, (MSG_WRONG_SHIFT, "y"))
+        self.assertEqual(protShiftParticlesCenterOfMass.shiftZ.get(), 32.0, (MSG_WRONG_SHIFT, "z"))
 
         protShiftVolPart = self.newProtocol(XmippProtShiftVolume,
                                             inputVol=protImportVol.outputVolume,
@@ -1414,11 +1430,11 @@ def testXmippShiftParticlesAndVolume(self):
         self.launchProtocol(protShiftVolPart)
         self.assertIsNotNone(protShiftVolPart.getFiles(),
                              "There was a problem with shift volume with particle shifts")
-        self.assertEqual(protShiftVolPart.outputVolume.getSamplingRate(), 7.08, (ERROR_SR, "volume"))
-        self.assertEqual(protShiftVolPart.outputVolume.getDim(), (64, 64, 64), (ERROR_DIM, "volume"))
-        self.assertEqual(protShiftVolPart.shiftX.get(), 2.0, (ERROR_SHIFT, "x"))
-        self.assertEqual(protShiftVolPart.shiftY.get(), 3.0, (ERROR_SHIFT, "y"))
-        self.assertEqual(protShiftVolPart.shiftZ.get(), 4.0, (ERROR_SHIFT, "z"))
+        self.assertEqual(protShiftVolPart.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protShiftVolPart.outputVolume.getDim(), (64, 64, 64), (MSG_WRONG_DIM, "volume"))
+        self.assertEqual(protShiftVolPart.shiftX.get(), 2.0, (MSG_WRONG_SHIFT, "x"))
+        self.assertEqual(protShiftVolPart.shiftY.get(), 3.0, (MSG_WRONG_SHIFT, "y"))
+        self.assertEqual(protShiftVolPart.shiftZ.get(), 4.0, (MSG_WRONG_SHIFT, "z"))
 
         protShiftVolCrop = self.newProtocol(XmippProtShiftVolume,
                                             inputVol=protImportVol.outputVolume,
@@ -1429,11 +1445,11 @@ def testXmippShiftParticlesAndVolume(self):
         self.launchProtocol(protShiftVolCrop)
         self.assertIsNotNone(protShiftVolCrop.getFiles(),
                              "There was a problem with shift crop volume")
-        self.assertEqual(protShiftVolCrop.outputVolume.getSamplingRate(), 7.08, (ERROR_SR, "volume"))
-        self.assertEqual(protShiftVolCrop.outputVolume.getDim(), (32, 32, 32), (ERROR_DIM, "volume"))
-        self.assertEqual(protShiftVolCrop.shiftX.get(), 5.0, (ERROR_SHIFT, "x"))
-        self.assertEqual(protShiftVolCrop.shiftY.get(), 5.0, (ERROR_SHIFT, "y"))
-        self.assertEqual(protShiftVolCrop.shiftZ.get(), 5.0, (ERROR_SHIFT, "z"))
+        self.assertEqual(protShiftVolCrop.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protShiftVolCrop.outputVolume.getDim(), (32, 32, 32), (MSG_WRONG_DIM, "volume"))
+        self.assertEqual(protShiftVolCrop.shiftX.get(), 5.0, (MSG_WRONG_SHIFT, "x"))
+        self.assertEqual(protShiftVolCrop.shiftY.get(), 5.0, (MSG_WRONG_SHIFT, "y"))
+        self.assertEqual(protShiftVolCrop.shiftZ.get(), 5.0, (MSG_WRONG_SHIFT, "z"))
 
         protShiftVolPad = self.newProtocol(XmippProtShiftVolume,
                                            inputVol=protImportVol.outputVolume,
@@ -1444,11 +1460,11 @@ def testXmippShiftParticlesAndVolume(self):
         self.launchProtocol(protShiftVolPad)
         self.assertIsNotNone(protShiftVolPad.getFiles(),
                              "There was a problem with shift pad volume")
-        self.assertEqual(protShiftVolCrop.outputVolume.getSamplingRate(), 7.08, (ERROR_SR, "volume"))
-        self.assertEqual(protShiftVolPad.outputVolume.getDim(), (80, 80, 80), (ERROR_DIM, "volume"))
-        self.assertEqual(protShiftVolPad.shiftX.get(), 5.0, (ERROR_SHIFT, "x"))
-        self.assertEqual(protShiftVolPad.shiftY.get(), 5.0, (ERROR_SHIFT, "y"))
-        self.assertEqual(protShiftVolPad.shiftZ.get(), 5.0, (ERROR_SHIFT, "z"))
+        self.assertEqual(protShiftVolCrop.outputVolume.getSamplingRate(), 7.08, (MSG_WRONG_SAMPLING, "volume"))
+        self.assertEqual(protShiftVolPad.outputVolume.getDim(), (80, 80, 80), (MSG_WRONG_DIM, "volume"))
+        self.assertEqual(protShiftVolPad.shiftX.get(), 5.0, (MSG_WRONG_SHIFT, "x"))
+        self.assertEqual(protShiftVolPad.shiftY.get(), 5.0, (MSG_WRONG_SHIFT, "y"))
+        self.assertEqual(protShiftVolPad.shiftZ.get(), 5.0, (MSG_WRONG_SHIFT, "z"))
 
 
 class TestXmippProjSubtraction(TestXmippBase):
@@ -1458,11 +1474,8 @@ def setUpClass(cls):
         setupTestProject(cls)
 
     def testXmippProjSub(self):
-
-        ERROR_SIZE_PART = "There was a problem with the size of output set of particles"
-        ERROR_DIM = "There was a problem with the dimensions of output "
-        ERROR_SR = "There was a problem with the sampling rate value of output "
-
+        # Create input data: phantom with two spheres and its projections (particles), phantom with one sphere
+        # (reference volume) and its mask
         protCreatePhantom2items = self.newProtocol(XmippProtPhantom,
                                                desc='80 80 80 0\nsph + 1 15 15 0 10\nsph + 5 -15 -15 0 10',
                                                sampling=1.0)
@@ -1475,13 +1488,7 @@ def testXmippProjSub(self):
                                              tiltStep=90.0)
         self.launchProtocol(protCreateGallery)
         self.assertIsNotNone(protCreateGallery.getFiles(),
-                             "There was a problem with create gallery")
-        protCreateMask = self.newProtocol(XmippProtCreateMask3D,
-                                          inputVolume=protCreatePhantom2items.outputVolume,
-                                          threshold=0.1)
-        self.launchProtocol(protCreateMask)
-        self.assertIsNotNone(protCreateMask.getFiles(),
-                             "There was a problem with the 3D mask of the 2 items phantom")
+                             MSG_WRONG_GALLERY)
         protCreatePhantom1item = self.newProtocol(XmippProtPhantom,
                                                   desc='80 80 80 0\nsph + 1 -15 -15 0 10',
                                                   sampling=1.0)
@@ -1494,42 +1501,31 @@ def testXmippProjSub(self):
         self.launchProtocol(protCreateMaskKeep)
         self.assertIsNotNone(protCreateMaskKeep.getFiles(),
                              "There was a problem with the 3D mask of the 1 item phantom")
+
+        # Subtraction of particles - reference volume with and without mask
         protSubtractProj = self.newProtocol(XmippProtSubtractProjection,
                                             particles=protCreateGallery.outputReprojections,
-                                            vol=protCreatePhantom2items.outputVolume,
-                                            maskVol=protCreateMask.outputMask,
-                                            mask=protCreateMaskKeep.outputMask,
-                                            saveFiles=True)
+                                            vol=protCreatePhantom1item.outputVolume)
         self.launchProtocol(protSubtractProj)
         self.assertIsNotNone(protSubtractProj.outputParticles,
                              "There was a problem with projection subtraction")
-        self.assertEqual(protSubtractProj.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProj.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProj.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProj.outputParticles.getSize(), 181, ERROR_SIZE_PART)
-
-        protCreateGalleryShift = self.newProtocol(XmippProtCreateGallery,
-                                                  inputVolume=protCreatePhantom2items.outputVolume,
-                                                  rotStep=15.0,
-                                                  tiltStep=90.0,
-                                                  shiftSigma=10)
-        self.launchProtocol(protCreateGalleryShift)
-        self.assertIsNotNone(protCreateGalleryShift.getFiles(),
-                             "There was a problem with create gallery shifted")
-        protSubtractProjShift = self.newProtocol(XmippProtSubtractProjection,
-                                                 particles=protCreateGalleryShift.outputReprojections,
-                                                 vol=protCreatePhantom2items.outputVolume,
-                                                 maskVol=protCreateMask.outputMask,
-                                                 mask=protCreateMaskKeep.outputMask,
-                                                 saveFiles=True)
-        self.launchProtocol(protSubtractProjShift)
-        self.assertIsNotNone(protSubtractProjShift.outputParticles,
-                             "There was a problem with projection subtraction shifted")
-        self.assertEqual(protSubtractProjShift.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
-        self.assertEqual(protSubtractProjShift.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjShift.outputParticles.getSize(), 181, ERROR_SIZE_PART)
-
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProj.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+        protSubtractProjMask = self.newProtocol(XmippProtSubtractProjection,
+                                                particles=protCreateGallery.outputReprojections,
+                                                vol=protCreatePhantom1item.outputVolume,
+                                                mask=protCreateMaskKeep.outputMask)
+        self.launchProtocol(protSubtractProjMask)
+        self.assertIsNotNone(protSubtractProjMask.outputParticles,
+                             "There was a problem with projection subtraction with mask")
+        self.assertEqual(protSubtractProjMask.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
+        self.assertEqual(protSubtractProjMask.outputParticles.getFirstItem().getDim(), (80, 80, 1),
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjMask.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+
+        # Add CTF and noise to particles (projections of the two spheres phantom) and perform the subtraction
         protSimulateCTF = self.newProtocol(XmippProtSimulateCTF,
                                            inputParticles=protCreateGallery.outputReprojections)
         self.launchProtocol(protSimulateCTF)
@@ -1537,56 +1533,48 @@ def testXmippProjSub(self):
                              "There was a problem with CTF simulation")
         protSubtractProjCTF = self.newProtocol(XmippProtSubtractProjection,
                                                particles=protSimulateCTF.outputParticles,
-                                               vol=protCreatePhantom2items.outputVolume,
-                                               maskVol=protCreateMask.outputMask,
-                                               mask=protCreateMaskKeep.outputMask,
-                                               saveFiles=True)
+                                               vol=protCreatePhantom1item.outputVolume)
         self.launchProtocol(protSubtractProjCTF)
         self.assertIsNotNone(protSubtractProjCTF.outputParticles,
                              "There was a problem with projection subtraction CTF")
-        self.assertEqual(protSubtractProjCTF.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProjCTF.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProjCTF.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjCTF.outputParticles.getSize(), 181, ERROR_SIZE_PART)
-
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjCTF.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
         protAddNoise = self.newProtocol(XmippProtAddNoiseParticles,
                                         input=protCreateGallery.outputReprojections,
-                                        gaussianStd=50.0)
+                                        gaussianStd=15.0)
         self.launchProtocol(protAddNoise)
         self.assertIsNotNone(protAddNoise.outputParticles,
                              "There was a problem with add noise protocol")
         protSubtractProjNoise = self.newProtocol(XmippProtSubtractProjection,
                                                  particles=protAddNoise.outputParticles,
-                                                 vol=protCreatePhantom2items.outputVolume,
-                                                 maskVol=protCreateMask.outputMask,
-                                                 mask=protCreateMaskKeep.outputMask,
-                                                 saveFiles=True)
+                                                 vol=protCreatePhantom1item.outputVolume)
         self.launchProtocol(protSubtractProjNoise)
         self.assertIsNotNone(protSubtractProjNoise.outputParticles,
                              "There was a problem with projection subtraction with noise")
-        self.assertEqual(protSubtractProjNoise.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProjNoise.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProjNoise.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjNoise.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoise.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
         protAddNoiseCTF = self.newProtocol(XmippProtAddNoiseParticles,
                                            input=protSimulateCTF.outputParticles,
-                                           gaussianStd=50.0)
+                                           gaussianStd=15.0)
         self.launchProtocol(protAddNoiseCTF)
         self.assertIsNotNone(protAddNoiseCTF.outputParticles,
                              "There was a problem with add noise to ctf particles protocol")
         protSubtractProjNoiseCTF = self.newProtocol(XmippProtSubtractProjection,
                                                     particles=protAddNoiseCTF.outputParticles,
-                                                    vol=protCreatePhantom2items.outputVolume,
-                                                    maskVol=protCreateMask.outputMask,
-                                                    mask=protCreateMaskKeep.outputMask,
-                                                    saveFiles=True)
+                                                    vol=protCreatePhantom1item.outputVolume)
         self.launchProtocol(protSubtractProjNoiseCTF)
         self.assertIsNotNone(protSubtractProjNoiseCTF.outputParticles,
                              "There was a problem with projection subtraction with noise and CTF")
-        self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+
+        # Create a new phantom with two *overlapping* spheres and project it to create the particles
         protCreatePhantom2Over = self.newProtocol(XmippProtPhantom,
                                                   desc='80 80 80 0\nsph + 1 5 5 0 10\nsph + 5 -5 -5 0 10',
                                                   sampling=1.0)
@@ -1594,19 +1582,12 @@ def testXmippProjSub(self):
         self.assertIsNotNone(protCreatePhantom2Over.getFiles(),
                              "There was a problem with phantom with 2 items overlap creation")
         protCreateGalleryOver = self.newProtocol(XmippProtCreateGallery,
-                                                 inputVolume=protCreatePhantom2Over.outputVolume,
-                                                 rotStep=15.0,
-                                                 tiltStep=90.0,
-                                                 shiftSigma=10)
+                                                 inputVolume=protCreatePhantom2Over.outputVolume)
         self.launchProtocol(protCreateGalleryOver)
         self.assertIsNotNone(protCreateGalleryOver.getFiles(),
                              "There was a problem with create gallery overlap")
-        protCreateMaskOver = self.newProtocol(XmippProtCreateMask3D,
-                                              inputVolume=protCreatePhantom2Over.outputVolume,
-                                              threshold=0.1)
-        self.launchProtocol(protCreateMaskOver)
-        self.assertIsNotNone(protCreateMaskOver.getFiles(),
-                             "There was a problem with the 3D mask of the 2 items overlap phantom")
+
+        # Create phantom with just one of the two overlapping spheres to use it as reference volume and mask
         protCreatePhantom1Over = self.newProtocol(XmippProtPhantom,
                                                   desc='80 80 80 0\nsph + 1 -5 -5 0 10',
                                                   sampling=1.0)
@@ -1619,19 +1600,29 @@ def testXmippProjSub(self):
         self.launchProtocol(protCreateMaskKeepOver)
         self.assertIsNotNone(protCreateMaskKeepOver.getFiles(),
                              "There was a problem with the 3D mask of the 1 item overlap phantom")
+
+        # Perform subtraction of overlapping particles with and without mask, noise and CTF
         protSubtractProjOver = self.newProtocol(XmippProtSubtractProjection,
                                                 particles=protCreateGalleryOver.outputReprojections,
-                                                vol=protCreatePhantom2Over.outputVolume,
-                                                maskVol=protCreateMaskOver.outputMask,
-                                                mask=protCreateMaskKeepOver.outputMask,
-                                                saveFiles=True)
+                                                vol=protCreatePhantom1Over.outputVolume,
+                                                mask=protCreateMaskKeepOver.outputMask)
         self.launchProtocol(protSubtractProjOver)
         self.assertIsNotNone(protSubtractProjOver.outputParticles,
                              "There was a problem with projection subtraction with overlap")
-        self.assertEqual(protSubtractProjOver.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProjOver.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProjOver.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjOver.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjOver.outputParticles.getSize(), 1647, (MSG_WRONG_SIZE, "particles"))
+        protSubtractProjOverNoMask = self.newProtocol(XmippProtSubtractProjection,
+                                                      particles=protCreateGalleryOver.outputReprojections,
+                                                      vol=protCreatePhantom1Over.outputVolume)
+        self.launchProtocol(protSubtractProjOverNoMask)
+        self.assertIsNotNone(protSubtractProjOverNoMask.outputParticles,
+                             "There was a problem with projection subtraction with overlap")
+        self.assertEqual(protSubtractProjOverNoMask.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
+        self.assertEqual(protSubtractProjOverNoMask.outputParticles.getFirstItem().getDim(), (80, 80, 1),
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjOverNoMask.outputParticles.getSize(), 1647, (MSG_WRONG_SIZE, "particles"))
         protSimulateCTFOver = self.newProtocol(XmippProtSimulateCTF,
                                                inputParticles=protCreateGalleryOver.outputReprojections)
         self.launchProtocol(protSimulateCTFOver)
@@ -1639,23 +1630,135 @@ def testXmippProjSub(self):
                              "There was a problem with overlap CTF simulation")
         protAddNoiseCTFOver = self.newProtocol(XmippProtAddNoiseParticles,
                                                input=protSimulateCTFOver.outputParticles,
-                                               gaussianStd=50.0)
+                                               gaussianStd=15.0)
         self.launchProtocol(protAddNoiseCTFOver)
         self.assertIsNotNone(protAddNoiseCTFOver.outputParticles,
                              "There was a problem with add noise to ctf overlap particles protocol")
         protSubtractProjNoiseCTFOver = self.newProtocol(XmippProtSubtractProjection,
                                                         particles=protAddNoiseCTFOver.outputParticles,
-                                                        vol=protCreatePhantom2Over.outputVolume,
-                                                        maskVol=protCreateMaskOver.outputMask,
-                                                        mask=protCreateMaskKeepOver.outputMask,
-                                                        saveFiles=True)
+                                                        vol=protCreatePhantom1Over.outputVolume)
         self.launchProtocol(protSubtractProjNoiseCTFOver)
         self.assertIsNotNone(protSubtractProjNoiseCTFOver.outputParticles,
                              "There was a problem with projection subtraction with noise and CTF overlap")
-        self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getSamplingRate(), 1.0, (ERROR_SR, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
         self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getFirstItem().getDim(), (80, 80, 1),
-                         (ERROR_DIM, "particles"))
-        self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getSize(), 1647, (MSG_WRONG_SIZE, "particles"))
+
+class TestXmippAlignVolumeAndParticles(TestXmippBase):
+
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+
+    def testXmippAlignVolumeAndParticles(self):
+        finalMatrix = np.array([[-4.32981701e-17, 1.29893286e-16, -1.00000000e+00, -0.00000000e+00],
+                           [-7.07108780e-01, 7.07104782e-01, 1.22464680e-16, -0.00000000e+00],
+                           [ 7.07104782e-01, 7.07108780e-01, 6.12323400e-17, 0.00000000e+00],
+                           [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])
+        # Create input data: phantom with two cylinders and its projections (particles)
+        protCreatePhantom1item1 = self.newProtocol(XmippProtPhantom,
+                                               desc='80 80 80 0\ncyl + 1 0 0 0 5 5 10 0 0 0',
+                                               sampling=1.0)
+        self.launchProtocol(protCreatePhantom1item1)
+        self.assertIsNotNone(protCreatePhantom1item1.getFiles(),
+                             "There was a problem with the first phantom creation")
+
+        protCreatePhantom1item2 = self.newProtocol(XmippProtPhantom,
+                                               desc='80 80 80 0\ncyl + 5 0 0 0 5 5 10 45 90 0',
+                                               sampling=1.0)
+        self.launchProtocol(protCreatePhantom1item2)
+        self.assertIsNotNone(protCreatePhantom1item2.getFiles(),
+                             "There was a problem with the second phantom creation")
+
+        protCreateGallery = self.newProtocol(XmippProtCreateGallery,
+                                             inputVolume=protCreatePhantom1item2.outputVolume,
+                                             rotStep=15.0,
+                                             tiltStep=90.0)
+        self.launchProtocol(protCreateGallery)
+        self.assertIsNotNone(protCreateGallery.getFiles(),
+                             MSG_WRONG_GALLERY)
+
+        protAlignVolumeParticles = self.newProtocol(XmippProtAlignVolumeParticles,
+                                                    inputReference=protCreatePhantom1item1.outputVolume,
+                                                    inputVolume=protCreatePhantom1item2.outputVolume,
+                                                    inputParticles=protCreateGallery.outputReprojections)
+
+        self.launchProtocol(protAlignVolumeParticles)
+        self.assertIsNotNone(protAlignVolumeParticles.outputVolume,
+                             "There was a problem with the alignment of the volume")
+        self.assertIsNotNone(protAlignVolumeParticles.outputParticles,
+                             "There was a problem with the alignment of the particles")
+        self.assertEqual(protAlignVolumeParticles.outputParticles.getSamplingRate(), 1.0, (MSG_WRONG_SAMPLING, "particles"))
+        self.assertEqual(protAlignVolumeParticles.outputParticles.getFirstItem().getDim(), (80, 80, 1),
+                         (MSG_WRONG_DIM, "particles"))
+        self.assertEqual(protAlignVolumeParticles.outputParticles.getSize(), 181, (MSG_WRONG_SIZE, "particles"))
+        self.assertEqual(protAlignVolumeParticles.outputParticles.getFirstItem().getTransform().getMatrix().all(),
+                         protCreateGallery.outputReprojections.getFirstItem().getTransform().getMatrix().all(),
+                         (MSG_WRONG_ALIGNMENT, "particles"))
+        self.assertEqual(protAlignVolumeParticles.outputParticles.getFirstItem().getTransform().getMatrix().all(),
+                         finalMatrix.all(), (MSG_WRONG_ALIGNMENT, "particles"))
+
+class TestXmippRotateVolume(TestXmippBase):
+    """This class checks if the protocol rotate volume in Xmipp works properly."""
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+
+    def testXmippAlignVolumeAndParticles(self):
+        # Create input data: phantom with one cylinder
+        protCreatePhantomRotated = self.newProtocol(XmippProtPhantom,
+                                               desc='80 80 80 0\ncyl + 5 0 0 0 5 5 10 0 90 0',
+                                               sampling=1.0)
+        self.launchProtocol(protCreatePhantomRotated)
+        self.assertIsNotNone(protCreatePhantomRotated.getFiles(),
+                             "There was a problem with the rotated phantom creation")
+
+        # First type of rotation (Align with Z)
+        protRotateVolume = self.newProtocol(XmippProtRotateVolume,
+                                            vol=protCreatePhantomRotated.outputVolume,
+                                            rotType=0,
+                                            dirParam=0)
+        self.launchProtocol(protRotateVolume)
+        self.assertIsNotNone(protRotateVolume.getFiles(),
+                             MSG_WRONG_ROTATION)
+
+        # Second type of rotation (rotate)
+        protRotateVolume2 = self.newProtocol(XmippProtRotateVolume,
+                                            vol=protCreatePhantomRotated.outputVolume,
+                                            rotType=1,
+                                            dirParam=1,
+                                            deg=90)
+        self.launchProtocol(protRotateVolume2)
+        self.assertIsNotNone(protRotateVolume2.getFiles(),
+                             MSG_WRONG_ROTATION)
+
+        # Create the referenced cylinder (without rotation)
+        # This new cylinder is the reference, meaning that it is created to check visually that
+        # the results of volume rotation protocols should look like this one
+        protCreatePhantomReference = self.newProtocol(XmippProtPhantom,
+                                                      desc='80 80 80 0\ncyl + 5 0 0 0 5 5 10 0 0 0',
+                                                      sampling=1.0)
+        self.launchProtocol(protCreatePhantomReference)
+        self.assertIsNotNone(protCreatePhantomReference.getFiles(),
+                             "There was a problem with the referenced phantom creation")
+        # First type of rotation checked (Align with Z)
+        self.assertEqual(protRotateVolume.rotType.get(), 0, "The phantom is not aligning with the Z axis")
+        self.assertEqual(protRotateVolume.dirParam.get(), 0, "The phantom is rotating in the wrong axis")
+        self.assertEqual(protRotateVolume.outputVolume.getDim(), protCreatePhantomReference.outputVolume.getDim(),
+                         (MSG_WRONG_DIM, "phantom (initially rotated)"))
+        self.assertEqual(protRotateVolume.outputVolume.getSamplingRate(),
+                         protCreatePhantomReference.outputVolume.getSamplingRate(),
+                         (MSG_WRONG_SAMPLING, "rotated phantom"))
+        # Second type of rotation checked (rotate)
+        self.assertEqual(protRotateVolume2.rotType.get(), 1, "The phantom is not rotating")
+        self.assertEqual(protRotateVolume2.dirParam.get(), 1, "The phantom is rotating in the wrong axis")
+        self.assertEqual(protRotateVolume2.deg.get(), 90, "The degree of rotation is wrong")
+        self.assertEqual(protRotateVolume2.outputVolume.getDim(), protCreatePhantomReference.outputVolume.getDim(),
+                         (MSG_WRONG_DIM, "phantom (initially rotated)"))
+        self.assertEqual(protRotateVolume2.outputVolume.getSamplingRate(),
+                         protCreatePhantomReference.outputVolume.getSamplingRate(),
+                         (MSG_WRONG_SAMPLING, "rotated phantom"))
 
 
 if __name__ == "__main__":
diff --git a/xmipp3/tests/test_protocols_zernike3d.py b/xmipp3/tests/test_protocols_zernike3d.py
index c3ba43a66..eef550c50 100644
--- a/xmipp3/tests/test_protocols_zernike3d.py
+++ b/xmipp3/tests/test_protocols_zernike3d.py
@@ -26,280 +26,280 @@
 # **************************************************************************
 
 
-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.sqrt(np.sum((clnm_gold[3] - clnm_test[3]) ** 2) / clnm_gold[3].size)
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
-                               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.sqrt(np.sum((clnm_gold[3] - clnm_test[3]) ** 2) / clnm_gold[3].size)
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_gold_def - sm_test_def) ** 2) / sm_gold_def.size)
-        self.assertAlmostEqual(diff_sm_def, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_cpu_gold_cons - sm_cpu_test_cons) ** 2) / sm_cpu_gold_cons.size)
-        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_gold_def - sm_test_def) ** 2) / sm_gold_def.size)
-        self.assertAlmostEqual(diff_sm_def, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
-                               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.sqrt(np.sum((sm_cpu_gold_cons - sm_cpu_test_cons) ** 2) / sm_cpu_gold_cons.size)
-        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.5,
-                               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.sqrt(np.sum((aa_gold_clnm - aa_test_clnm) ** 2) / aa_gold_clnm.size)
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
-                               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.sqrt(np.sum((aa_gold_clnm - aa_test_clnm) ** 2) / aa_gold_clnm.size)
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
-                               msg="Unexpected Zernike coefficients")
\ No newline at end of file
+# 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.sqrt(np.sum((clnm_gold[3] - clnm_test[3]) ** 2) / clnm_gold[3].size)
+#         self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
+#                                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.sqrt(np.sum((clnm_gold[3] - clnm_test[3]) ** 2) / clnm_gold[3].size)
+#         self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_gold_def - sm_test_def) ** 2) / sm_gold_def.size)
+#         self.assertAlmostEqual(diff_sm_def, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
+#         self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_cpu_gold_cons - sm_cpu_test_cons) ** 2) / sm_cpu_gold_cons.size)
+#         self.assertAlmostEqual(diff_sm_cons, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
+#         self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_gold_def - sm_test_def) ** 2) / sm_gold_def.size)
+#         self.assertAlmostEqual(diff_sm_def, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_cpu_gold_corr - sm_cpu_test_corr) ** 2) / sm_cpu_gold_corr.size)
+#         self.assertAlmostEqual(diff_sm_corr, 0, delta=0.5,
+#                                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.sqrt(np.sum((sm_cpu_gold_cons - sm_cpu_test_cons) ** 2) / sm_cpu_gold_cons.size)
+#         self.assertAlmostEqual(diff_sm_cons, 0, delta=0.5,
+#                                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.sqrt(np.sum((aa_gold_clnm - aa_test_clnm) ** 2) / aa_gold_clnm.size)
+#         self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
+#                                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.sqrt(np.sum((aa_gold_clnm - aa_test_clnm) ** 2) / aa_gold_clnm.size)
+#         self.assertAlmostEqual(diff_clnm, 0, delta=0.5,
+#                                msg="Unexpected Zernike coefficients")
\ No newline at end of file