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

Issue:
Revert "Co fixes"
Partially reverts I2PC/scipion-em-xmipp#471

Number of threads is used by the protocol and can affect the performance significantly.
See documentation of `cuda_reconstruct_fourier`


Discussion / hints:

Fix patch:
diff --git a/setup.py b/setup.py
index 25bd031ea..d06a58c61 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 """A setuptools based setup module.
 
 See:
-https://packaging.python.org/en/latest/distributing.html
+https://packaging.python.org/en/latest/
 https://github.com/pypa/sampleproject
 """
 
@@ -46,7 +46,7 @@
     # For a discussion on single-sourcing the version across setup.py and the
     # project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='21.9.0',  # Required
+    version='22.01.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 905252184..cf45b64e4 100644
--- a/xmipp3/__init__.py
+++ b/xmipp3/__init__.py
@@ -39,7 +39,7 @@
 
 _logo = "xmipp_logo.png"
 _references = ['delaRosaTrevin2013', 'Sorzano2013']
-_currentVersion = '3.21.09.0'
+_currentVersion = '3.22.01.0'
 
 
 class Plugin(pwem.Plugin):
diff --git a/xmipp3/convert/convert.py b/xmipp3/convert/convert.py
index 2a0b2353c..be597af8c 100644
--- a/xmipp3/convert/convert.py
+++ b/xmipp3/convert/convert.py
@@ -40,6 +40,7 @@
 import numpy as np
 
 from pyworkflow.utils import replaceBaseExt
+from pyworkflow.utils.path import cleanPath
 from pyworkflow.object import ObjectWrap, String, Float, Integer, Object
 from pwem.constants import (NO_INDEX, ALIGN_NONE, ALIGN_PROJ, ALIGN_2D,
                             ALIGN_3D)
@@ -1710,3 +1711,11 @@ def createParamPhantomFile(filename, dimX, partSetMd, phFlip=False,
     str += "_applyShift 1\n_noisePixelLevel    '0 0'\n"
     f.write(str)
     f.close()
+
+def convertToMrc(self, fnVol, Ts, deleteVol=False):
+    fnMrc = fnVol.replace(".vol",".mrc")
+    self.runJob("xmipp_image_convert","-i %s -o %s -t vol"%(fnVol, fnMrc), numberOfMpi=1)
+    self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(fnMrc, Ts), numberOfMpi=1)
+    if deleteVol:
+        cleanPath(fnVol)
+    return fnMrc
\ No newline at end of file
diff --git a/xmipp3/protocols.conf b/xmipp3/protocols.conf
index fd944e360..07946b2aa 100644
--- a/xmipp3/protocols.conf
+++ b/xmipp3/protocols.conf
@@ -78,6 +78,7 @@ Protocols SPA = [
 	{"tag": "protocol", "value": "XmippProtConvertPdb",  "text": "default"},
 	{"tag": "section", "text": "more", "openItem": "False", "children": [
 	{"tag": "protocol", "value": "XmippProtRCT", "text": "default"},
+    {"tag": "protocol", "value": "XmippProtPhantom", "text": "default"},
 	{"tag": "protocol", "value": "XmippProtShiftParticles", "text": "default"}
 	]}
 	]},
@@ -148,7 +149,9 @@ Protocols SPA = [
 			{"tag": "protocol", "value": "XmippProtImageOperateVolumes", "text": "default"},
 		    {"tag": "protocol", "value": "XmippProtVolSubtraction", "text": "default"},
 		    {"tag": "protocol", "value": "XmippProtVolAdjust", "text": "default"},
-            {"tag": "protocol", "value": "XmippProtVolConsensus", "text": "default"}
+            {"tag": "protocol", "value": "XmippProtVolConsensus", "text": "default"},
+            {"tag": "protocol", "value": "XmippProtShiftVolume", "text": "default"},
+            {"tag": "protocol", "value": "XmippProtSubtractProjection", "text": "default"}
 		]}
 	]}]
 
@@ -180,26 +183,17 @@ Random Conical Tilt = [
 	{"tag": "protocol", "value": "XmippProtRCT",                   "text": "Random Conical Tilt"}]
 
 Tomography = [
-	{"tag": "section", "text": "Imports", "icon": "bookmark.png", "children": []},
-	{"tag": "section", "text": "Movies/Micrographs", "children": []},
-	{"tag": "section", "text": "Reconstruction", "children": []},
     {"tag": "section", "text": "Particles", "children": [
-		{"tag": "protocol_group", "text": "Picking", "openItem": "False", "children": []},
-		{"tag": "protocol_group", "text": "Extract", "openItem": "False", "children": []}
+		{"tag": "protocol_group", "text": "Preprocess", "openItem": "False", "children": [
+                    {"tag": "protocol", "value": "XmippProtCreateMask3D",   "text": "default"},
+                    {"tag": "protocol", "value": "XmippProtMaskVolumes",    "text": "default"}
+                ]}
     ]},
-    {"tag": "section", "text": "Preprocessing", "children": [
-        {"tag": "protocol", "value": "XmippProtCropResizeVolumes",   "text": "default"},
-        {"tag": "protocol", "value": "XmippProtPreprocessVolumes",   "text": "default"},
-        {"tag": "protocol", "value": "XmippProtFilterVolumes",  "text": "default"},
-        {"tag": "protocol", "value": "XmippProtCreateMask3D",   "text": "default"},
-        {"tag": "protocol", "value": "XmippProtMaskVolumes",    "text": "default"},
-        {"tag": "protocol", "value": "XmippProtSubtomoProject", "text": "default"}
-        ]},
-    {"tag": "section", "text": "Subtomogram averaging", "children": []},
-    {"tag": "section", "text": "Postprocessing", "children": [
-        {"tag": "protocol", "value": "XmippProtVolSubtraction", "text": "default"},
+    {"tag": "section", "text": "Subtomogram averaging", "children": [
+    {"tag": "section", "text": "PostProcessing", "children": [
         {"tag": "protocol", "value": "XmippProtVolConsensus", "text": "default"},
-        {"tag": "protocol_group", "text": "Denoise", "openItem": "False", "children": []},
-	    {"tag": "protocol_group", "text": "Segmentation", "openItem": "False", "children": []}
-    ]}]
+        {"tag": "protocol", "value": "XmippProtSubtomoProject", "text": "default"}
+    ]}
+    ]}
+ ]
 
diff --git a/xmipp3/protocols/__init__.py b/xmipp3/protocols/__init__.py
index c753c8be9..25fa08b60 100644
--- a/xmipp3/protocols/__init__.py
+++ b/xmipp3/protocols/__init__.py
@@ -77,6 +77,7 @@
 from .protocol_particle_boxsize import XmippProtParticleBoxsize
 from .protocol_particle_pick import XmippProtParticlePicking
 from .protocol_particle_pick_pairs import XmippProtParticlePickingPairs
+from .protocol_phantom_create import XmippProtPhantom
 from .protocol_preprocess_micrographs import XmippProtPreprocessMicrographs
 from .protocol_projmatch import XmippProtProjMatch
 from .protocol_random_conical_tilt import XmippProtRCT
@@ -99,9 +100,11 @@
 from .protocol_screen_deepConsensus import XmippProtScreenDeepConsensus, XmippProtDeepConsSubSet
 from .protocol_screen_deeplearning import XmippProtScreenDeepLearning
 from .protocol_shift_particles import XmippProtShiftParticles
+from .protocol_shift_volume import XmippProtShiftVolume
 from .protocol_simulate_ctf import XmippProtSimulateCTF
 from .protocol_solid_angles import XmippProtSolidAngles
 from .protocol_split_volume import XmippProtSplitvolume
+from .protocol_subtract_projection import XmippProtSubtractProjection
 from .protocol_trigger_data import XmippProtTriggerData
 from .protocol_validate_nontilt import XmippProtValidateNonTilt
 from .protocol_validate_overfitting import XmippProtValidateOverfitting
@@ -131,4 +134,5 @@
 from .protocol_structure_map import XmippProtStructureMap
 from .protocol_pdb_deform_zernike3d import XmippProtDeformPDBZernike3D
 from .protocol_volume_adjust_sub import XmippProtVolAdjust, XmippProtVolSubtraction
-from .protocol_volume_consensus import XmippProtVolConsensus
\ No newline at end of file
+from .protocol_volume_consensus import XmippProtVolConsensus
+from .protocol_classes_2d_mapping import XmippProtCL2DMap
diff --git a/xmipp3/protocols/protocol_align_volume.py b/xmipp3/protocols/protocol_align_volume.py
index 858b1d418..44ad19269 100644
--- a/xmipp3/protocols/protocol_align_volume.py
+++ b/xmipp3/protocols/protocol_align_volume.py
@@ -191,11 +191,14 @@ def alignVolumeStep(self, refFn, inVolFn, outVolFn, maskArgs, alignArgs, volId):
             self.runJob("xmipp_volume_align", args)
     
     def createOutputStep(self):        
+        Ts = self.inputReference.get().getSamplingRate()
+
         vols = []
         idx=1
         for vol in self._iterInputVolumes():
             outVol = Volume()
-            outVol.setLocation(self._getExtraPath("vol%02d.mrc"%idx))
+            fnOutVol = self._getExtraPath("vol%02d.mrc"%idx)
+            outVol.setLocation(fnOutVol)
             outVol.setObjComment(vol.getObjComment())
             #set transformation matrix             
             fhInputTranMat = self._getExtraPath('transformation-matrix_vol%06d.txt'%idx)
@@ -205,16 +208,20 @@ def createOutputStep(self):
             transform.setMatrix(transformationMat)
             outVol.setTransform(transform)            
             vols.append(outVol)
+
+            # Set the sampling rate in the mrc header
+            self.runJob("xmipp_image_header", "-i %s --sampling_rate %f"%(fnOutVol, Ts))
+
             idx+=1
-                        
+
         if len(vols) > 1:
             volSet = self._createSetOfVolumes()
-            volSet.setSamplingRate(self.inputReference.get().getSamplingRate())
+            volSet.setSamplingRate(Ts)
             for vol in vols:
                 volSet.append(vol)
             outputArgs = {'outputVolumes': volSet}
         else:
-            vols[0].setSamplingRate(self.inputReference.get().getSamplingRate())
+            vols[0].setSamplingRate(Ts)
             outputArgs = {'outputVolume': vols[0]}
             
         self._defineOutputs(**outputArgs)
@@ -283,7 +290,7 @@ def _iterInputVolumes(self):
                 break
             itemId = item.getObjId()
             if isinstance(item, Volume):
-                item.outputName = self._getExtraPath('output_vol%06d.vol' % itemId)
+                item.outputName = self._getExtraPath('output_vol%06d.mrc' % itemId)
                 # If item is a Volume and label is empty
                 if not item.getObjLabel():
                     # Volume part of a set
@@ -294,7 +301,7 @@ def _iterInputVolumes(self):
                 yield item
             elif isinstance(item, SetOfVolumes):
                 for vol in item:
-                    vol.outputName = self._getExtraPath('output_vol%06d_%03d.vol' % (itemId, vol.getObjId()))
+                    vol.outputName = self._getExtraPath('output_vol%06d_%03d.mrc' % (itemId, vol.getObjId()))
                     # If set item label is empty
                     if not vol.getObjLabel():
                         # if set label is not empty use it
diff --git a/xmipp3/protocols/protocol_align_volume_and_particles.py b/xmipp3/protocols/protocol_align_volume_and_particles.py
index dba4b4992..188a0d4ca 100644
--- a/xmipp3/protocols/protocol_align_volume_and_particles.py
+++ b/xmipp3/protocols/protocol_align_volume_and_particles.py
@@ -28,6 +28,7 @@
 import pyworkflow.protocol.params as params
 from pyworkflow import VERSION_2_0
 from pyworkflow.protocol import STEPS_PARALLEL
+from pyworkflow.protocol.constants import LEVEL_ADVANCED
 
 from pwem.protocols import ProtAlignVolume
 from pwem.emlib.image import ImageHandler
@@ -83,6 +84,9 @@ def _defineParams(self, form):
                       label="Symmetry group",
                       help='See %s page for a description of the symmetries '
                            'accepted by Xmipp' % SYM_URL)
+        form.addParam('wrap', params.BooleanParam, default=False,
+                      label='Wrap', expertLevel=LEVEL_ADVANCED,
+                      help='Wrap the input volume when aligning to the reference')
         
         group1 = form.addGroup('Mask')
         group1.addParam('applyMask', params.BooleanParam, default=False, 
@@ -150,7 +154,7 @@ def convertStep(self):
     def alignVolumeStep(self, maskArgs):
 
         fhInputTranMat = self._getExtraPath('transformation-matrix.txt')
-        outVolFn = self._getExtraPath("inputVolumeAligned.vol")
+        outVolFn = self._getExtraPath("inputVolumeAligned.mrc")
       
         args = "--i1 %s --i2 %s --apply %s" % \
                (self.fnRefVol, self.fnInputVol, outVolFn)
@@ -159,7 +163,9 @@ def alignVolumeStep(self, maskArgs):
             args += " --frm"
         else:
             args += " --local"
-        args += " --copyGeo %s" % fhInputTranMat        
+        args += " --copyGeo %s" % fhInputTranMat
+        if not self.wrap:
+            args += ' --dontWrap'
         self.runJob("xmipp_volume_align", args)
         cleanPath(self.fnRefVol)
         cleanPath(self.fnInputVol)
@@ -192,7 +198,9 @@ def alignParticlesStep(self):
 
     def createOutputStep(self):   
 
-        outVolFn = self._getExtraPath("inputVolumeAligned.vol")
+        outVolFn = self._getExtraPath("inputVolumeAligned.mrc")
+        Ts = self.inputVolume.get().getSamplingRate()
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(outVolFn,Ts))
         outVol = Volume()
         outVol.setLocation(outVolFn)
         #set transformation matrix             
@@ -202,7 +210,7 @@ def createOutputStep(self):
         transform = Transform()
         transform.setMatrix(transformationMat)
         outVol.setTransform(transform)
-        outVol.setSamplingRate(self.inputVolume.get().getSamplingRate())
+        outVol.setSamplingRate(Ts)
 
         outputArgs = {'outputVolume': outVol}
         self._defineOutputs(**outputArgs)
diff --git a/xmipp3/protocols/protocol_analyze_local_ctf.py b/xmipp3/protocols/protocol_analyze_local_ctf.py
index 7e0884105..5ca07c749 100644
--- a/xmipp3/protocols/protocol_analyze_local_ctf.py
+++ b/xmipp3/protocols/protocol_analyze_local_ctf.py
@@ -110,7 +110,7 @@ def analyzeDefocus(self):
             mdBlock = emlib.MetaData()
             for xi, yi, deltafi, parti in zip(xbyId,ybyId,meanDefocusbyId,particleIdsbyMicId):
                 objId = mdBlock.addObject()
-                mdBlock.setValue(emlib.MDL_ITEM_ID,long(parti),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)
diff --git a/xmipp3/protocols/protocol_apply_deformation_zernike3d.py b/xmipp3/protocols/protocol_apply_deformation_zernike3d.py
index fbe719722..398fa01a4 100644
--- a/xmipp3/protocols/protocol_apply_deformation_zernike3d.py
+++ b/xmipp3/protocols/protocol_apply_deformation_zernike3d.py
@@ -1,5 +1,6 @@
 # # *************************************************
 # # * This protocol will be offered in future releases (more testing is needed)
+# # * TODO: If this protocol is added, check that sampling rate is properly set in header of mrc
 # # *************************************************
 
 # # **************************************************************************
@@ -76,7 +77,7 @@
 #             with open(file, 'w') as fid:
 #                 fid.write(' '.join(map(str, basisParams)) + "\n")
 #                 fid.write(' '.join(map(str, coeffs)) + "\n")
-#             outFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.vol' % idx
+#             outFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.mrc' % idx
 #             params = ' -i %s --clnm %s -o %s' % \
 #                      (self.inputVol.get().getFileName(), file, self._getExtraPath(outFile))
 #             self.runJob("xmipp_volume_apply_deform_sph", params)
@@ -101,7 +102,7 @@
 #
 #     def _updateClass(self, item):
 #         representative = item.getRepresentative()
-#         volumeFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.vol' \
+#         volumeFile = pwutils.removeBaseExt(self.inputVol.get().getFileName()) + '_%d_deformed.mrc' \
 #                      % (item.getObjId() + 1)
 #         volumeFile = self._getExtraPath(volumeFile)
 #         representative.setSamplingRate(self.samplingRate_Volume)
diff --git a/xmipp3/protocols/protocol_classes_2d_mapping.py b/xmipp3/protocols/protocol_classes_2d_mapping.py
new file mode 100644
index 000000000..e56b61b8a
--- /dev/null
+++ b/xmipp3/protocols/protocol_classes_2d_mapping.py
@@ -0,0 +1,399 @@
+# ******************************************************************************
+# *
+# * Authors:     David Herreros Calero (dherreros@cnb.csic.es)
+# *
+# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
+# *
+# * This program is free software; you can redistribute it and/or modify
+# * it under the terms of the GNU General Public License as published by
+# * the Free Software Foundation; either version 2 of the License, or
+# * (at your option) any later version.
+# *
+# * This program is distributed in the hope that it will be useful,
+# * but WITHOUT ANY WARRANTY; without even the implied warranty of
+# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# * GNU General Public License for more details.
+# *
+# * You should have received a copy of the GNU General Public License
+# * along with this program; if not, write to the Free Software
+# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+# * 02111-1307  USA
+# *
+# *  All comments concerning this program package may be sent to the
+# *  e-mail address 'scipion@cnb.csic.es'
+# *
+# ******************************************************************************
+
+
+import os.path
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.offsetbox import OffsetImage, AnnotationBbox
+from matplotlib.widgets import Slider, RectangleSelector, Button
+from matplotlib.cm import get_cmap, ScalarMappable
+
+from pyworkflow.object import Float
+from pyworkflow.utils.properties import Message
+from pyworkflow.gui.dialog import askYesNo
+import pyworkflow.protocol.params as param
+
+from pwem.emlib.image import ImageHandler
+import pwem.emlib.metadata as md
+from pwem.protocols import ProtAnalysis2D
+from pwem.objects import SetOfClasses2D
+
+from xmipp3.convert import xmippToLocation, writeSetOfClasses2D
+
+
+class XmippProtCL2DMap(ProtAnalysis2D):
+    """ Create a low dimensional mapping from a SetOfClasses2D with interactive selection of classes.
+    Use mouse left-click to select/deselect classes individually or mouse right-click to select/deselect
+    several classes."""
+    
+    _label = '2D classes mapping'
+    red_methods = ['PCA', 'LTSA', 'DM', 'LLTSA', 'LPP', 'kPCA', 'pPCA', 'LE', 'HLLE', 'SPE', 'NPE']
+    distances = ['Euclidean', 'Correlation']
+    CLASSES_MD = 'classes.xmd'
+    
+    def __init__(self, **args):
+        ProtAnalysis2D.__init__(self, **args)
+
+    #--------------------------- DEFINE param functions ------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('inputClasses', param.PointerParam,
+                      label="Input 2D classes",
+                      important=True, pointerClass='SetOfClasses2D',
+                      help='Select the input classes to be mapped.')
+        form.addParam('intSel', param.BooleanParam,
+                      default=True,
+                      label="Interactive class selection?")
+        form.addSection(label='Mapping')
+        form.addParam('method', param.EnumParam,
+                      choices=self.red_methods, default=0,
+                      label='Dimension reduction method')
+        form.addParam('distance', param.EnumParam,
+                      choices=self.distances, default=1,
+                      label='Distance metric to compare images')
+
+    #--------------------------- INSERT steps functions ------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('convertStep')
+        self._insertFunctionStep('computeMappingStep')
+        self._insertFunctionStep('interactiveSelStep', interactive=True)
+
+    #--------------------------- STEPS functions -------------------------------
+    def convertStep(self):
+        metadata = self._getExtraPath(self.CLASSES_MD)
+        writeSetOfClasses2D(self.inputClasses.get(),
+                            metadata,
+                            writeParticles=True)
+
+    def computeMappingStep(self):
+        metadata = self._getExtraPath(self.CLASSES_MD)
+        params = dict(metadata=metadata, method=self.red_methods[self.method.get()],
+                      metric=self.distances[self.distance.get()])
+        args = '-i classes@{metadata} -m {method} --distance {metric}'.format(**params)
+
+        self.runJob("xmipp_transform_dimred", args)
+
+    def interactiveSelStep(self):
+        selection_file = self._getExtraPath('selected_ids.txt')
+        metadata = self._getExtraPath(self.CLASSES_MD)
+        self.classes = self.inputClasses.get()
+        self.mdOut = md.MetaData(metadata)
+        if self.intSel.get():
+            img_paths = np.unique(np.asarray([rep.getFileName() for rep in self.classes.iterRepresentatives()]))
+            pos = self.mdOut.getColumnValues(md.MDL_DIMRED)
+            img_ids = self.mdOut.getColumnValues(md.MDL_REF)
+            occupancy = self.mdOut.getColumnValues(md.MDL_CLASS_COUNT)
+            occupancy = np.asarray(occupancy)
+            pos = np.vstack(pos)
+
+            # Read selected ids from previous runs (if they exist)
+            if os.path.isfile(selection_file):
+                self.selection = np.loadtxt(selection_file)
+                self.selection = [int(self.selection)] if self.selection.size == 1 else \
+                                  self.selection.astype(int).tolist()
+            else:
+                self.selection = None
+
+            view = ScatterImageMarker(pos=pos, img_paths=img_paths, ids=img_ids, occupancy=occupancy,
+                                      prevsel=self.selection)
+            view.initializePlot()
+
+            self.selection = view.selected_ids
+
+            # Save selected ids for interactive mode
+            np.savetxt(selection_file, np.asarray(self.selection))
+
+            if askYesNo(Message.TITLE_SAVE_OUTPUT, Message.LABEL_SAVE_OUTPUT, None):
+                self._createOutputStep()
+        else:
+            self.selection = [self.mdOut.getValue(md.MDL_REF, row.getObjId())
+                              for row in md.iterRows(self.mdOut)]
+            self._createOutputStep()
+
+
+    def _createOutputStep(self):
+        self._loadClassesInfo()
+        suffix = self._getOutputSuffix()
+        selected_classes = self._createSetOfClasses2D(self.classes.getImages(), suffix=suffix)
+        selected_classes.copyInfo(self.classes)
+        selected_classes.appendFromClasses(self.classes, updateClassCallback=self._updateClass,
+                                           filterClassFunc=self.isSelected)
+
+        result = {'selectedClasses2D_' + suffix: selected_classes}
+        self._defineOutputs(**result)
+        self._defineSourceRelation(self.inputClasses, selected_classes)
+
+    def isSelected(self, cls):
+        if cls.getObjId() in self.selection:
+            return True
+        else:
+            return False
+
+    def _updateClass(self, item):
+        classId = item.getObjId()
+        index, fn, _, c = self._classesInfo[classId]
+        item.setAlignment2D()
+        rep = item.getRepresentative()
+        rep.setLocation(index, fn)
+        rep.setSamplingRate(self.inputClasses.get().getSamplingRate())
+        rep._c_x = Float(c[0])
+        rep._c_y = Float(c[1])
+
+    def _loadClassesInfo(self):
+        """ Read some information about the produced 2D classes
+        from the metadata file.
+        """
+        self._classesInfo = {}  # store classes info, indexed by class id
+
+        for classNumber, row in enumerate(md.iterRows(self.mdOut)):
+            index, fn = xmippToLocation(row.getValue(md.MDL_IMAGE))
+            c = row.getValue(md.MDL_DIMRED)
+            # Store info indexed by id, we need to store the row.clone() since
+            # the same reference is used for iteration
+            self._classesInfo[classNumber + 1] = (index, fn, row.clone(), c)
+
+    #--------------------------- INFO functions --------------------------------
+    def _summary(self):
+        summary = []
+        if self.getOutputsSize() >= 1:
+            totalClasses = self.inputClasses.get().getSize()
+            for key, outClasses in self.iterOutputAttributes():
+                summary.append("*Output %s*" % key.split('_')[-1])
+                summary.append("A total of %d classes out of %d were selected"
+                               % (outClasses.getSize(), totalClasses))
+        else:
+            summary.append("No classes selected yet.")
+        return summary
+
+    def _methods(self):
+        methods = []
+        methods.append("*Dimensionality reduction method:* %s" % self.red_methods[self.method.get()])
+        methods.append("*Classes comparison metric:* %s" % self.distances[self.distance.get()])
+        return methods
+    
+    #--------------------------- UTILS functions -------------------------------
+    def _getOutputSuffix(self):
+        maxCounter = -1
+        for attrName, _ in self.iterOutputAttributes(SetOfClasses2D):
+            suffix = attrName.replace('selectedClasses2D_', '')
+            try:
+                counter = int(suffix)
+            except ValueError:
+                counter = 1  # when there is not number assume 1
+            maxCounter = max(counter, maxCounter)
+
+        return str(maxCounter+1) if maxCounter > 0 else '1'  # empty if not outputs
+
+
+class ScatterImageMarker(object):
+
+    def __init__(self, img_paths, pos, ids, occupancy=None, prevsel=None):
+        self.running = True
+        self.ids = ids
+        self.occupancy = occupancy / np.sum(occupancy) if occupancy is not None else None
+        self.selected_ids = prevsel if prevsel is not None else []
+        self.pos = pos
+        self.readImages(img_paths)
+        self.zoom = 30 / np.amax(self.images[0].shape)
+        self.pad = 0.6
+        self.artists = []
+        with plt.style.context('seaborn-darkgrid'):
+            self.fig, self.ax = plt.subplots(figsize=(10, 8))
+        self.fig.patch.set_facecolor('whitesmoke')
+        lim_x_low, lim_x_high = -1.5 * np.amax(np.abs(self.pos[:, 0])), 1.5 * np.amax(np.abs(self.pos[:, 0]))
+        lim_y_low, lim_y_high = -1.5 * np.amax(np.abs(self.pos[:, 1])), 1.5 * np.amax(np.abs(self.pos[:, 1]))
+        self.ax.set_xlim(lim_x_low, lim_x_high)
+        self.ax.set_ylim(lim_y_low, lim_y_high)
+        plt.title('Interactive class selector', fontweight="bold", fontsize=15)
+        plt.setp(self.ax.get_yticklabels(), fontweight="bold")
+        plt.setp(self.ax.get_xticklabels(), fontweight="bold")
+        plt.rcParams["font.weight"] = "bold"
+
+    def readImages(self, img_paths):
+        if len(img_paths) == 1:
+            # Numpy does not convert a list of size one to an array but to the actual object type of
+            # element inside the list. We need to consider this case as it will raise an error in Python
+            # when iterating over a non-array object type
+            self.images = np.squeeze(ImageHandler().read(img_paths[0]).getData())
+        else:
+            self.images = [np.squeeze(ImageHandler().read(img_path).getData()) for img_path in img_paths]
+
+    def imScatter(self, image, x, y, imid, edge_color=None):
+        image = OffsetImage(image, zoom=self.zoom, cmap=plt.cm.gray)
+        x, y = np.atleast_1d(x, y)
+        for x0, y0 in zip(x, y):
+            ab = AnnotationBbox(image, (x0, y0), xycoords='data', frameon=True)
+            ab.patch.set_boxstyle("Round, pad={}".format(self.pad))
+            if imid in self.selected_ids:
+                ab.patch.set_facecolor('palegreen')
+            else:
+                ab.patch.set_facecolor("lightgray")
+            ab.patch.set_alpha(0.5)
+            ab.patch.set_linewidth(2)
+            if edge_color is not None:
+                ab.patch.set_edgecolor(edge_color)
+            elif edge_color is None and imid in self.selected_ids:
+                ab.patch.set_edgecolor('darkgreen')
+            self.artists.append(self.ax.add_artist(ab))
+
+    def is_window_closed(self, event):
+        self.running = False
+
+    def updatePatch(self, ind):
+        selected_artist = self.artists[ind]
+        selected_id = self.ids[ind]
+        patch = selected_artist.patch
+        if selected_id in self.selected_ids:
+            self.selected_ids.remove(selected_id)
+            patch.set_facecolor("lightgray")
+            if self.occupancy is None:
+                patch.set_edgecolor('black')
+        else:
+            self.selected_ids.append(selected_id)
+            patch = self.artists[ind].patch
+            patch.set_facecolor('palegreen')
+            if self.occupancy is None:
+                patch.set_edgecolor('darkgreen')
+        self.fig.canvas.draw_idle()
+
+    def setPickCallback(self):
+        # Picking callback
+        def onPickImage(event):
+            ind = event.ind[0]
+            self.updatePatch(ind)
+
+        self.fig.canvas.mpl_connect('pick_event', onPickImage)
+
+    def setSliderCallback(self):
+        # Slider callback
+        slider = plt.axes([0.28, 0.01, 0.3, 0.03], facecolor='lavender') if self.occupancy is not None else \
+            plt.axes([0.355, 0.01, 0.3, 0.03], facecolor='lavender')
+        self.size_slider = Slider(ax=slider,
+                                  label='Image zoom',
+                                  valstep=0.1,
+                                  valmin=0,
+                                  valmax=2,
+                                  valinit=self.zoom,
+                                  color='springgreen')
+        def updateSize(val):
+            self.pad = val * self.pad / self.zoom
+            self.zoom = val
+            for image, artist in zip(self.images, self.artists):
+                image = OffsetImage(image, zoom=self.zoom, cmap=plt.cm.gray)
+                artist.offsetbox = image
+                artist.patch.set_boxstyle("Round, pad={}".format(self.pad))
+            self.fig.canvas.draw_idle()
+
+        self.size_slider.on_changed(updateSize)
+
+    def setRectangleSelector(self):
+        # Rectangle Selector
+        def images_select(eclick, erelease):
+            x1, y1 = eclick.xdata, eclick.ydata
+            x2, y2 = erelease.xdata, erelease.ydata
+            x_sorted, y_sorted = [x1, x2], [y1, y2]
+            x_sorted.sort()
+            y_sorted.sort()
+            x_inrange = np.logical_and(self.pos[:, 0] >= x_sorted[0], self.pos[:, 0] <= x_sorted[1])
+            y_inrange = np.logical_and(self.pos[:, 1] >= y_sorted[0], self.pos[:, 1] <= y_sorted[1])
+            ids_inrange = np.argwhere(np.logical_and(x_inrange, y_inrange)).flatten()
+            _, = [self.updatePatch(ind) for ind in ids_inrange]
+
+        def toggle_selector(event):
+            """
+            This method must be left empty. I will be completed afterwards with the object
+            to be modifed during a callback as required by Matplotlib.
+            """
+            pass
+
+        rectprops = dict(facecolor='cyan', edgecolor='gray',
+                         alpha=0.2, fill=True)
+        toggle_selector.RS = RectangleSelector(self.ax, images_select,
+                                               drawtype='box', useblit=True,
+                                               button=3,  # use only right click
+                                               minspanx=5, minspany=5,
+                                               spancoords='data',
+                                               rectprops=rectprops,
+                                               interactive=False)
+        self.fig.canvas.mpl_connect('key_press_event', toggle_selector)
+
+    def setSelectionButtons(self):
+        # Selection Buttons
+        def selectAll(event):
+            self.selected_ids = self.ids.copy()
+            for ind in self.ids:
+                patch = self.artists[ind].patch
+                patch.set_facecolor('palegreen')
+                if self.occupancy is None:
+                    patch.set_edgecolor('darkgreen')
+            self.fig.canvas.draw_idle()
+
+        def selectNone(event):
+            self.selected_ids = []
+            for ind in self.ids:
+                patch = self.artists[ind].patch
+                patch.set_facecolor("lightgray")
+                if self.occupancy is None:
+                    patch.set_edgecolor('black')
+            self.fig.canvas.draw_idle()
+
+        axprev = plt.axes([0.65, 0.015, 0.15, 0.035], facecolor='lavender')
+        axnext = plt.axes([0.81, 0.015, 0.15, 0.035], facecolor='lavender')
+        self.bnall = Button(axnext, 'Select all')
+        self.bnall.on_clicked(selectAll)
+        self.bnone = Button(axprev, 'Select none')
+        self.bnone.on_clicked(selectNone)
+
+    def initializePlot(self):
+        if self.occupancy is not None:
+            cmap = get_cmap('cool')
+            # Iterate over all class images and create the matplotlib representation that will be used
+            # to render that image in the scatter plot point position (we move the resulting empty list
+            # to a junk variable)
+            _ = [self.imScatter(image, x, y, imid, cmap(occupancy))
+                 for x, y, image, occupancy, imid in zip(self.pos[:, 0], self.pos[:, 1],
+                                                         self.images, self.occupancy, self.ids)]
+            cb = self.fig.colorbar(ScalarMappable(cmap=cmap), ax=self.ax, extend='both')
+            cb.set_label("Class occupancy", fontweight="bold", labelpad=15)
+        else:
+            # Iterate over all class images and create the matplotlib representation that will be used
+            # to render that image in the scatter plot point position (we move the resulting empty list
+            # to a junk variable)
+            _ = [self.imScatter(image, x, y, imid)
+                 for x, y, image, imid in zip(self.pos[:, 0], self.pos[:, 1], self.images, self.ids)]
+        self.ax.scatter(self.pos[:, 0], self.pos[:, 1], alpha=0, picker=True)  # Probably set pickradius param?
+
+        # Callback initialization
+        self.setPickCallback()
+        self.setSliderCallback()
+        self.setRectangleSelector()
+        self.setSelectionButtons()
+
+        # Wait until interactive plot is closed
+        self.fig.canvas.mpl_connect('close_event', self.is_window_closed)
+        while self.running:
+            plt.pause(.001)
diff --git a/xmipp3/protocols/protocol_compare_reprojections.py b/xmipp3/protocols/protocol_compare_reprojections.py
index 0000c68dd..26ff60a99 100644
--- a/xmipp3/protocols/protocol_compare_reprojections.py
+++ b/xmipp3/protocols/protocol_compare_reprojections.py
@@ -72,7 +72,7 @@ def _defineParams(self, form):
                       label='Ignore CTF',
                       help='By ignoring the CTF you will create projections more similar to what a person expects, '
                            'while by using the CTF you will create projections more similar to what the microscope sees')
-        form.addParam('evaluateResiduals', BooleanParam, default=False, expertLevel=LEVEL_ADVANCED,
+        form.addParam('doEvaluateResiduals', BooleanParam, default=False, expertLevel=LEVEL_ADVANCED,
                       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')
@@ -102,7 +102,7 @@ def _insertAllSteps(self):
             anglesFn=self.imgsFn
 
         self._insertFunctionStep("produceResiduals", vol.getFileName(), anglesFn, vol.getSamplingRate())
-        if self.evaluateResiduals.get():
+        if self.doEvaluateResiduals.get():
             self._insertFunctionStep("evaluateResiduals")
         self._insertFunctionStep("createOutputStep")
 
@@ -130,12 +130,12 @@ def produceResiduals(self, fnVol, fnAngles, Ts):
         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)
-        if self.evaluateResiduals:
+        if self.doEvaluateResiduals:
             args+=" --oresiduals %s"%self._getExtraPath("residuals.stk")
         if self.ignoreCTF:
             args+=" --ignoreCTF"
         if self.optimizeGray:
-            args+="--optimizeGray --max_gray_scale 0.95 "
+            args+=" --optimizeGray --max_gray_scale 0.95 "
         self.runJob("xmipp_angular_continuous_assign2", args)
         fnNewParticles=self._getExtraPath("images.stk")
         if os.path.exists(fnNewParticles):
@@ -225,7 +225,7 @@ def _createItemMatrix(self, particle, row):
                            emlib.MDL_CONTINUOUS_Y,
                            emlib.MDL_CORRELATION_IDX, emlib.MDL_CORRELATION_MASK,
                            emlib.MDL_CORRELATION_WEIGHT, emlib.MDL_IMED)
-        if self.evaluateResiduals:
+        if self.doEvaluateResiduals:
             setXmippAttributes(particle, row,
                                emlib.MDL_ZSCORE_RESVAR, emlib.MDL_ZSCORE_RESMEAN,
                                emlib.MDL_ZSCORE_RESCOV)
@@ -241,7 +241,7 @@ def __setXmippImage(label):
         
         __setXmippImage(emlib.MDL_IMAGE)
         __setXmippImage(emlib.MDL_IMAGE_REF)
-        if self.evaluateResiduals:
+        if self.doEvaluateResiduals:
             __setXmippImage(emlib.MDL_IMAGE_RESIDUAL)
             __setXmippImage(emlib.MDL_IMAGE_COVARIANCE)
 
diff --git a/xmipp3/protocols/protocol_ctf_consensus.py b/xmipp3/protocols/protocol_ctf_consensus.py
index 8f5bc929a..fcb29bd7d 100644
--- a/xmipp3/protocols/protocol_ctf_consensus.py
+++ b/xmipp3/protocols/protocol_ctf_consensus.py
@@ -32,7 +32,7 @@
 from cmath import rect, phase
 from math import radians, degrees
 
-from pyworkflow import VERSION_2_0
+from pyworkflow import VERSION_3_0
 from pwem.objects import SetOfCTF, SetOfMicrographs
 from pyworkflow.object import Set, Integer, Pointer
 import pyworkflow.protocol.params as params
@@ -56,12 +56,12 @@ class XmippProtCTFConsensus(ProtCTFMicrographs):
     the agreement with a secondary CTF for the same set of micrographs.
     """
     _label = 'ctf consensus'
-    _lastUpdateVersion = VERSION_2_0
+    _lastUpdateVersion = VERSION_3_0
 
     def __init__(self, **args):
         ProtCTFMicrographs.__init__(self, **args)
         self._freqResol = {}
-        self.stepsExecutionMode = params.STEPS_SERIAL
+        self.stepsExecutionMode = params.STEPS_PARALLEL
 
     def _defineParams(self, form):
         form.addSection(label='Input')
@@ -83,7 +83,7 @@ def _defineParams(self, form):
         line.addParam('maxDefocus', params.FloatParam,
                       default=40000, label='Max')
 
-        form.addParam('useAstigmatism', params.BooleanParam, default=True,
+        form.addParam('useAstigmatism', params.BooleanParam, default=False,
                       label='Use Astigmatism for selection',
                       help='Use this button to decide if carry out the '
                            'selection taking into account or not the '
@@ -94,6 +94,17 @@ def _defineParams(self, form):
                            'Angstroms. If the evaluated CTF has a '
                            'larger Astigmatism, it will be discarded.')
 
+        form.addParam('useAstigmatismPercentage', params.BooleanParam, default=True,
+                      label='Use Astigmatism percentage for selection',
+                      help='Use this button to decide if carry out the '
+                           'selection taking into account or not the '
+                           'astigmatism value.')
+        form.addParam('astigmatismPer', params.FloatParam, default=0.10,
+                      label='Astigmatism percentage', condition="useAstigmatismPercentage",
+                      help='Maximum value allowed for astigmatism '
+                           'percentage (|defocus_U-defocus_V|/mean_defocus). If the evaluated CTF has a '
+                           'larger Astigmatism Percentage, it will be discarded.')
+
         form.addParam('useResolution', params.BooleanParam, default=True,
                       label='Use Resolution for selection',
                       help='Use this button to decide if carry out the '
@@ -176,32 +187,31 @@ def _defineParams(self, form):
                            'If *No*, only the primary metadata (plus consensus '
                            'scores) will be in the resulting CTF.')
 
-        form.addParallelSection(threads=0, mpi=0)
+        form.addParallelSection(threads=4, mpi=1)
 
 # --------------------------- INSERT steps functions -------------------------
     def _insertAllSteps(self):
-        self.finished = False
-        self.insertedDict = {}
-        self.processedDict = []
-        self.outputDict = []
-        self.allCtf1 = []
-        self.allCtf2 = []
-        self.initializeRejDict()
-        self.setSecondaryAttributes()
-
-        self.ctfFn1 = self.inputCTF.get().getFileName()
+        self.initializeParams()
         if self.calculateConsensus:
             self.ctfFn2 = self.inputCTF2.get().getFileName()
-            ctfSteps = self._checkNewInput()
-        else:
-            ctfSteps = self._insertNewSelectionSteps(self.insertedDict,
-                                                     self.inputCTF.get())
+            self.allCtf2 = {}
+
+        ctfSteps = self._checkNewInput()
         self._insertFunctionStep('createOutputStep',
                                  prerequisites=ctfSteps, wait=True)
 
     def createOutputStep(self):
         pass
 
+    def initializeParams(self):
+        self.finished = False
+        self.insertedDict = {}
+        self.initializeRejDict()
+        self.setSecondaryAttributes()
+        self.ctfFn1 = self.inputCTF.get().getFileName()
+        self.allCtf1 = {}
+        pwutils.makePath(self._getExtraPath('DONE'))
+
     def _getFirstJoinStepName(self):
         # This function will be used for streaming, to check which is
         # the first function that need to wait for all ctfs
@@ -215,38 +225,43 @@ def _getFirstJoinStep(self):
                 return s
         return None
 
-    def _insertNewCtfsSteps(self, setOfCtf1, setOfCtf2, insDict):
+
+    def _insertNewCtfsSteps(self, newIDs1, newIDs2, insertedDict):
         deps = []
-        ctf1Dict = setOfCtf1.getObjDict(includeBasic=True)
-        ctf2Dict = setOfCtf2.getObjDict(includeBasic=True)
-        discrepId = self._insertFunctionStep("ctfDiscrepancyStep",
-                                             ctf1Dict, ctf2Dict,
-                                             prerequisites=[])
-        deps.append(discrepId)
-        if len(setOfCtf1) > len(setOfCtf2):
-            setOfCtf = setOfCtf2
-        else:
-            setOfCtf = setOfCtf1
-        for ctf in setOfCtf:
-            ctfId = ctf.getObjId()
-            if ctfId not in insDict:
-                stepId = self._insertFunctionStep('selectCtfStep', ctfId,
-                                                  prerequisites=[discrepId])
+
+        newIDs = list(set(newIDs1).intersection(set(newIDs2)))
+        md1 = emlib.MetaData()
+        md2 = emlib.MetaData()
+
+        for ctfID in newIDs:
+            if ctfID not in insertedDict:
+                ctf1 = self.allCtf1.get(ctfID)
+                ctf2 = self.allCtf2.get(ctfID)
+                try:
+                    self._ctfToMd(ctf1, md1)
+                    self._ctfToMd(ctf2, md2)
+                    self._freqResol[ctfID] = emlib.errorMaxFreqCTFs2D(md1, md2)
+                except TypeError as exc:
+                    print("Error reading ctf for id:%s. %s" % (ctfID, exc))
+                    self._freqResol[ctfID] = 9999
+
+                stepId = self._insertFunctionStep('selectCtfStep', ctfID,
+                                                  prerequisites=[])
                 deps.append(stepId)
-                insDict[ctfId] = stepId
+                insertedDict[ctfID] = stepId
 
         return deps
 
-    def _insertNewSelectionSteps(self, insertedDict, inputCtfs):
+    def _insertNewSelectionSteps(self, insertedDict, newIDs):
         deps = []
         # For each ctf insert the step to process it
-        for ctf in inputCtfs:
-            ctfId = ctf.getObjId()
-            if ctfId not in insertedDict:
-                stepId = self._insertFunctionStep('selectCtfStep', ctfId,
+        for ctfID in newIDs:
+
+            if ctfID not in insertedDict:
+                stepId = self._insertFunctionStep('selectCtfStep', ctfID,
                                                   prerequisites=[])
                 deps.append(stepId)
-                insertedDict[ctfId] = stepId
+                insertedDict[ctfID] = stepId
         return deps
 
     def _stepsCheck(self):
@@ -259,47 +274,39 @@ def _checkNewInput(self):
             self.lastCheck = getattr(self, 'lastCheck', datetime.now())
             mTime = max(datetime.fromtimestamp(os.path.getmtime(self.ctfFn1)),
                         datetime.fromtimestamp(os.path.getmtime(self.ctfFn2)))
+            self.debug('Last check: %s, modification: %s'
+                       % (pwutils.prettyTime(self.lastCheck),
+                          pwutils.prettyTime(mTime)))
             # If the input movies.sqlite have not changed since our last check,
             # it does not make sense to check for new input data
-            if self.lastCheck > mTime and hasattr(self, 'SetOfCtf1'):
+            if self.lastCheck > mTime and (hasattr(self, 'outputCTF') or hasattr(self, "outputCTFDiscarded")):
                 return None
+
             ctfsSet1 = self._loadInputCtfSet(self.ctfFn1)
             ctfsSet2 = self._loadInputCtfSet(self.ctfFn2)
-            if len(self.allCtf1) > 0:
-                newCtf1 = [ctf.clone() for ctf in
-                           ctfsSet1.iterItems(orderBy='creation',
-                                              where='creation>"' + str(
-                                                  self.checkCtf1) + '"')]
-            else:
-                newCtf1 = [ctf.clone() for ctf in ctfsSet1]
-            self.allCtf1 = self.allCtf1 + newCtf1
-            if len(newCtf1) > 0:
-                for ctf in ctfsSet1.iterItems(orderBy='creation',
-                                              direction='DESC'):
-                    self.checkCtf1 = ctf.getObjCreation()
-                    break
-            if len(self.allCtf2) > 0:
-                newCtf2 = [ctf.clone() for ctf in
-                           ctfsSet2.iterItems(orderBy='creation',
-                                              where='creation>"' + str(
-                                                  self.checkCtf2) + '"')]
-            else:
-                newCtf2 = [ctf.clone() for ctf in ctfsSet2]
-            self.allCtf2 = self.allCtf2 + newCtf2
-            if len(newCtf2) > 0:
-                for ctf in ctfsSet2.iterItems(orderBy='creation',
-                                              direction='DESC'):
-                    self.checkCtf2 = ctf.getObjCreation()
-                    break
+
+            ctfDict1 = {ctf.getObjId(): ctf.clone() for ctf
+                        in ctfsSet1.iterItems()}
+
+            ctfDict2 = {ctf.getObjId(): ctf.clone() for ctf
+                        in ctfsSet2.iterItems()}
+
+            newIds1 = [idCTF1 for idCTF1 in ctfDict1.keys() if idCTF1 not in self.insertedDict]
+            self.allCtf1.update(ctfDict1)
+
+            newIds2 = [idCTF2 for idCTF2 in ctfDict2.keys() if idCTF2 not in self.insertedDict]
+            self.allCtf2.update(ctfDict2)
+
             self.lastCheck = datetime.now()
             self.isStreamClosed = ctfsSet1.isStreamClosed() and \
                                   ctfsSet2.isStreamClosed()
             ctfsSet1.close()
             ctfsSet2.close()
+
             outputStep = self._getFirstJoinStep()
-            if len(set(self.allCtf1)) > len(set(self.processedDict)) and \
-               len(set(self.allCtf2)) > len(set(self.processedDict)):
-                fDeps = self._insertNewCtfsSteps(ctfsSet1, ctfsSet2,
+            if len(set(self.allCtf1)) > len(set(self.insertedDict)) and \
+               len(set(self.allCtf2)) > len(set(self.insertedDict)):
+                fDeps = self._insertNewCtfsSteps(newIds1, newIds2,
                                                  self.insertedDict)
                 if outputStep is not None:
                     outputStep.addPrerequisites(*fDeps)
@@ -311,26 +318,30 @@ def _checkNewInput(self):
             self.debug('Last check: %s, modification: %s'
                        % (pwutils.prettyTime(self.lastCheck),
                           pwutils.prettyTime(mTime)))
+            # If the input ctfs.sqlite have not changed since our last check,
+            # it does not make sense to check for new input data
+            if self.lastCheck > mTime and (hasattr(self, 'outputCTF') or hasattr(self, "outputCTFDiscarded")):
+                return None
 
             # Open input ctfs.sqlite and close it as soon as possible
             ctfSet = self._loadInputCtfSet(self.ctfFn1)
+            ctfDict = {ctf.getObjId(): ctf.clone() for ctf
+                        in ctfSet.iterItems()}
+
+            newIds = [idCTF for idCTF in ctfDict.keys() if idCTF not in self.insertedDict]
+            self.allCtf1.update(ctfDict)
+
             self.isStreamClosed = ctfSet.isStreamClosed()
-            self.allCtf1 = [m.clone() for m in ctfSet]
             ctfSet.close()
 
-            # If the input ctfs.sqlite have not changed since our last check,
-            # it does not make sense to check for new input data
-            if self.lastCheck > mTime and hasattr(self, 'allCtf1'):
-                return None
-
             self.lastCheck = now
             newCtf = any(ctf.getObjId() not in
-                         self.insertedDict for ctf in self.allCtf1)
+                         self.insertedDict for ctf in self.allCtf1.values())
             outputStep = self._getFirstJoinStep()
 
             if newCtf:
                 fDeps = self._insertNewSelectionSteps(self.insertedDict,
-                                                      self.allCtf1)
+                                                      newIds)
                 if outputStep is not None:
                     outputStep.addPrerequisites(*fDeps)
                 self.updateSteps()
@@ -351,6 +362,7 @@ def _checkNewOutput(self):
                            if ctfId not in doneListAccepted]
         newDoneDiscarded = [ctfId for ctfId in ctfListIdDiscarded
                             if ctfId not in doneListDiscarded]
+
         firstTimeAccepted = len(doneListAccepted) == 0
         firstTimeDiscarded = len(doneListDiscarded) == 0
         allDone = len(doneListAccepted) + len(doneListDiscarded) +\
@@ -436,8 +448,8 @@ def fillOutput(self, ctfSet, micSet, newDone, label):
 
                 if self.calculateConsensus:
                     ctf2 = inputCtfSet2[ctfId]
-                    setAttribute(ctf, '_consensus_resolution',
-                                 self._freqResol[ctfId])
+                    conRes = self._freqResol[ctfId]
+                    setAttribute(ctf, '_consensus_resolution', conRes)
                     setAttribute(ctf, '_ctf2_defocus_diff',
                                  max(abs(ctf.getDefocusU()-ctf2.getDefocusU()),
                                      abs(ctf.getDefocusV()-ctf2.getDefocusV())))
@@ -481,6 +493,10 @@ def fillOutput(self, ctfSet, micSet, newDone, label):
                 setAttribute(ctf, '_astigmatism',
                              abs(ctf.getDefocusU() - ctf.getDefocusV()))
 
+                # percentage _astigmatism always but after consensus if so
+                astigmatismPer = abs(ctf.getDefocusU() - ctf.getDefocusV())/(0.5 * (ctf.getDefocusU() + ctf.getDefocusV()))
+                setAttribute(ctf, '_astigmatismPercentage', astigmatismPer)
+
                 ctfSet.append(ctf)
                 micSet.append(mic)
                 self._writeCertainDoneList(ctfId, label)
@@ -534,41 +550,13 @@ def _ctfToMd(self, ctf, ctfMd):
         ctfRow = Row()
         xmipp3.convert.ctfModelToRow(ctf, ctfRow)
         xmipp3.convert.micrographToRow(ctf.getMicrograph(), ctfRow,
-                                alignType=xmipp3.convert.ALIGN_NONE)
+                                       alignType=xmipp3.convert.ALIGN_NONE)
         ctfRow.addToMd(ctfMd)
 
-    def ctfDiscrepancyStep(self, met1Dict, met2Dict):
-        # TODO must be same micrographs
-        # move to a single step, each step takes 5 sec while the function
-        # takes 0.03 sec
-        # convert to md
-        method1 = SetOfCTF()
-        method1.setAttributesFromDict(met1Dict, setBasic=True, ignoreMissing=True)
-        method2 = SetOfCTF()
-        method2.setAttributesFromDict(met2Dict, setBasic=True, ignoreMissing=True)
-        md1 = emlib.MetaData()
-        md2 = emlib.MetaData()
-
-        for ctf1 in method1:  # reference CTF
-            ctfId = ctf1.getObjId()
-            if ctfId in self.processedDict:
-                continue
-            for ctf2 in method2:
-                ctfId2 = ctf2.getObjId()
-                if ctfId2 != ctfId:
-                    continue
-                self.processedDict.append(ctfId)
-                try:
-                    self._ctfToMd(ctf1, md1)
-                    self._ctfToMd(ctf2, md2)
-                    self._freqResol[ctfId] = emlib.errorMaxFreqCTFs2D(md1, md2)
-                except TypeError as exc:
-                    print("Error reading ctf for id:%s. %s" % (ctfId, exc))
-                    self._freqResol[ctfId] = 9999
-
     def initializeRejDict(self):
         self.discDict = {'defocus': 0,
                          'astigmatism': 0,
+                         'astigmatismPer': 0,
                          'singleResolution': 0,
                          '_xmipp_ctfCritFirstZero': 0,
                          '_xmipp_ctfCritfirstZeroRatio': 0,
@@ -586,6 +574,15 @@ def selectCtfStep(self, ctfId):
         # Depending on the flags selected by the user, we set the values of
         # the params to compare with
 
+        doneFn = self._getCTFDone(ctfId)
+
+        if self.isContinued() and self._isCTFDone(ctfId):
+            self.info("Skipping CTF with ID: %s, seems to be done" % ctfId)
+            return
+
+        # Clean old finished files
+        pwutils.cleanPath(doneFn)
+
         def compareValue(ctf, label, comp, crit):
             """ Returns True if the ctf.label NOT complain the crit by comp
             """
@@ -605,19 +602,15 @@ def compareValue(ctf, label, comp, crit):
 
         minDef, maxDef = self._getDefociValues()
         maxAstig = self._getMaxAstisgmatism()
+        maxAstigPer = self._getMaxAstigmatismPer()
         minResol = self._getMinResol()
 
-        # TODO: Change this way to get the ctf.
-        ctf = self.inputCTF.get()[ctfId]
-
-        # FIXME: this is a workaround to skip errors, but it should be treat at checkNewInput
-        if ctf is None:
-            return
-
+        ctf = self.allCtf1.get(ctfId)
 
         defocusU = ctf.getDefocusU()
         defocusV = ctf.getDefocusV()
         astigm = abs(defocusU - defocusV)
+        astigmPer = abs(defocusU - defocusV)/((defocusU+defocusV)/2)
         resol = self._getCtfResol(ctf)
 
         defRangeCrit = (defocusU < minDef or defocusU > maxDef or
@@ -629,11 +622,15 @@ def compareValue(ctf, label, comp, crit):
         if astigCrit:
             self.discDict['astigmatism'] += 1
 
+        astigPer = (astigmPer > maxAstigPer)
+        if astigPer:
+            self.discDict['astigmatismPer'] += 1
+
         singleResolCrit = resol > minResol
         if singleResolCrit:
             self.discDict['singleResolution'] += 1
 
-        firstCondition = defRangeCrit or astigCrit or singleResolCrit
+        firstCondition = defRangeCrit or astigCrit or singleResolCrit or astigPer
 
         consResolCrit = False
         if self.calculateConsensus:
@@ -685,7 +682,11 @@ def compareValue(ctf, label, comp, crit):
 
         for k, v in self.discDict.items():
             setattr(self, "rejBy"+k, Integer(v))
+
         self._store()
+        # Mark this ctf as finished
+        open(doneFn, 'w').close()
+
 
     def _readDoneList(self):
         """ Read from a file the id's of the items that have been done. """
@@ -700,11 +701,19 @@ def _readDoneList(self):
     def _getAllDone(self):
         return self._getExtraPath('DONE_all.TXT')
 
-    def _writeDoneList(self, partList):
+    def _writeDoneList(self, idList):
         """ Write to a text file the items that have been done. """
         with open(self._getAllDone(), 'a') as f:
-            for part in partList:
-                f.write('%d\n' % part.getObjId())
+            for id in idList:
+                f.write('%d\n' % id)
+
+    def _isCTFDone(self, id):
+        """ A mic is done if the marker file exists. """
+        return os.path.exists(self._getCTFDone(id))
+
+    def _getCTFDone(self, id):
+        """ Return the file that is used as a flag of termination. """
+        return self._getExtraPath('DONE', 'ctf_%06d.TXT' % id)
 
     def _citations(self):
         return ['Marabini2014a']
@@ -729,7 +738,7 @@ def addDiscardedStr(label):
             obj = getattr(self, "rejBy%s" % label, Integer(0))
             number = obj.get()
             return "" if number == 0 else "  (%d discarded)" % number
-        if any([self.useDefocus, self.useAstigmatism, self.useResolution]):
+        if any([self.useDefocus, self.useAstigmatism, self.useResolution, self.useAstigmatismPercentage]):
             message.append("*General Criteria*:")
         if self.useDefocus:
             message.append(" - _Defocus_. Range: %.0f - %.0f %s"
@@ -741,6 +750,11 @@ def addDiscardedStr(label):
                            % (self.astigmatism,
                               addDiscardedStr('astigmatism')))
 
+        if self.useAstigmatismPercentage:
+            message.append(" - _AstigmatismPer_. Threshold: %.3f %s"
+                           % (self.astigmatismPer,
+                              addDiscardedStr('astigmatismPer')))
+
         if self.useResolution:
             message.append(" - _Resolution_. Threshold: %.0f %s"
                            % (self.resolution,
@@ -827,17 +841,6 @@ def _getCtfResol(self, ctf):
         else:
             return 0
 
-    # def _readDoneListDiscarded(self):
-    #     """ Read from a text file the id's of the items
-    #     that have been done. """
-    #     DiscardedFile = self._getDiscardedDone()
-    #     DiscardedList = []
-    #     # Check what items have been previously done
-    #     if os.path.exists(DiscardedFile):
-    #         with open(DiscardedFile) as f:
-    #             DiscardedList += [int(line.strip()) for line in f]
-    #     return DiscardedList
-
     def _readCertainDoneList(self, label):
         """ Read from a text file the id's of the items
         that have been done. """
@@ -855,12 +858,6 @@ def _writeCertainDoneList(self, ctfId, label):
         with open(doneFile, 'a') as f:
             f.write('%d\n' % ctfId)
 
-    # def _writeDoneListAccepted(self, ctfId):
-    #     """ Write to a text file the items that have been done. """
-    #     AcceptedFile = self._getAcceptedDone()
-    #     with open(AcceptedFile, 'a') as f:
-    #         f.write('%d\n' % ctfId)
-
     def _getCertainDone(self, label):
         return self._getExtraPath('DONE_'+label+'.TXT')
 
@@ -906,6 +903,12 @@ def _getDefociValues(self):
         else:
             return self.minDefocus.get(), self.maxDefocus.get()
 
+    def _getMaxAstigmatismPer(self):
+        if not self.useAstigmatismPercentage:
+            return 0.1
+        else:
+            return self.astigmatismPer.get()
+
     def _getMaxAstisgmatism(self):
         if not self.useAstigmatism:
             return 10000000
diff --git a/xmipp3/protocols/protocol_extract_particles.py b/xmipp3/protocols/protocol_extract_particles.py
index ceb5b6618..a11717982 100644
--- a/xmipp3/protocols/protocol_extract_particles.py
+++ b/xmipp3/protocols/protocol_extract_particles.py
@@ -32,11 +32,12 @@
 
 import pwem.emlib.metadata as md
 import pyworkflow.utils as pwutils
+from pyworkflow.object import Integer
 from pyworkflow.protocol.constants import (STEPS_PARALLEL, LEVEL_ADVANCED,
                                            STATUS_FINISHED)
 import pyworkflow.protocol.params as params
 from pwem.protocols import ProtExtractParticles
-from pwem.objects import Particle, Integer
+from pwem.objects import Particle
 
 from xmipp3.base import XmippProtocol
 from xmipp3.convert import (micrographToCTFParam, writeMicCoordinates,
@@ -101,7 +102,7 @@ def _definePreprocessParams(self, form):
                            'itself may be affected so that a higher value may '
                            'be preferable.')
 
-        form.addParam('doInvert', params.BooleanParam, default=None,
+        form.addParam('doInvert', params.BooleanParam, default=True,
                       label='Invert contrast', 
                       help='Invert the contrast if your particles are black '
                            'over a white background.  Xmipp, Spider, Relion '
@@ -109,7 +110,7 @@ def _definePreprocessParams(self, form):
                            'background. Frealign (up to v9.07) requires black '
                            'particles over a white background')
         
-        form.addParam('doFlip', params.BooleanParam, default=None,
+        form.addParam('doFlip', params.BooleanParam, default=False,
                       label='Phase flipping',
                       help='Use the information from the CTF to compensate for '
                            'phase reversals.\n'
diff --git a/xmipp3/protocols/protocol_helical_parameters.py b/xmipp3/protocols/protocol_helical_parameters.py
index 7967b99fd..fb62e2c52 100644
--- a/xmipp3/protocols/protocol_helical_parameters.py
+++ b/xmipp3/protocols/protocol_helical_parameters.py
@@ -80,8 +80,16 @@ def _insertAllSteps(self):
         self._insertFunctionStep('symmetrize')
         self._insertFunctionStep('createOutput')
         self.fnVol = getImageLocation(self.inputVolume.get())
-        self.fnVolSym=self._getPath('volume_symmetrized.vol')
+        self.fnVolSym=self._getPath('volume_symmetrized.mrc')
         [self.height,_,_]=self.inputVolume.get().getDim()
+
+    def _getFileName(self, key, **kwargs):
+        if key=="fine":
+            return self._getExtraPath('fineParams.xmd')
+        elif key=="coarse":
+            return self._getExtraPath('coarseParams.xmd')
+        else:
+            return ""
     
     #--------------------------- STEPS functions --------------------------------------------
     def copyInput(self):
@@ -102,7 +110,7 @@ def coarseSearch(self):
         self.runCoarseSearch(self.fnVolSym,self.dihedral.get(),float(self.heightFraction.get()),
                              float(self.z0.get()),float(self.zF.get()),float(self.zStep.get()),
                              float(self.rot0.get()),float(self.rotF.get()),float(self.rotStep.get()),
-                             self.numberOfThreads.get(),self._getExtraPath('coarseParams.xmd'),
+                             self.numberOfThreads.get(),self._getFileName('coarse'),
                              int(self.cylinderInnerRadius.get()),int(self.cylinderOuterRadius.get()),int(self.height),
                              self.inputVolume.get().getSamplingRate(), Cn)
 
@@ -110,8 +118,8 @@ def fineSearch(self):
         Cn = "c1"
         if self.additionalCn:
             Cn=self.Cn.get()
-        self.runFineSearch(self.fnVolSym, self.dihedral.get(), self._getExtraPath('coarseParams.xmd'),
-                           self._getExtraPath('fineParams.xmd'),
+        self.runFineSearch(self.fnVolSym, self.dihedral.get(), self._getFileName('coarse'),
+                           self._getFileName('fine'),
                            float(self.heightFraction.get()),float(self.z0.get()),float(self.zF.get()),
                            float(self.rot0.get()),float(self.rotF.get()),
                            int(self.cylinderInnerRadius.get()),int(self.cylinderOuterRadius.get()),int(self.height),
@@ -121,20 +129,21 @@ def symmetrize(self):
         Cn = "c1"
         if self.additionalCn:
             Cn=self.Cn.get()
-        self.runSymmetrize(self.fnVolSym, self.dihedral.get(), self._getExtraPath('fineParams.xmd'), self.fnVolSym,
+        self.runSymmetrize(self.fnVolSym, self.dihedral.get(), self._getFileName('fine'), self.fnVolSym,
                            float(self.heightFraction.get()),
                            self.cylinderInnerRadius.get(), self.cylinderOuterRadius.get(), self.height,
                            self.inputVolume.get().getSamplingRate(), Cn)
 
     def createOutput(self):
         volume = Volume()
-        volume.setFileName(self._getPath('volume_symmetrized.vol'))
-        #volume.setSamplingRate(self.inputVolume.get().getSamplingRate())
+        volume.setFileName(self.fnVolSym)
+        Ts = self.inputVolume.get().getSamplingRate()
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(self.fnVolSym,Ts))
         volume.copyInfo(self.inputVolume.get())
         self._defineOutputs(outputVolume=volume)
         self._defineTransformRelation(self.inputVolume, self.outputVolume)
         
-        md = MetaData(self._getExtraPath('fineParams.xmd'))
+        md = MetaData(self._getFileName('fine'))
         objId = md.firstObject()
         self._defineOutputs(deltaRot=pwobj.Float(md.getValue(MDL_ANGLE_ROT, objId)),
                             deltaZ=pwobj.Float(md.getValue(MDL_SHIFT_Z, objId)))
diff --git a/xmipp3/protocols/protocol_mltomo.py b/xmipp3/protocols/protocol_mltomo.py
index d23a86e47..817c522b6 100644
--- a/xmipp3/protocols/protocol_mltomo.py
+++ b/xmipp3/protocols/protocol_mltomo.py
@@ -81,7 +81,7 @@ def _createFilenameTemplates(self):
                   'ref': self._getExtraPath('results/mltomo_ref.xmd'),
                   'fsc_it': self.extraIter + '.fsc',
                   #'fsc': self._getExtraPath('results/mltomo.fsc'),
-                  'volume': self.extraIter + '_ref%(ref3d)06d.vol'
+                  'volume': self.extraIter + '_ref%(ref3d)06d.mrc'
                   }
         self._updateFilenamesDict(myDict)
 
@@ -432,7 +432,7 @@ def runMLTomo(self):
     def createOutput(self):
         # output files:
         #   mltomo_ref.xmd contains all info for output 3D classes
-        #   mltomo_refXXXXXX.vol output volume - 3D class
+        #   mltomo_refXXXXXX.mrc output volume - 3D class
         #   mltomo_img.xmd contains alignment metadata for all vols
         #   mltomo.fsc
         outputGlobalMdFn = self._getFileName('ref')
diff --git a/xmipp3/protocols/protocol_movie_average.py b/xmipp3/protocols/protocol_movie_average.py
index cff022adc..eea7c2641 100644
--- a/xmipp3/protocols/protocol_movie_average.py
+++ b/xmipp3/protocols/protocol_movie_average.py
@@ -151,7 +151,8 @@ def _processMovie(self, movie):
 
         self.averageMovie(movie, inputMd, outputMicFn, self.binFactor.get(),
                           roi, self.inputMovies.get().getDark(), gainFn,
-                          splineOrder=self.INTERP_MAP[self.splineOrder.get()])
+                          splineOrder=self.INTERP_MAP[self.splineOrder.get()], 
+                          outxmd="\"\"") # we're not interested in the metadata
 
         self._storeSummary(movie)
     
diff --git a/xmipp3/protocols/protocol_multireference_alignability.py b/xmipp3/protocols/protocol_multireference_alignability.py
index e3ef216ba..ad6f615c7 100644
--- a/xmipp3/protocols/protocol_multireference_alignability.py
+++ b/xmipp3/protocols/protocol_multireference_alignability.py
@@ -54,6 +54,9 @@ class XmippProtMultiRefAlignability(ProtAnalysis3D):
     """
     _label = 'multireference alignability'
 
+    INPUTARG = "-i %s"
+    OUTPUTARG = " -o %s"
+
     def __init__(self, *args, **kwargs):
         ProtAnalysis3D.__init__(self, *args, **kwargs)
 
@@ -136,6 +139,14 @@ def _defineParams(self, form):
 
         form.addParallelSection(threads=1, mpi=1)
 
+    def _getFileName(self, key, **kwargs):
+        if key=="volume":
+            return self._getExtraPath("volume.mrc")
+        elif key=="reference_particles":
+            return self._getPath("reference_particles.xmd")
+        else:
+            return ""
+
     # --------------------------- INSERT steps functions --------------------------------------------
 
     def _insertAllSteps(self):
@@ -197,8 +208,8 @@ def convertInputStep(self, particlesId):
                             self._getPath('input_particles.xmd'))
 
         if self.doWiener.get():
-            params = '  -i %s' % self._getPath('input_particles.xmd')
-            params += '  -o %s' % self._getExtraPath(
+            params = self.INPUTARG % self._getPath('input_particles.xmd')
+            params += self.OUTPUTARG % self._getExtraPath(
                 'corrected_ctf_particles.stk')
             params += '  --save_metadata_stack %s' % self._getExtraPath(
                 'corrected_ctf_particles.xmd')
@@ -221,11 +232,11 @@ def convertInputStep(self, particlesId):
         newTs, newXdim = self._getModifiedSizeAndSampling()
 
         if self.doWiener.get():
-            params =  '  -i %s' % self._getExtraPath('corrected_ctf_particles.xmd')
+            params =  self.INPUTARG % self._getExtraPath('corrected_ctf_particles.xmd')
         else :
-            params =  '  -i %s' % self._getPath('input_particles.xmd')
+            params =  self.INPUTARG % self._getPath('input_particles.xmd')
             
-        params +=  '  -o %s' % self._getExtraPath('scaled_particles.stk')
+        params +=  self.OUTPUTARG % self._getExtraPath('scaled_particles.stk')
         params +=  '  --save_metadata_stack %s' % self._getExtraPath('scaled_particles.xmd')
         params +=  '  --fourier %d' % newXdim
         
@@ -233,15 +244,14 @@ def convertInputStep(self, particlesId):
         
         from pwem.emlib.image import ImageHandler
         img = ImageHandler()
-        img.convert(self.inputVolumes.get(), self._getExtraPath("volume.vol"))
+        img.convert(self.inputVolumes.get(), self._getFileName("volume"))
         Xdim = self.inputVolumes.get().getDim()[0]
         if Xdim != newXdim:
             self.runJob("xmipp_image_resize", "-i %s --dim %d" % \
-                        (self._getExtraPath("volume.vol"),
-                         newXdim), numberOfMpi=1)
+                        (self._getFileName("volume"), newXdim), numberOfMpi=1)
 
     def _getCommonParams(self):
-        params = '  -i %s' % self._getExtraPath('scaled_particles.xmd')
+        params = self.INPUTARG % self._getExtraPath('scaled_particles.xmd')
         params += ' --sym %s' % self.symmetryGroup.get()
         params += ' --dontReconstruct'
         params += ' --useForValidation %0.3f' % (self.numOrientations.get() - 1)
@@ -249,7 +259,7 @@ def _getCommonParams(self):
         return params
 
     def _getCommonParamsRef(self):
-        params = '  -i %s' % self._getPath('reference_particles.xmd')
+        params = self.INPUTARG % self._getFileName("reference_particles")
         params += ' --sym %s' % self.symmetryGroup.get()
         params += ' --dontReconstruct'
         params += ' --useForValidation %0.3f' % (self.numOrientations.get() - 1)
@@ -285,9 +295,9 @@ def phantomProject(self):
                                self.inputParticles.get().isPhaseFlipped(),
                                self.doWiener.get()))
         f.close()
-        param = ' -i %s' % self._getExtraPath("volume.vol")
+        param = self.INPUTARG % self._getFileName("volume")
         param += ' --params %s' % self._getExtraPath('params.txt')
-        param += ' -o %s' % self._getPath('reference_particles.xmd')
+        param += self.OUTPUTARG % self._getFileName("reference_particles")
         param += ' --sampling_rate % 0.3f' % newTs
         param += ' --method fourier'
                 
@@ -297,7 +307,7 @@ def phantomProject(self):
         self.runJob('xmipp_phantom_project', 
                     param, numberOfMpi=1,numberOfThreads=1)
 
-        param = ' -i %s' % self._getPath('reference_particles.stk')
+        param = self.INPUTARG % self._getPath('reference_particles.stk')
         param += ' --mask circular %d' % R
         self.runJob('xmipp_transform_mask', param, numberOfMpi=nproc,
                     numberOfThreads=nT)
@@ -307,7 +317,7 @@ def projectionLibraryStep(self, volName, volDir, angularSampling):
         # Generate projections from this reconstruction
         nproc = self.numberOfMpi.get()
         nT = self.numberOfThreads.get()
-        volName = self._getExtraPath("volume.vol")
+        volName = self._getFileName("volume")
         makePath(volDir)
         fnGallery = (volDir + '/gallery.stk')
         params = '-i %s -o %s --sampling_rate %f --sym %s --method fourier 1 0.25 bspline --compute_neighbors --angular_distance %f --experimental_images %s --max_tilt_angle %f --min_tilt_angle %f' \
@@ -351,9 +361,9 @@ def significantStep(self, volName, volDir, anglesPath, params):
                 os.environ["CUDA_VISIBLE_DEVICES"] = GpuListAux
 
             if anglesPath == 'exp_particles.xmd':
-                params = '  -i %s' % self._getExtraPath('scaled_particles.xmd')
+                params = self.INPUTARG % self._getExtraPath('scaled_particles.xmd')
             elif anglesPath == 'ref_particles.xmd':
-                params = '  -i %s' % self._getPath('reference_particles.xmd')
+                params = self.INPUTARG % self._getFileName("reference_particles")
             params += ' --keepBestN %f' % (self.numOrientations.get() - 1)
             params += ' -r  %s' % fnGallery
             params += ' -o  %s' % self._getTmpPath(anglesPath)
@@ -367,12 +377,12 @@ def alignabilityStep(self, volName, volDir, sym):
 
         makePath(volDir)
         inputFile = self._getPath('input_particles.xmd')
-        inputFileRef = self._getPath('reference_particles.xmd')
+        inputFileRef = self._getFileName("reference_particles")
         aFile = self._getTmpPath('exp_particles.xmd')
         aFileRef = self._getTmpPath('ref_particles.xmd')
         aFileGallery = (volDir + '/gallery.doc')
 
-        params = '  -i %s' % inputFile
+        params = self.INPUTARG % inputFile
         params += ' -i2 %s' % inputFileRef
         params += '  --volume %s' % volName
         params += '  --angles_file %s' % aFile
@@ -395,7 +405,7 @@ def neighbourhoodDirectionStep(self, volName, volDir, sym):
 
         params = '  --i1 %s' % self._getPath('input_particles.xmd')
         params += ' --i2 %s' % aFileGallery
-        params += ' -o %s' % neighbours
+        params += self.OUTPUTARG % neighbours
         params += ' --dist %s' % (self.angDist.get() + 1)
         params += ' --sym %s' % sym
 
@@ -409,9 +419,9 @@ def angularAccuracyStep(self, volName, volDir, indx):
 
         neighbours = (volDir + '/neighbours.xmd')
 
-        params = ' -i %s' % volName
+        params = self.INPUTARG % volName
         params += ' --i2 %s' % neighbours
-        params += ' -o %s' % (
+        params += self.OUTPUTARG % (
                 volDir + '/pruned_particles_alignability_accuracy.xmd')
 
         self.runJob('xmipp_angular_accuracy_pca', params, numberOfMpi=nproc,
@@ -473,7 +483,7 @@ def createOutputStep(self):
         cleanPattern(self._getExtraPath("scaled_particles.*"))
         cleanPattern(self._getExtraPath("reference_particles.*"))
         cleanPattern(self._getExtraPath("corrected_ctf_particles.*"))
-        cleanPattern(self._getExtraPath("volume.vol"))
+        cleanPattern(self._getFileName("volume"))
         cleanPattern(self._getExtraPath("params.txt"))
 
     # --------------------------- INFO functions --------------------------------------------
diff --git a/xmipp3/protocols/protocol_normalize_strain.py b/xmipp3/protocols/protocol_normalize_strain.py
index 30e4877b3..32b62bbed 100644
--- a/xmipp3/protocols/protocol_normalize_strain.py
+++ b/xmipp3/protocols/protocol_normalize_strain.py
@@ -66,7 +66,7 @@ def normalize(self,what):
         for prot in self.inputRuns:
             protId=prot.get().getObjId()
             protDir=prot.get()._getPath('')
-            fnVol=os.path.join(protDir,"extra","result_%s.vol"%what)
+            fnVol=os.path.join(protDir,"extra","result_%s.mrc"%what)
             V.read(fnVol)
             _, _, minVal, maxVal = V.computeStats()
             minAll = min(minAll,minVal)
@@ -79,8 +79,8 @@ def normalize(self,what):
             fnRoot=os.path.relpath(os.path.join(protDir,"extra","result"),self._getPath(''))
             scriptFile = self._getPath('%d_result_%s_chimera.cxc'%(protId,what))
             fhCmd = open(scriptFile, 'w')
-            fhCmd.write("open %s\n" % (fnRoot+"_final.vol"))
-            fhCmd.write("open %s\n" % (fnRoot+"_%s.vol"%what))
+            fhCmd.write("open %s\n" % (fnRoot+"_final.mrc"))
+            fhCmd.write("open %s\n" % (fnRoot+"_%s.mrc"%what))
             fhCmd.write("vol #1 hide\n")
             fhCmd.write("scolor #0 volume #1 cmap rainbow cmapRange %f,%f reverseColors True\n"%(minAll,maxAll))
             fhCmd.close()
diff --git a/xmipp3/protocols/protocol_phantom_create.py b/xmipp3/protocols/protocol_phantom_create.py
new file mode 100644
index 000000000..97fa9094b
--- /dev/null
+++ b/xmipp3/protocols/protocol_phantom_create.py
@@ -0,0 +1,76 @@
+# -*- coding: utf-8 -*-
+# **************************************************************************
+# *
+# * Authors:     Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *
+# *  BCU, 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.objects.data import Volume
+from pwem.protocols import EMProtocol
+from pyworkflow.protocol.params import FloatParam, TextParam
+
+
+class XmippProtPhantom(EMProtocol):
+    """ Create phantom volume from a feature description file using xmipp_phantom_create """
+
+    _label = 'phantom volume'
+
+    # --------------------------- DEFINE param functions ------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('desc', TextParam, label='Create phantom',
+                      default='40 40 40 0\ncyl + 1 0 0 0 15 15 2 0 0 0\nsph + 1 0 0 5 2\ncyl + 1 0 0 -5 2 2 10 0 90 0\n'
+                              'sph + 1 0 -5 5 2',
+                      help="create a phantom description: x y z backgroundValue geometry(cyl, sph...) +(superimpose) "
+                           "desnsityValue origin radius height rot tilt psi. See more information in "
+                           "https://web.archive.org/web/20180813105422/http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/FileFormats#Phantom_metadata_file")
+        form.addParam('sampling', FloatParam, label='Sampling rate', default=4)
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('createPhantomsStep')
+        self._insertFunctionStep('createOutputStep')
+
+    # --------------------------- STEPS functions --------------------------------------------
+    def createPhantomsStep(self):
+        desc = self.desc.get()
+        fnDescr = self._getExtraPath("phantom.descr")
+        fhDescr = open(fnDescr, 'w')
+        fhDescr.write(desc)
+        fhDescr.close()
+        self.runJob("xmipp_phantom_create", " -i %s -o %s" % (fnDescr, self._getExtraPath("phantom.vol")))
+
+    def createOutputStep(self):
+        outputVol = Volume()
+        outputVol.setLocation(self._getExtraPath("phantom.vol"))
+        outputVol.setSamplingRate(self.sampling.get())
+        self._defineOutputs(outputVolume=outputVol)
+
+    # --------------------------- INFO functions --------------------------------------------
+    def _summary(self):
+        summary = []
+        if not hasattr(self, 'outputVolume'):
+            summary.append("Output phantom not ready yet.")
+        else:
+            summary.append("Phantoms created")
+        return summary
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py b/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
index 0abab045d..b2f0202f2 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
@@ -179,7 +179,7 @@ def convertInputStep(self):
         pass
         
     def _getNoisyOutputPath(self, fnvol):
-        fnNoisy = self._getExtraPath(removeExt(basename(fnvol)) + '_Noisy.vol')
+        fnNoisy = self._getExtraPath(removeExt(basename(fnvol)) + '_Noisy.mrc')
         return fnNoisy
 
     def _addNoisetoVolumeStep(self, kindNoise, noiseParams, vol):
@@ -188,6 +188,7 @@ def _addNoisetoVolumeStep(self, kindNoise, noiseParams, vol):
         params = " -i %s --type %s %s -o %s" % (fnvol, kindNoise, noiseParams, 
                                                     fnNoisy)
         self.runJob('xmipp_transform_add_noise', params, numberOfMpi=1)
+        self.runJob('xmipp_image_header', '-i %s --sampling_rate %f'%(fnvol, vol.getSamplingRate()), numberOfMpi=1)
 
     def addNoiseStep(self):
         kindNoise, noiseParams = self._getTypeOfNoise()
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_create_mask3d.py b/xmipp3/protocols/protocol_preprocess/protocol_create_mask3d.py
index 4d01a1064..63884f71a 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_create_mask3d.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_create_mask3d.py
@@ -192,7 +192,7 @@ def _defineParams(self, form):
 
     #--------------------------- INSERT steps functions ------------------------
     def _insertAllSteps(self):
-        self.maskFile = self._getPath('mask.vol')
+        self.maskFile = self._getPath('mask.mrc')
         
         if self.source == SOURCE_VOLUME:
             self._insertFunctionStep('createMaskFromVolumeStep')
@@ -287,10 +287,12 @@ def createOutputStep(self):
         volMask.setFileName(self.maskFile)
         
         if self.source==SOURCE_VOLUME:
-            volMask.setSamplingRate(self.inputVolume.get().getSamplingRate())
+            Ts = self.inputVolume.get().getSamplingRate()
         else:
-            volMask.setSamplingRate(self.samplingRate.get())
-        
+            Ts = self.samplingRate.get()
+        volMask.setSamplingRate(Ts)
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(self.maskFile,Ts))
+
         self._defineOutputs(outputMask=volMask)
         
         if self.source==SOURCE_VOLUME:
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_preprocess.py b/xmipp3/protocols/protocol_preprocess/protocol_preprocess.py
index ce149072d..dabe2e8b4 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_preprocess.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_preprocess.py
@@ -418,7 +418,7 @@ def _defineProcessParams(self, form):
         form.addParam('doSymmetrize', BooleanParam, default=False,
                       label="Symmetrize", 
                       help='Symmetrize the input model.')
-        form.addParam('symmetryGroup', TextParam, default='i1',
+        form.addParam('symmetryGroup', StringParam, default='i1',
                       label="Symmetry group", condition='doSymmetrize',
                       help='See [[http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/Symmetry][Symmetry]] '
                       'for a description of the symmetry groups format.'
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_process.py b/xmipp3/protocols/protocol_preprocess/protocol_process.py
index 0f347f530..b0a32033d 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_process.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_process.py
@@ -152,7 +152,9 @@ def createOutputStep(self):
             volClass = volInput.getClass()
             vol = volClass() # Create an instance with the same class of input 
             vol.copyInfo(volInput)
-            vol.setLocation(1, self.outputStk)
+            vol.setLocation(self.outputStk)
+            if self.outputStk.endswith(".mrc"):
+                self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(self.outputStk,volInput.getSamplingRate()))
             if volInput.hasOrigin():
                 vol.setOrigin(volInput.getOrigin())
             self._postprocessOutput(vol)
@@ -179,7 +181,7 @@ def _defineFilenames(self):
         """ Prepare the files to process """
         if self._isSingleInput():
             self.inputFn = getImageLocation(self.inputVolumes.get())
-            self.outputStk = self._getExtraPath("output_volume.vol")
+            self.outputStk = self._getExtraPath("output_volume.mrc")
         else:
             self.inputFn = self._getTmpPath('input_volumes.xmd')
             self.outputStk = self._getExtraPath("output_volumes.stk")
diff --git a/xmipp3/protocols/protocol_projmatch/projmatch_steps.py b/xmipp3/protocols/protocol_projmatch/projmatch_steps.py
index 8bf1ddf73..21975df73 100644
--- a/xmipp3/protocols/protocol_projmatch/projmatch_steps.py
+++ b/xmipp3/protocols/protocol_projmatch/projmatch_steps.py
@@ -38,7 +38,7 @@
 from pyworkflow.utils import getMemoryAvailable, removeExt, cleanPath, makePath, copyFile
 
 from pwem import emlib
-from xmipp3.convert import createClassesFromImages
+from xmipp3.convert import createClassesFromImages, convertToMrc
 from xmipp3.base import isMdEmpty
 
 
@@ -714,12 +714,12 @@ def runFilterVolumeStep(self, iterN, refN, constantToAddToFiltration):
                   }
         self.runJob("xmipp_transform_filter", args % params)
 
-
 def runCreateOutputStep(self):
     ''' Create standard output results_images, result_classes'''
     #creating results files
     imgSet = self.inputParticles.get()
     lastIter = self.numberOfIterations.get()
+    Ts = imgSet.getSamplingRate()
     if self.numberOfReferences != 1:
         inDocfile = self._getFileName('docfileInputAnglesIters', iter=lastIter)
         ClassFnTemplate = '%(rootDir)s/reconstruction_Ref3D_%(ref)03d.vol'
@@ -738,10 +738,11 @@ def runCreateOutputStep(self):
         classes.appendFromClasses(clsSet)
         
         volumes = self._createSetOfVolumes()
-        volumes.setSamplingRate(imgSet.getSamplingRate())
+        volumes.setSamplingRate(Ts)
         
         for refN in self.allRefs():
             volFn = self._getFileName('reconstructedFileNamesIters', iter=lastIter, ref=refN)
+            volFn = convertToMrc(self, volFn, Ts, True)
             vol = Volume()
             vol.setFileName(volFn)
             volumes.append(vol)
@@ -760,9 +761,13 @@ def runCreateOutputStep(self):
         halfMap2 = self._getFileName('reconstructedFileNamesItersSplit2',
                                      iter=lastIter, ref=1)
 
+        volFn = convertToMrc(self, volFn, Ts, True)
+        halfMap1 = convertToMrc(self, halfMap1, Ts, True)
+        halfMap2 = convertToMrc(self, halfMap2, Ts, True)
+
         vol = Volume()
         vol.setFileName(volFn)
-        vol.setSamplingRate(imgSet.getSamplingRate())
+        vol.setSamplingRate(Ts)
         vol.setHalfMaps([halfMap1, halfMap2])
         self._defineOutputs(outputVolume=vol)
         self._defineSourceRelation(self.inputParticles, vol)
diff --git a/xmipp3/protocols/protocol_random_conical_tilt.py b/xmipp3/protocols/protocol_random_conical_tilt.py
index e83915c15..dea2dfe5e 100644
--- a/xmipp3/protocols/protocol_random_conical_tilt.py
+++ b/xmipp3/protocols/protocol_random_conical_tilt.py
@@ -30,13 +30,14 @@
 from pyworkflow.object import String
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 import pyworkflow.protocol.params as params
+from pyworkflow.utils.path import cleanPath
 
 from pwem.protocols import ProtInitialVolume
 from pwem.objects import Volume, SetOfParticles
 from pwem.constants import ALIGN_2D
 from pwem import emlib
 
-from xmipp3.convert import getImageLocation, alignmentToRow
+from xmipp3.convert import getImageLocation, alignmentToRow, convertToMrc
 
 
 class XmippProtRCT(ProtInitialVolume):
@@ -297,15 +298,17 @@ def createOutputStep(self):
         self._defineSourceRelation(self.inputParticlesTiltPair, self.volumesSet)
 
     def _appendOutputVolume(self, volumeOut):
+        fnMrc = convertToMrc(self, volumeOut, self.sampling, True)
         vol = Volume()
-        vol.setFileName(volumeOut)
+        vol.setFileName(fnMrc)
         vol.setSamplingRate(self.sampling)
         self.volumesSet.append(vol)
 
         if self.doFilter.get():
             volumeFilterOut = volumeOut.replace('.vol', '_filtered.vol')
+            fnMrc = convertToMrc(self, volumeFilterOut, self.sampling, True)
             volf = Volume()
-            volf.setFileName(volumeFilterOut)
+            volf.setFileName(fnMrc)
             volf.setSamplingRate(self.sampling)
             self.volumesFilterSet.append(volf)
 
diff --git a/xmipp3/protocols/protocol_ransac.py b/xmipp3/protocols/protocol_ransac.py
index 2788f1c1a..c131f38e0 100644
--- a/xmipp3/protocols/protocol_ransac.py
+++ b/xmipp3/protocols/protocol_ransac.py
@@ -305,7 +305,10 @@ def reconstructStep(self, fnRoot):
     def resizeStep(self,fnRoot,Xdim):
         if os.path.exists(fnRoot+".vol"):
             self.runJob("xmipp_image_resize","-i %s.vol -o %s.vol --dim %d %d" %(fnRoot,fnRoot,Xdim,Xdim))
-        
+            self.runJob("xmipp_image_convert","-i %s.vol -o %s.mrc -t vol" %(fnRoot,fnRoot))
+            cleanPath("%s.vol"%fnRoot)
+            self.runJob("xmipp_image_header", "-i %s.mrc --sampling_rate %f" %\
+                        (fnRoot, self.inputSet.get().getSamplingRate()))
      
     def getCorrThreshStep(self):
         corrVector = []
@@ -436,7 +439,7 @@ def scoreFinalVolumes(self):
                 else:
                     avg=0.0
                 id=mdOut.addObject()
-                mdOut.setValue(emlib.MDL_IMAGE,fnRoot+".vol",id)
+                mdOut.setValue(emlib.MDL_IMAGE,fnRoot+".mrc",id)
                 mdOut.setValue(emlib.MDL_VOLUME_SCORE_SUM,float(sum),id)
                 mdOut.setValue(emlib.MDL_VOLUME_SCORE_SUM_TH,float(thresholdedSum),id)
                 mdOut.setValue(emlib.MDL_VOLUME_SCORE_MEAN,float(avg),id)
@@ -474,7 +477,7 @@ def createOutputStep(self):
         self._counter = 0
         readSetOfVolumes(fn, volumesSet, postprocessImageRow=self._postprocessVolume)
         
-        # Set a meanful comment
+        # Set a meaningful comment
 #         for vol in volumesSet:
 #             vol.setObjComment('ransac volume %02d' % vol.getObjId())
 #             volumesSet.update(vol)
diff --git a/xmipp3/protocols/protocol_reconstruct_fourier.py b/xmipp3/protocols/protocol_reconstruct_fourier.py
index 7b72c13f5..4424f86f3 100644
--- a/xmipp3/protocols/protocol_reconstruct_fourier.py
+++ b/xmipp3/protocols/protocol_reconstruct_fourier.py
@@ -102,7 +102,7 @@ def _createFilenameTemplates(self):
         """ Centralize how files are called for iterations and references. """
         myDict = {
             'input_xmd': self._getExtraPath('input_particles.xmd'),
-            'output_volume': self._getPath('output_volume.vol')
+            'output_volume': self._getPath('output_volume.mrc')
             }
         self._updateFilenamesDict(myDict)
 
diff --git a/xmipp3/protocols/protocol_reconstruct_heterogeneous.py b/xmipp3/protocols/protocol_reconstruct_heterogeneous.py
index 594b79f1d..a5026c15c 100644
--- a/xmipp3/protocols/protocol_reconstruct_heterogeneous.py
+++ b/xmipp3/protocols/protocol_reconstruct_heterogeneous.py
@@ -996,18 +996,14 @@ def _updateParticle(self, particle, row):
     def _updateClass(self, item):
         classId = item.getObjId()
         item.setAlignmentProj()
-        #item.setSamplingRate(self.Ts)
-        #item.getRepresentative().setFileName(
-        #    join(self.fnLastDir, "volume%02d.mrc" % classId))
+        fnMrc = join(self.fnLastDir, "volume%02d.mrc" % classId)
         if self.scaleFactor!=1.0:
             newXdim=self.inputParticles.get().getDimensions()[0]
-            self.runJob("xmipp_image_resize", "-i %s -o %s --fourier %d" % (
-                join(self.fnLastDir, "volume%02d.mrc" % classId),
-                join(self.fnLastDir, "volume%02d.mrc" % classId), newXdim),
-                        numberOfMpi=self.numberOfMpi.get() * self.numberOfThreads.get())
-        item.setSamplingRate(self.inputParticles.get().getSamplingRate())
-        item.getRepresentative().setFileName(
-            join(self.fnLastDir, "volume%02d.mrc" % classId))
+            self.runJob("xmipp_image_resize", "-i %s --fourier %d" % (fnMrc, newXdim), numberOfMpi=1)
+        Ts = self.inputParticles.get().getSamplingRate()
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(fnMrc,Ts), numberOfMpi=1)
+        item.setSamplingRate(Ts)
+        item.getRepresentative().setFileName(fnMrc)
 
     def _saveVolumesToProcess(self, volumesToProcess):
         fn = open(self._getExtraPath('volumesToProcess.txt'), 'w')
diff --git a/xmipp3/protocols/protocol_reconstruct_significant.py b/xmipp3/protocols/protocol_reconstruct_significant.py
index 54a045a61..6f5ab8d78 100644
--- a/xmipp3/protocols/protocol_reconstruct_significant.py
+++ b/xmipp3/protocols/protocol_reconstruct_significant.py
@@ -74,7 +74,7 @@ def _defineParams(self, form):
                       help='Select the input classes2D from the project.\n'
                            'It should be a SetOfClasses2D class with  class '
                            'representative')
-        form.addParam('symmetryGroup', TextParam, default='c1',
+        form.addParam('symmetryGroup', StringParam, default='c1',
                       label="Symmetry group",
                       help='See [[http://xmipp.cnb.csic.es/twiki/bin/view/Xmipp/Symmetry][Symmetry]]'
                            'for a description of the symmetry groups format.'
@@ -374,6 +374,9 @@ def significantStep(self, iterNumber, alpha):
         if not self.keepIntermediate:
             cleanPath(prevVolFn, iterDir)
 
+        if self.thereisRefVolume:
+            cleanPath(self._getExtraPath('filteredVolume.vol'))
+
     def convertInputStep(self, classesFn):
         inputSet = self.inputSet.get()
 
@@ -413,7 +416,8 @@ def convertInputStep(self, classesFn):
         # previously
         if self.thereisRefVolume:
             fnFilVol = self._getExtraPath('filteredVolume.vol')
-            copy(self.refVolume.get().getFileName(), fnFilVol)
+            self.runJob("xmipp_image_convert", "-i %s -o %s -t vol" % (self.refVolume.get().getFileName(), fnFilVol),
+                        numberOfMpi=1)
             # TsVol = self.refVolume.get().getSamplingRate()
             if self.useMaxRes:
                 if self.newXdim != Xdim:
@@ -440,18 +444,23 @@ def convertInputStep(self, classesFn):
 
     def createOutputStep(self):
         lastIter = self.getLastIteration(1)
+        Ts = self.inputSet.get().getSamplingRate()
 
         # To recover the original size of the volume if it was changed
-        volFn = self.getIterVolume(lastIter)
+        fnVol = self.getIterVolume(lastIter)
         Xdim = self.inputSet.get().getDimensions()[0]
         if self.useMaxRes and self.newXdim != Xdim:
             self.runJob('xmipp_image_resize', "-i %s --fourier %d" %
-                        (volFn, Xdim), numberOfMpi=1)
+                        (fnVol, Xdim), numberOfMpi=1)
+        fnMrc = fnVol.replace(".vol",".mrc")
+        self.runJob("xmipp_image_convert","-i %s -o %s -t vol"%(fnVol,fnMrc),numberOfMpi=1)
+        cleanPath(fnVol)
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(fnMrc,Ts),numberOfMpi=1)
 
         vol = Volume()
         vol.setObjComment('significant volume 1')
-        vol.setLocation(volFn)
-        vol.setSamplingRate(self.inputSet.get().getSamplingRate())
+        vol.setLocation(fnMrc)
+        vol.setSamplingRate(Ts)
         self._defineOutputs(outputVolume=vol)
         self._defineSourceRelation(self.inputSet, vol)
 
diff --git a/xmipp3/protocols/protocol_reconstruct_swarm.py b/xmipp3/protocols/protocol_reconstruct_swarm.py
index f9f342142..4263774e6 100644
--- a/xmipp3/protocols/protocol_reconstruct_swarm.py
+++ b/xmipp3/protocols/protocol_reconstruct_swarm.py
@@ -335,9 +335,13 @@ def createOutput(self):
         TsOrig=self.inputParticles.get().getSamplingRate()
         XdimOrig=self.inputParticles.get().getDimensions()[0]
         fnAvg = self._getExtraPath("volumeAvg.vol")
+        fnAvgMrc = self._getExtraPath("volumeAvg.mrc")
         self.runJob("xmipp_image_resize","-i %s --dim %d"%(fnAvg,XdimOrig),numberOfMpi=1)
+        self.runJob("xmipp_image_convert","-i %s -o %s -t vol"%(fnAvg,fnAvgMrc),numberOfMpi=1)
+        cleanPath(fnAvg)
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(fnAvgMrc,TsOrig),numberOfMpi=1)
         volume=Volume()
-        volume.setFileName(fnAvg)
+        volume.setFileName(fnAvgMrc)
         volume.setSamplingRate(TsOrig)
         self._defineOutputs(outputVolume=volume)
         self._defineSourceRelation(self.inputParticles.get(),volume)
@@ -348,8 +352,13 @@ def createOutput(self):
         volSet.setSamplingRate(Ts)
         for i in range(self.inputVolumes.get().getSize()):
             fnVol = self._getExtraPath("volume%03d_best.vol"%i)
+            fnMrc = self._getExtraPath("volume%03d_best.mrc"%i)
+            self.runJob("xmipp_image_convert", "-i %s -o %s" % (fnVol, fnMrc), numberOfMpi=1)
+            self.runJob("xmipp_image_header", "-i %s --sampling_rate %f" % (fnMrc, TsOrig), numberOfMpi=1)
+            cleanPath(fnVol)
+
             vol=Volume()
-            vol.setFileName(fnVol)
+            vol.setFileName(fnMrc)
             vol.setSamplingRate(Ts)
             volSet.append(vol)
         self._defineOutputs(outputVolumes=volSet)
diff --git a/xmipp3/protocols/protocol_resolution_bfactor.py b/xmipp3/protocols/protocol_resolution_bfactor.py
index 79b4271af..786bea3eb 100644
--- a/xmipp3/protocols/protocol_resolution_bfactor.py
+++ b/xmipp3/protocols/protocol_resolution_bfactor.py
@@ -55,6 +55,7 @@ class XmippProtbfactorResolution(ProtAnalysis3D):
     _lastUpdateVersion = VERSION_2_0
 
     def __init__(self, **args):
+        self.vol = ''
         ProtAnalysis3D.__init__(self, **args)
 
     # --------------------------- DEFINE param functions ----------------------
@@ -104,6 +105,11 @@ def _defineParams(self, form):
                       help='The local resolution per residue can be estimated using'
                            'the mean (by default - No) or the median (yes)')
 
+        form.addParam('centered', BooleanParam, default=True,
+                      label="is the atomic centered",
+                      help='True if the atomic model centered in midle of'
+                           ' the local resolution map')
+
         form.addParallelSection(threads=4, mpi=0)
 
     # --------------------------- INSERT steps functions --------------------------------------------
@@ -148,6 +154,8 @@ def matchingBfactorLocalResolution(self):
         params += ' --sampling %f' % self.localResolutionMap.get().getSamplingRate()
         if self.medianEstimation.get():
             params += ' --useMedian '
+        if self.centered.get():
+            params += ' --centered '
         params += ' -o %s' % self._getExtraPath()
 
         self.runJob('xmipp_resolution_pdb_bfactor', params)
diff --git a/xmipp3/protocols/protocol_resolution_fso.py b/xmipp3/protocols/protocol_resolution_fso.py
index 06f1ae04c..d2972beb7 100644
--- a/xmipp3/protocols/protocol_resolution_fso.py
+++ b/xmipp3/protocols/protocol_resolution_fso.py
@@ -97,6 +97,12 @@ def _defineParams(self, form):
                            ' 3DFSC as anisotropic filter to obtain a directional'
                            'filtered map.')
 
+#        form.addParam('filterHalfMaps', BooleanParam, default=True,
+#                      condition="estimate3DFSC",
+#                      label="Anisotropic filtering of Half Maps",
+#                      help='Use the 3DFSCD map as anisotropic filter to directionally'
+#                           'filter the half maps')
+
         form.addParam('threshold', FloatParam, expertLevel=LEVEL_ADVANCED,
                       default=0.143,
                       label="FSC Threshold",
@@ -158,19 +164,6 @@ def convertInputStep(self):
             if (extMask == '.mrc') or (extMask == '.map'):
                 self.maskFn = self.maskFn + ':mrc'
 
-        """
-        if self.halfVolumesFile:
-            self.vol1Fn, self.vol2Fn = self.inputHalves.get().getHalfMaps().split(',')
-        else:
-            self.vol1Fn = self.mrc_convert(self.half1.get().getFileName(),
-                                  self._getTmpPath('half1.mrc'))
-            self.vol2Fn = self.mrc_convert(self.half2.get().getFileName(),
-                                  self._getTmpPath('half2.mrc'))
-        if (self.mask.hasValue()):
-            self.maskFn = self.mrc_convert(self.mask.get().getFileName(),
-                                  self._getExtraPath('mask.mrc'))
-        """
-
 
     def FSOestimationStep(self):
         import os
@@ -193,6 +186,10 @@ def FSOestimationStep(self):
 
         if self.estimate3DFSC.get():
             params += ' --threedfsc_filter'
+            #if self.filterHalfMaps.get():
+            #    params += ' --threedfsc_filter_halves'
+
+
 
         params += ' --threshold %s' % self.threshold.get()
         params += ' --threads %s' % self.numberOfThreads.get()
@@ -204,24 +201,31 @@ def createOutputStep(self):
 
         if self.estimate3DFSC:
             volume = Volume()
-            volume.setFileName(self._getExtraPath("3dFSC.mrc"))
             volume.setFileName(self._getExtraPath("filteredMap.mrc"))
+
+            #if self.filterHalfMaps:
+            #    half1filt = Volume()
+            #    half2filt = Volume()
+            #    fnhalfMap1 = self._getExtraPath("dirfiltered_half1.mrc")
+            #    fnhalfMap2 = self._getExtraPath("dirfiltered_half2.mrc")
+            #    half1filt.setFileName(fnhalfMap1)
+            #    half2filt.setFileName(fnhalfMap2)
+
             if self.halfVolumesFile:
                 volume.setSamplingRate(self.inputHalves.get().getSamplingRate())
-                self._defineOutputs(fsc3D=volume)
+                #if self.filterHalfMaps:
+                #    volume.setHalfMaps([fnhalfMap1, fnhalfMap2])
                 self._defineOutputs(directionalFilteredMap=volume)
                 self._defineSourceRelation(self.inputHalves, volume)
             else:
                 volume.setSamplingRate(self.half1.get().getSamplingRate())
-                self._defineOutputs(fsc3D=volume)
+                #if self.filterHalfMaps:
+                #    volume.setHalfMaps([fnhalfMap1, fnhalfMap2])
                 self._defineOutputs(directionalFilteredMap=volume)
                 self._defineSourceRelation(self.half1, volume)
 
 
 
-
-
-
     # --------------------------- INFO functions ------------------------------
 
     def _methods(self):
diff --git a/xmipp3/protocols/protocol_rotational_symmetry.py b/xmipp3/protocols/protocol_rotational_symmetry.py
index b97c06e78..34a2b9ceb 100644
--- a/xmipp3/protocols/protocol_rotational_symmetry.py
+++ b/xmipp3/protocols/protocol_rotational_symmetry.py
@@ -58,7 +58,7 @@ def _defineParams(self, form):
                       label='Symmetry order',
                       help="3 for a three-fold symmetry axis, "
                            "4 for a four-fold symmetry axis, ...")
-        form.addParam('searchMode', EnumParam,
+        form.addParam('searchMode', EnumParam,label='Search mode',
                       choices=['Global','Local','Global+Local'],
                       default=self.GLOBAL_LOCAL_SEARCH)
 
@@ -102,7 +102,7 @@ def _insertAllSteps(self):
         self._insertFunctionStep('symmetrize')
         self._insertFunctionStep('createOutput')
         self.fnVol = getImageLocation(self.inputVolume.get())
-        self.fnVolSym=self._getPath('volume_symmetrized.vol')
+        self.fnVolSym=self._getPath('volume_symmetrized.mrc')
         [self.height,_,_]=self.inputVolume.get().getDim()
     
     #--------------------------- STEPS functions -------------------------------
@@ -150,6 +150,9 @@ def symmetrize(self):
                     (self.fnVolSym,self.symOrder.get()))
 
     def createOutput(self):
+        Ts = self.inputVolume.get().getSamplingRate()
+        self.runJob("xmipp_image_header","-i %s --sampling_rate %f" %(self.fnVolSym,Ts))
+
         volume = Volume()
         volume.setFileName(self.fnVolSym)
         volume.copyInfo(self.inputVolume.get())
diff --git a/xmipp3/protocols/protocol_shift_particles.py b/xmipp3/protocols/protocol_shift_particles.py
index 249e8f4e8..162834ddb 100644
--- a/xmipp3/protocols/protocol_shift_particles.py
+++ b/xmipp3/protocols/protocol_shift_particles.py
@@ -26,13 +26,15 @@
 # *
 # **************************************************************************
 
-from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, EnumParam
+from pwem.emlib.image import ImageHandler as ih
+from pwem.objects import Volume
+from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, EnumParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 import pyworkflow.object as pwobj
 from pwem import ALIGN_3D, ALIGN_2D
 import pwem.emlib.metadata as md
 from pwem.protocols import EMProtocol
-from xmipp3.convert import xmippToLocation, writeSetOfParticles
+from xmipp3.convert import xmippToLocation, writeSetOfParticles, readSetOfParticles
 
 
 class XmippProtShiftParticles(EMProtocol):
@@ -46,13 +48,19 @@ def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputParticles', PointerParam, pointerClass='SetOfParticles', label="Particles",
                       help='Select the SetOfParticles with transformation matrix to be shifted.')
+        form.addParam('option', BooleanParam, label='Select position in volume?', default='True',
+                      help='Select the position where the particles will be shifted in a volume displayed in a wizard.')
         form.addParam('inputVol', PointerParam, pointerClass='Volume', label="Volume", allowsNull=True,
-                      help='Volume to select the point (by clicking in the wizard for selecting the new center) that '
-                           'will be the new center of the particles.')
-        form.addParam('x', IntParam, label="x", help='Use the wizard to select by clicking in the volume the new '
-                                                     'center for the shifted particles')
-        form.addParam('y', IntParam, label="y")
-        form.addParam('z', IntParam, label="z")
+                      condition='option', help='Volume to select the point (by clicking in the wizard for selecting the'
+                                               ' new center) that will be the new center of the particles.')
+        form.addParam('xin', FloatParam, label="x", condition='option',
+                      help='Use the wizard to select the new center for the shifted particles by shift+click on the '
+                           'blue point a drag it to the desired location while pressing shift.')
+        form.addParam('yin', FloatParam, label="y", condition='option')
+        form.addParam('zin', FloatParam, label="z", condition='option')
+        form.addParam('inputMask', PointerParam, pointerClass='VolumeMask', label="Volume mask", allowsNull=True,
+                      condition='not option', help='3D mask to compute the center of mass, the particles will be '
+                                                   'shifted to the computed center of mass')
         form.addParam('boxSizeBool', BooleanParam, label='Use original box size for the shifted particles?',
                       default='True', help='Use input particles box size for the shifted particles.')
         form.addParam('boxSize', IntParam, label='Final box size', condition='not boxSizeBool',
@@ -60,7 +68,7 @@ def _defineParams(self, form):
         form.addParam('inv', BooleanParam, label='Inverse', expertLevel=LEVEL_ADVANCED, default='True',
                       help='Use inverse transformation matrix')
         form.addParam('interp', EnumParam, default=0, choices=['Linear', 'Spline'], expertLevel=LEVEL_ADVANCED,
-                      display=EnumParam.DISPLAY_HLIST , label='Interpolation',
+                      display=EnumParam.DISPLAY_HLIST, label='Interpolation',
                       help='Linear: Use bilinear/trilinear interpolation\nSpline: Use spline interpolation')
 
     # --------------------------- INSERT steps functions --------------------------------------------
@@ -77,14 +85,31 @@ def convertStep(self):
 
     def shiftStep(self):
         """call xmipp program to shift the particles"""
-        program = "xmipp_transform_geometry"
         centermd = self._getExtraPath("center_particles.xmd")
+        args = '-i "%s" -o "%s" ' % (self._getExtraPath("input_particles.xmd"), centermd)
+        if self.option:
+            self.x = self.xin.get()
+            self.y = self.yin.get()
+            self.z = self.zin.get()
+        else:
+            fnvol = self.inputMask.get().getFileName()
+            if fnvol.endswith('.mrc'):
+                fnvol += ':mrc'
+            vol = Volume()
+            vol.setFileName(fnvol)
+            vol = ih().read(vol.getFileName())
+            masscenter = vol.centerOfMass()
+            self.x = masscenter[0]
+            self.y = masscenter[1]
+            self.z = masscenter[2]
+
+        args += '--shift_to %f %f %f ' % (self.x, self.y, self.z)
+        program = "xmipp_transform_geometry"
         if not self.interp.get():
             interp = 'linear'
         else:
             interp = 'spline'
-        args = '-i "%s" -o "%s" --shift_to %f %f %f --apply_transform --dont_wrap --interp %s' % \
-               (self._getExtraPath("input_particles.xmd"), centermd, self.x.get(), self.y.get(), self.z.get(), interp)
+        args += ' --apply_transform --dont_wrap --interp %s' % interp
         if self.inv.get():
             args += ' --inverse'
         self.runJob(program, args)
@@ -104,12 +129,11 @@ def createOutputStep(self):
             outputmd = self._getExtraPath("center_particles.xmd")
         else:
             outputmd = self._getExtraPath("crop_particles.xmd")
-        outputSet.copyItems(inputParticles, updateItemCallback=self._updateItem,
-                            itemDataIterator=md.iterRows(outputmd))
+        readSetOfParticles(outputmd, outputSet)
         self._defineOutputs(outputParticles=outputSet)
-        self._defineOutputs(shiftX=pwobj.Float(self.x.get()),
-                            shiftY=pwobj.Float(self.y.get()),
-                            shiftZ=pwobj.Float(self.z.get()))
+        self._defineOutputs(shiftX=pwobj.Float(self.x),
+                            shiftY=pwobj.Float(self.y),
+                            shiftZ=pwobj.Float(self.z))
         self._defineSourceRelation(inputParticles, outputSet)
 
     # --------------------------- INFO functions --------------------------------------------
@@ -118,27 +142,21 @@ def _summary(self):
         if not hasattr(self, 'outputParticles'):
             summary.append("Output particles not ready yet.")
         else:
+            if self.option:
+                option = "User defined shift"
+            else:
+                option = "Shift to center of mass"
             if not self.interp.get():
                 interp = 'linear'
             else:
                 interp = 'spline'
-            summary.append("%d particles shifted\ninterpolation: %s" % (self.inputParticles.get().getSize(), interp))
+            summary.append("%s\ninterpolation: %s" % (option, interp))
             if self.inv.get():
                 summary.append("inverse matrix applied")
         return summary
 
-    def _methods(self):
-        return ["%d particles shifted to x = %d, y = %d, z = %d." % (self.inputParticles.get().getSize(),
-                                                                     self.x.get(), self.y.get(), self.z.get())]
-
     def _validate(self):
         for part in self.inputParticles.get().iterItems():
             if not part.hasTransform():
                 validatemsg = ['Please provide particles which have transformation matrix.']
                 return validatemsg
-
-    # --------------------------- UTLIS functions --------------------------------------------
-    def _updateItem(self, item, row):
-        newFn = row.getValue(md.MDL_IMAGE)
-        newLoc = xmippToLocation(newFn)
-        item.setLocation(newLoc)
diff --git a/xmipp3/protocols/protocol_shift_volume.py b/xmipp3/protocols/protocol_shift_volume.py
new file mode 100644
index 000000000..cf8b16c08
--- /dev/null
+++ b/xmipp3/protocols/protocol_shift_volume.py
@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+# **************************************************************************
+# *
+# * Authors:     Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *
+# *  BCU, 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, FloatParam, BooleanParam, IntParam
+import pyworkflow.object as pwobj
+from pwem.protocols import EMProtocol
+from pwem.objects import Volume
+
+
+class XmippProtShiftVolume(EMProtocol):
+    """ This protocol shifts a volume according to the input shifts"""
+
+    _label = 'shift volume'
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('inputVol', PointerParam, pointerClass='Volume', label="Volume", help='Volume to shift')
+        form.addParam('shiftBool', BooleanParam, label='Use the same shifts as for the particles?', default='True',
+                      help='Use output shifts of protocol "shift particles" which should be executed previously')
+        form.addParam('inputProtocol', PointerParam, pointerClass='XmippProtShiftParticles', allowsNull=True,
+                      label="Shift particles protocol", condition='shiftBool')
+        COND = 'not shiftBool'
+        form.addParam('x', FloatParam, label="x", condition=COND, allowsNull=True)
+        form.addParam('y', FloatParam, label="y", condition=COND, allowsNull=True)
+        form.addParam('z', FloatParam, label="z", condition=COND, allowsNull=True)
+        form.addParam('boxSizeBool', BooleanParam, label='Use original box size for the shifted volume?',
+                      default='True', help='Use input volume box size for the shifted volume.')
+        form.addParam('boxSize', IntParam, label='Final box size', condition='not boxSizeBool',
+                      help='Box size of the shifted volume.')
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('shiftStep')
+        self._insertFunctionStep('createOutputStep')
+
+    # --------------------------- STEPS functions --------------------------------------------
+    def shiftStep(self):
+        fnVol = self.inputVol.get().getFileName()
+        if not self.boxSizeBool:
+            box = self.boxSize.get()
+            self.runJob('xmipp_transform_window', '-i "%s" -o "%s" --size %d %d %d' %
+                        (fnVol, self._getExtraPath("resized_volume.mrc"), box, box, box))
+            fnVol = self._getExtraPath("resized_volume.mrc")
+        if self.shiftBool:
+            shiftprot = self.inputProtocol.get()
+            self.shiftx = shiftprot.shiftX.get()
+            self.shifty = shiftprot.shiftY.get()
+            self.shiftz = shiftprot.shiftZ.get()
+        else:
+            self.shiftx = self.x.get()
+            self.shifty = self.y.get()
+            self.shiftz = self.z.get()
+        program = "xmipp_transform_geometry"
+        args = '-i %s -o %s --shift %f %f %f --dont_wrap' % \
+               (fnVol, self._getExtraPath("shift_volume.mrc"), self.shiftx, self.shifty, self.shiftz)
+        self.runJob(program, args)
+
+    def createOutputStep(self):
+        out_vol = Volume()
+        in_vol = self.inputVol.get()
+        out_vol.setSamplingRate(in_vol.getSamplingRate())
+        out_vol.setFileName(self._getExtraPath("shift_volume.mrc"))
+        self._defineOutputs(outputVolume=out_vol)
+        self._defineOutputs(shiftX=pwobj.Float(self.shiftx),
+                            shiftY=pwobj.Float(self.shifty),
+                            shiftZ=pwobj.Float(self.shiftz))
+        self._defineSourceRelation(in_vol, out_vol)
+
+    # --------------------------- INFO functions --------------------------------------------
+    def _summary(self):
+        summary = []
+        if not hasattr(self, 'outputVolume'):
+            summary.append("Output volume not ready yet.")
+        else:
+            if self.shiftBool:
+                summary.append("Volume shift as particles in %s" % self.inputProtocol.get())
+            else:
+                summary.append("User defined shift")
+        return summary
diff --git a/xmipp3/protocols/protocol_solid_angles.py b/xmipp3/protocols/protocol_solid_angles.py
index f63263117..6a6a8e8cd 100644
--- a/xmipp3/protocols/protocol_solid_angles.py
+++ b/xmipp3/protocols/protocol_solid_angles.py
@@ -30,7 +30,7 @@
 
 import pyworkflow.protocol.params as params
 from pyworkflow import VERSION_1_1
-from pyworkflow.utils import makePath, cleanPattern, moveFile
+from pyworkflow.utils import makePath, cleanPattern, moveFile, cleanPath
 from pwem.emlib.image import ImageHandler
 from pwem.constants import ALIGN_PROJ
 from pwem.objects import Image, Volume
@@ -540,9 +540,14 @@ def createOutputStep(self):
             volumesSet = self._createSetOfVolumes()
             volumesSet.setSamplingRate(newTs)
             for i in range(2):
+                fnVol = self._getExtraPath("split_v%d.vol" % (i + 1))
+                fnMrc = self._getExtraPath("split_v%d.mrc" % (i + 1))
+                self.runJob("xmipp_image_convert", "-i %s -o %s -t vol" % (fnVol, fnMrc), numberOfMpi=1)
+                self.runJob("xmipp_image_header", "-i %s --sampling_rate %f" % (fnMrc, newTs), numberOfMpi=1)
+                cleanPath(fnVol)
+
                 vol = Volume()
-                vol.setLocation(1,
-                                self._getExtraPath("split_v%d.vol" % (i + 1)))
+                vol.setLocation(1,fnMrc)
                 volumesSet.append(vol)
 
             self._defineOutputs(outputVolumes=volumesSet)
diff --git a/xmipp3/protocols/protocol_split_volume.py b/xmipp3/protocols/protocol_split_volume.py
index 872171060..2b4c480b7 100644
--- a/xmipp3/protocols/protocol_split_volume.py
+++ b/xmipp3/protocols/protocol_split_volume.py
@@ -29,6 +29,7 @@
 
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 from pyworkflow.protocol.params import PointerParam, FloatParam, IntParam, StringParam
+from pyworkflow.utils.path import cleanPath
 from pwem.protocols import ProtClassify3D
 from pwem.objects import Volume
 from pwem.emlib.image import ImageHandler
@@ -74,11 +75,17 @@ def convertInputStep(self, inputParticlesId):
 
     def createOutput(self):
         inputParticles = self.directionalClasses.get()
+        Ts = inputParticles.getSamplingRate()
         volumesSet = self._createSetOfVolumes()
-        volumesSet.setSamplingRate(inputParticles.getSamplingRate())
+        volumesSet.setSamplingRate(Ts)
         for i in range(2):
             vol = Volume()
-            vol.setLocation(1, self._getExtraPath("split_v%d.vol"%(i+1)))
+            fnVol = self._getExtraPath("split_v%d.vol"%(i+1))
+            fnMrc = self._getExtraPath("split_v%d.mrc"%(i+1))
+            self.runJob("xmipp_image_convert","-i %s -o %s -t vol"%(fnVol,fnMrc), numberOfMpi=1)
+            self.runJob("xmipp_image_header","-i %s --sampling_rate %f"%(fnMrc, Ts), numberOfMpi=1)
+            cleanPath(fnVol)
+            vol.setLocation(1, fnMrc)
             volumesSet.append(vol)
         
         self._defineOutputs(outputVolumes=volumesSet)
diff --git a/xmipp3/protocols/protocol_split_volume_hierarchical_cluster.py b/xmipp3/protocols/protocol_split_volume_hierarchical_cluster.py
index f436e6fd5..2030f7713 100644
--- a/xmipp3/protocols/protocol_split_volume_hierarchical_cluster.py
+++ b/xmipp3/protocols/protocol_split_volume_hierarchical_cluster.py
@@ -31,7 +31,7 @@
 
 import pyworkflow.protocol.params as params
 from pyworkflow import VERSION_2_0
-from pyworkflow.utils.path import makePath, cleanPattern, moveFile
+from pyworkflow.utils.path import makePath, cleanPattern, moveFile, cleanPath
 from pwem.emlib.image import ImageHandler
 from pwem.constants import ALIGN_PROJ
 from pwem.objects import Image, Volume
@@ -766,12 +766,15 @@ def createOutputStep(self):
             volumesSet.setSamplingRate(origTs)
             for i in range(2):
                 vol = Volume()
+                fnVol = self._getExtraPath("split%d.vol" % (i + 1))
                 if origTs != lastTs:
                     newXdim = inputParticles.getXDim()
-                    self.runJob("xmipp_image_resize", "-i %s --dim %d"
-                                % (self._getExtraPath("split%d.vol" % (i + 1)),
-                                   newXdim), numberOfMpi=1)
-                vol.setLocation(1, self._getExtraPath("split%d.vol" % (i + 1)))
+                    self.runJob("xmipp_image_resize", "-i %s --dim %d" % (fnVol, newXdim), numberOfMpi=1)
+                fnMrc = self._getExtraPath("split%d.mrc" % (i + 1))
+                self.runJob("xmipp_image_convert", "-i %s -o %s -t vol" % (fnVol, fnMrc), numberOfMpi=1)
+                self.runJob("xmipp_image_header", "-i %s --sampling_rate %f" % (fnMrc, origTs), numberOfMpi=1)
+                cleanPath(fnVol)
+                vol.setLocation(1, fnMrc)
                 volumesSet.append(vol)
 
             self._defineOutputs(outputVolumes=volumesSet)
diff --git a/xmipp3/protocols/protocol_structure_map_zernike3d.py b/xmipp3/protocols/protocol_structure_map_zernike3d.py
index 3aa77ffa0..c3d476a89 100644
--- a/xmipp3/protocols/protocol_structure_map_zernike3d.py
+++ b/xmipp3/protocols/protocol_structure_map_zernike3d.py
@@ -84,8 +84,8 @@ class XmippProtStructureMapZernike3D(ProtAnalysis3D):
     """ Protocol for structure mapping based on Zernike3D. """
     _label = 'struct map - Zernike3D'
     _lastUpdateVersion = VERSION_2_0
-    OUTPUT_SUFFIX = '_%d_crop.vol'
-    ALIGNED_VOL = 'vol%dAligned.vol'
+    OUTPUT_SUFFIX = '_%d_crop.mrc'
+    ALIGNED_VOL = 'vol%dAligned.mrc'
 
     def __init__(self, **args):
         ProtAnalysis3D.__init__(self, **args)
@@ -241,11 +241,11 @@ def deformStep(self, inputVolFn, refVolFn, i, j, step_id):
                                                             + self.OUTPUT_SUFFIX % (i + 1)))
         else:
             fnOut_aux = self._getTmpPath(self.ALIGNED_VOL % (i + 1))
-        refVolFn = self._getTmpPath("reference_%d.vol" % step_id)
-        fnOut = self._getTmpPath("input_%d.vol" % step_id)
+        refVolFn = self._getTmpPath("reference_%d.mrc" % step_id)
+        fnOut = self._getTmpPath("input_%d.mrc" % step_id)
         copyFile(refVolFn_aux, refVolFn)
         copyFile(fnOut_aux, fnOut)
-        fnOut2 = self._getTmpPath('vol%dDeformedTo%d.vol' % (i + 1, j + 1))
+        fnOut2 = self._getTmpPath('vol%dDeformedTo%d.mrc' % (i + 1, j + 1))
 
         params = ' -i %s -r %s -o %s --l1 %d --l2 %d --sigma "%s" --oroot %s --regularization %f' %\
                  (fnOut, refVolFn, fnOut2, self.l1.get(), self.l2.get(), self.sigma.get(),
@@ -266,7 +266,7 @@ def deformStep(self, inputVolFn, refVolFn, i, j, step_id):
         cleanPath(fnOut2)
 
     def deformationMatrix(self, volList):
-        cleanPath(self._getTmpPath("*.vol"))
+        cleanPath(self._getTmpPath("*.mrc"))
         numVol = len(volList)
         self.distanceMatrix = np.zeros((numVol, numVol))
         for i in range(numVol):
diff --git a/xmipp3/protocols/protocol_subtract_projection.py b/xmipp3/protocols/protocol_subtract_projection.py
new file mode 100644
index 000000000..5843416a5
--- /dev/null
+++ b/xmipp3/protocols/protocol_subtract_projection.py
@@ -0,0 +1,146 @@
+# -*- 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 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.protocols import EMProtocol
+from xmipp3.convert import writeSetOfParticles
+
+
+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. """
+
+    _label = 'subtract projection'
+    INPUT_PARTICLES = "input_particles.xmd"
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    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.')
+        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')
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+    def _insertAllSteps(self):
+        self._insertFunctionStep('convertStep')
+        self._insertFunctionStep('subtractionStep')
+        self._insertFunctionStep('createOutputStep')
+
+    # --------------------------- STEPS functions --------------------------------------------
+    def convertStep(self):
+        writeSetOfParticles(self.particles.get(), self._getExtraPath(self.INPUT_PARTICLES))
+        
+    def subtractionStep(self):
+        vol = self.vol.get().clone()
+        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()
+        self.runJob(program, 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)))
+        self._defineOutputs(outputParticles=outputSet)
+        self._defineSourceRelation(inputSet, outputSet)
+
+    # --------------------------- INFO functions --------------------------------------------
+    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())
+        return summary
+
+    def _methods(self):
+        methods = []
+        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())
+        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_validate_fscq.py b/xmipp3/protocols/protocol_validate_fscq.py
index 99cd62eda..8dde147c4 100644
--- a/xmipp3/protocols/protocol_validate_fscq.py
+++ b/xmipp3/protocols/protocol_validate_fscq.py
@@ -369,76 +369,69 @@ def assignPdbStep(self):
         params += ' --radius 0.8' 
         params += ' --md %s' % self._getFileName(MD2_MEANS) 
         self.runJob('xmipp_pdb_label_from_volume', params)    
-        
-    def createOutputStep(self):    
-        
-        volume=Volume()
-        volume.setFileName(self._getFileName(RESTA_FILE_MRC))
-        volume.setSamplingRate(self.inputVolume.get().getSamplingRate())
-        volume.setOrigin(self.inputVolume.get().getOrigin(True).clone())
-        self._defineOutputs(fscq_Volume=volume)
-        self._defineTransformRelation(self.inputVolume, volume)
-        
-        #mean values of FSC-Q
-              
+
+    def _getMetrics(self):
+        """ Internal method to compute some metrics. """
+        # mean values of FSC-Q
+
         mtd = md.MetaData()
         mtd.read(self._getFileName(MD_MEANS))
-            
-        mean = mtd.getValue(MDL_VOLUME_SCORE1,1)
-        meanA = mtd.getValue(MDL_VOLUME_SCORE2,1)  
-        
-        #means value for map divided by resolution (FSC-Qr)
+
+        mean = mtd.getValue(MDL_VOLUME_SCORE1, 1)
+        meanA = mtd.getValue(MDL_VOLUME_SCORE2, 1)
+
+        # means value for map divided by resolution (FSC-Qr)
         mtd2 = md.MetaData()
         mtd2.read(self._getFileName(MD2_MEANS))
-            
-        mean2 = mtd2.getValue(MDL_VOLUME_SCORE1,1)
-        meanA2 = mtd2.getValue(MDL_VOLUME_SCORE2,1)
-              
-        #Setting the mean fsc-q for the summary
-        self.mean = Float(mean)
-        self._store(self) 
-        self.meanA = Float(meanA)
-        self._store(self) 
-        
-        self.mean2 = Float(mean2)
-        self._store(self) 
-        self.meanA2 = Float(meanA2)
-        self._store(self)
 
-        
-        #statistic from fnal pdb with fsc-q
-        #Number of atoms greater or less than 0.5
-        total_atom=0
-        fscq_greater=0
-        fscq_less=0
+        mean2 = mtd2.getValue(MDL_VOLUME_SCORE1, 1)
+        meanA2 = mtd2.getValue(MDL_VOLUME_SCORE2, 1)
+
+        # statistic from fnal pdb with fsc-q
+        # Number of atoms greater or less than 0.5
+        total_atom = 0
+        fscq_greater = 0
+        fscq_less = 0
         with open(self._getFileName(PDB_VALUE_FILE)) as f:
             lines_data = f.readlines()
-            for j,lin in enumerate(lines_data):
-                
-                if ( lin.startswith('ATOM') or lin.startswith('HETATM') ):
-                    
-                    total_atom=total_atom+1
+            for j, lin in enumerate(lines_data):
+
+                if (lin.startswith('ATOM') or lin.startswith('HETATM')):
+
+                    total_atom = total_atom + 1
                     fscq_atom = float(lin[54:60])
-                    
-                    if (fscq_atom>0.5):
-                        fscq_greater=fscq_greater+1
-                        
-                    if (fscq_atom<-0.5):
-                        fscq_less=fscq_less+1
-                        
-        porc_greater=(fscq_greater*100)/total_atom
-        porc_less=(fscq_less*100)/total_atom
- 
-        self.total_atom = Integer(total_atom)
-        self._store(self)
-        self.fscq_greater=Integer(fscq_greater)
-        self._store(self)
-        self.fscq_less=Integer(fscq_less)
-        self._store(self)        
-        self.porc_greater=Float(porc_greater)
-        self._store(self)
-        self.porc_less=Float(porc_less)
-        self._store(self)           
+
+                    if (fscq_atom > 0.5):
+                        fscq_greater = fscq_greater + 1
+
+                    if (fscq_atom < -0.5):
+                        fscq_less = fscq_less + 1
+
+        porc_greater = (fscq_greater * 100) / total_atom
+        porc_less = (fscq_less * 100) / total_atom
+
+        return {
+            'mean': Float(mean),
+            'meanA': Float(meanA),
+            'mean2': Float(mean2),
+            'meanA2': Float(meanA2),
+            'total_atom': Integer(total_atom),
+            'fscq_greater': Integer(fscq_greater),
+            'fscq_less': Integer(fscq_less),
+            'porc_greater': Float(porc_greater),
+            'porc_less': Float(porc_less)
+        }
+
+    def createOutputStep(self):
+        metrics = self._getMetrics()
+        volume = Volume()
+        volume.setFileName(self._getFileName(RESTA_FILE_MRC))
+        volume.setSamplingRate(self.inputVolume.get().getSamplingRate())
+        volume.setOrigin(self.inputVolume.get().getOrigin(True).clone())
+
+        self._defineOutputs(fscq_Volume=volume, **metrics)
+        self._defineTransformRelation(self.inputVolume, volume)
+        
 
     # --------------------------- INFO functions ------------------------------
 
diff --git a/xmipp3/protocols/protocol_validate_overfitting.py b/xmipp3/protocols/protocol_validate_overfitting.py
index e4575241c..1527a7a26 100644
--- a/xmipp3/protocols/protocol_validate_overfitting.py
+++ b/xmipp3/protocols/protocol_validate_overfitting.py
@@ -326,8 +326,6 @@ def reconstructionStep(self, numberOfImages, fractionCounter,
                 params += ' --sym %s' % self.symmetryGroup.get()
                 params += ' --max_resolution %0.3f' % self.maxRes
                 params += ' --padding 2'
-                params += ' --thr 1'
-                # params += ' --thr %d' % self.numberOfThreads.get()
                 params += ' --sampling %f' % Ts
 
                 if self.useGpu.get():
diff --git a/xmipp3/protocols/protocol_volume_consensus.py b/xmipp3/protocols/protocol_volume_consensus.py
index 62d457413..05aab7485 100644
--- a/xmipp3/protocols/protocol_volume_consensus.py
+++ b/xmipp3/protocols/protocol_volume_consensus.py
@@ -25,6 +25,7 @@
 # *
 # **************************************************************************
 
+from os.path import exists
 from pyworkflow.protocol.params import MultiPointerParam
 from pwem.objects import Volume, Transform
 from pwem.protocols import ProtInitialVolume
@@ -40,7 +41,6 @@ class XmippProtVolConsensus(ProtInitialVolume):
 
     # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
-
         form.addSection(label='Input')
         form.addParam('vols', MultiPointerParam, pointerClass='Volume', label="Volumes",
                       help='Select the volumes for the consensus.')
@@ -69,11 +69,15 @@ def createOutputStep(self):
         outVol = Volume()
         outVol.setSamplingRate(self.vols[0].get().getSamplingRate())
         outVol.setFileName(self._getExtraPath("consensus_volume.mrc"))
-        outVol2 = Volume()
-        outVol2.setSamplingRate(self.vols[0].get().getSamplingRate())
-        outVol2.setFileName(self._getExtraPath("consensus_volume_diff.mrc"))
-        self._defineOutputs(outputVolume=outVol)
-        self._defineOutputs(outputVolumeDiff=outVol2)
+        if not exists(self._getExtraPath("consensus_volume.mrc")):
+            raise NoOutputGenerated("Consensus volume NOT generated, please check input volumes to ensure they have "
+                                    "equal box size, voxel size and origin.")
+        else:
+            outVol2 = Volume()
+            outVol2.setSamplingRate(self.vols[0].get().getSamplingRate())
+            outVol2.setFileName(self._getExtraPath("consensus_volume_diff.mrc"))
+            self._defineOutputs(outputVolume=outVol)
+            self._defineOutputs(outputVolumeDiff=outVol2)
 
     def createChimeraScript(self):
         fnRoot = "extra/"
@@ -104,5 +108,20 @@ def _methods(self):
                            (self.vols.getSize(), self.vols[0].get().getSamplingRate()))
         return methods
 
+    def _validate(self):
+        errors = []
+        voxel_size = []
+        for i, vol in enumerate(self.vols):
+            voxel_size.append(vol.get().getSamplingRate())
+        result = all(element == voxel_size[0] for element in voxel_size)
+        if not result:
+            errors.append('Pixel size should be the same for all input volumes.')
+        return errors
+
     def _citations(self):
         return ['Fernandez-Gimenez2021']
+
+
+class NoOutputGenerated(Exception):
+    """No output generation error"""
+    pass
diff --git a/xmipp3/protocols/protocol_volume_deform_zernike3d.py b/xmipp3/protocols/protocol_volume_deform_zernike3d.py
index 306607bcd..27003f6da 100644
--- a/xmipp3/protocols/protocol_volume_deform_zernike3d.py
+++ b/xmipp3/protocols/protocol_volume_deform_zernike3d.py
@@ -82,23 +82,24 @@ def _defineParams(self, form):
     def _createFilenameTemplates(self):
         """ Centralize how files are called """
         myDict = {
-            'fnRefVol': self._getExtraPath('ref_volume.vol'),
-            'fnInputVol': self._getExtraPath('input_volume.vol'),
-            'fnInputFilt': self._getExtraPath('input_volume_filt.vol'),
-            'fnRefFilt': self._getExtraPath('ref_volume_filt.vol'),
-            'fnOutVol': self._getExtraPath('vol1DeformedTo2.vol')
+            'fnRefVol': self._getExtraPath('ref_volume.mrc'),
+            'fnInputVol': self._getExtraPath('input_volume.mrc'),
+            'fnInputFilt': self._getExtraPath('input_volume_filt.mrc'),
+            'fnRefFilt': self._getExtraPath('ref_volume_filt.mrc'),
+            'fnOutVol': self._getExtraPath('vol1DeformedTo2.mrc')
                  }
         self._updateFilenamesDict(myDict)
 
     # --------------------------- INSERT steps functions -----------------------
     def _insertAllSteps(self):
         self._createFilenameTemplates()
-        self._insertFunctionStep("convertStep")
-        self._insertFunctionStep("deformStep")
-        self._insertFunctionStep("createOutputStep")
+        self._insertFunctionStep(self.convertInputStep)
+        self._insertFunctionStep(self.deformStep)
+        self._insertFunctionStep(self.convertOutputStep)
+        self._insertFunctionStep(self.createOutputStep)
 
     # --------------------------- STEPS functions ------------------------------
-    def convertStep(self):
+    def convertInputStep(self):
         fnInputVol = self._getFileName('fnInputVol')
         fnRefVol = self._getFileName('fnRefVol')
 
@@ -153,6 +154,10 @@ def deformStep(self):
                 params = params + ' --thr %d' % self.numberOfThreads.get()
             self.runJob("xmipp_volume_deform_sph", params)
 
+    def convertOutputStep(self):
+        fnOutVol = self._getFileName('fnOutVol')
+        params = ' -i %s --sampling_rate %f' % (fnOutVol, self.newTs)
+        self.runJob("xmipp_image_header", params)
 
     def createOutputStep(self):
         if self.targetResolution.get() != 3.0:
diff --git a/xmipp3/protocols/protocol_volume_strain.py b/xmipp3/protocols/protocol_volume_strain.py
index 8d99f611d..a6f8a7c73 100644
--- a/xmipp3/protocols/protocol_volume_strain.py
+++ b/xmipp3/protocols/protocol_volume_strain.py
@@ -58,6 +58,9 @@ def _defineParams(self, form):
                         'If no symmetry is present, give c1')
     
     #--------------------------- INSERT steps functions --------------------------------------------
+    def _getFileName(self, fnRoot, key, **kwargs):
+        return "%s_%s.mrc"%(fnRoot,key)
+
     def _insertAllSteps(self):
         fnVol0 = self.inputVolume0.get().getFileName()
         fnVolF = self.inputVolumeF.get().getFileName()
@@ -77,49 +80,60 @@ def calculateStrain(self, fnVol0, fnVolF, fnMask):
     
     def prepareOutput(self):
         volDim = self.inputVolume0.get().getDim()[0]
+        Ts=self.inputVolume0.get().getSamplingRate()
         fnRoot=self._getExtraPath('result')
-        self.runJob("xmipp_image_convert", "-i %s_initial.raw#%d,%d,%d,0,float -o %s_initial.vol"%
-                    (fnRoot,volDim,volDim,volDim,fnRoot))
-        self.runJob("xmipp_transform_mirror","-i %s_initial.vol --flipX"%fnRoot)
-        self.runJob("xmipp_image_convert", "-i %s_final.raw#%d,%d,%d,0,float -o %s_final.vol"%
-                    (fnRoot,volDim,volDim,volDim,fnRoot))
-        self.runJob("xmipp_transform_mirror","-i %s_final.vol --flipX"%fnRoot)
-        self.runJob("xmipp_image_convert", "-i %s_initialDeformedToFinal.raw#%d,%d,%d,0,float -o %s_initialDeformedToFinal.vol"%
-                    (fnRoot,volDim,volDim,volDim,fnRoot))
-        self.runJob("xmipp_transform_mirror","-i %s_initialDeformedToFinal.vol --flipX"%fnRoot)
-        self.runJob("xmipp_image_convert", "-i %s_strain.raw#%d,%d,%d,0,float -o %s_strain.vol"%
-                    (fnRoot,volDim,volDim,volDim,fnRoot))
-        self.runJob("xmipp_transform_mirror","-i %s_strain.vol --flipX"%fnRoot)
-        self.runJob("xmipp_image_convert", "-i %s_localrot.raw#%d,%d,%d,0,float -o %s_localrot.vol"%
-                    (fnRoot,volDim,volDim,volDim,fnRoot))
-        self.runJob("xmipp_transform_mirror","-i %s_localrot.vol --flipX"%fnRoot)
+
+        def symmetrize(key):
+            self.runJob("xmipp_transform_symmetrize", "-i %s --sym %s --dont_wrap" % \
+                        (self._getFileName(fnRoot, key), self.symmetryGroup.get()))
+
+        def changeSamplingRate(key):
+            self.runJob("xmipp_image_header", "-i %s --sampling_rate %f" % (self._getFileName(fnRoot, key), Ts))
+
+        def convert(key):
+            self.runJob("xmipp_image_convert", "-i %s_%s.raw#%d,%d,%d,0,float -o %s" %
+                        (fnRoot, key, volDim, volDim, volDim, self._getFileName(fnRoot, key)))
+            self.runJob("xmipp_transform_mirror", "-i %s --flipX" % self._getFileName(fnRoot, key))
+            changeSamplingRate(key)
+
+        convert("initial")
+        convert("final")
+        convert("initialDeformedToFinal")
+        convert("strain")
+        convert("localrot")
+
         self.runJob("rm","-f "+self._getExtraPath('result_*.raw'))
         if self.symmetryGroup!="c1":
-            self.runJob("xmipp_transform_symmetrize","-i %s --sym %s --dont_wrap"%(fnRoot+"_strain.vol",self.symmetryGroup.get()))
-            self.runJob("xmipp_transform_symmetrize","-i %s --sym %s --dont_wrap"%(fnRoot+"_localrot.vol",self.symmetryGroup.get()))
-    
+            symmetrize("strain")
+            symmetrize("localrot")
+            changeSamplingRate("strain")
+            changeSamplingRate("localrot")
+
     def createChimeraScript(self):
         fnRoot = "extra/result"
         scriptFile = self._getPath('result') + '_strain_chimera.cmd'
+
+        openStr = "open %s\n"
+
         fhCmd = open(scriptFile, 'w')
-        fhCmd.write("open %s\n" % (fnRoot+"_final.vol"))
-        fhCmd.write("open %s\n" % (fnRoot+"_strain.vol"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"final"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"strain"))
         fhCmd.write("vol #1 hide\n")
         fhCmd.write("scolor #0 volume #1 cmap rainbow reverseColors True\n")
         fhCmd.close()
 
         scriptFile = self._getPath('result') + '_localrot_chimera.cmd'
         fhCmd = open(scriptFile, 'w')
-        fhCmd.write("open %s\n" % (fnRoot+"_final.vol"))
-        fhCmd.write("open %s\n" % (fnRoot+"_localrot.vol"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"final"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"localrot"))
         fhCmd.write("vol #1 hide\n")
         fhCmd.write("scolor #0 volume #1 cmap rainbow reverseColors True\n")
         fhCmd.close()
 
         scriptFile = self._getPath('result') + '_morph_chimera.cmd'
         fhCmd = open(scriptFile, 'w')
-        fhCmd.write("open %s\n" % (fnRoot+"_initial.vol"))
-        fhCmd.write("open %s\n" % (fnRoot+"_final.vol"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"initial"))
+        fhCmd.write(openStr % self._getFileName(fnRoot,"final"))
         fhCmd.write("vol #0 hide\n")
         fhCmd.write("vol #1 hide\n")
         fhCmd.write("vop morph #0,1 frames 50\n")
diff --git a/xmipp3/viewers/viewer_resolution_bfactor.py b/xmipp3/viewers/viewer_resolution_bfactor.py
index 3bb6895da..7c9a75e19 100644
--- a/xmipp3/viewers/viewer_resolution_bfactor.py
+++ b/xmipp3/viewers/viewer_resolution_bfactor.py
@@ -29,7 +29,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pyworkflow.protocol.params import LabelParam
+from pyworkflow.protocol.params import LabelParam, FloatParam
 from pyworkflow.viewer import ProtocolViewer, DESKTOP_TKINTER
 
 from pwem import emlib
@@ -53,20 +53,35 @@ def __init__(self, *args, **kwargs):
     def _defineParams(self, form):
         form.addSection(label='Visualization')
 
-        line = form.addLine("Color scale options:",
+        groupColorScaleOptions = form.addGroup('Color scale options and representation')
+
+        line = groupColorScaleOptions.addLine("Range normalized bfactor:",
                             help="Options to define the color scale limits and color set")
 
-        line.addParam('highest', FloatParam, default=self.maxNormRes,
+        line.addParam('highestBF', FloatParam, allowsNull =True,
+                      label="Highest",
+                      help="Highest value for the scale")
+
+        line.addParam('lowestBF', FloatParam, allowsNull =True,
+                      label="Lowest",
+                      help="lowest value of the scale.")
+
+        line2 = groupColorScaleOptions.addLine("Range normalized local resolution:",
+                                              help="Options to define the color scale limits and color set")
+
+        line2.addParam('highestLR', FloatParam, allowsNull =True,
                       label="Highest",
                       help="Highest value for the scale")
 
-        line.addParam('lowest', FloatParam, default=self.minNormRes,
+        line2.addParam('lowestLR', FloatParam, allowsNull =True,
                       label="Lowest",
                       help="lowest value of the scale.")
 
-        line.addParam('doShowColorBands', LabelParam,
+        groupColorScaleOptions.addParam('doShowColorBands', LabelParam,
                       label="Show bfactor-local resolution comparison")
 
+
+
     def _getVisualizeDict(self):
         return {'doShowColorBands': self._showColorBands,
                 }
@@ -89,18 +104,26 @@ def _showColorBands(self, param=None):
         lr = np.array(lr)
         bf = np.array(bf)
         r = np.array(r)
+
+        lowBF = self.lowestBF if self.lowestBF.get() else bf[0]
+        highBF = self.highestBF if self.highestBF.get() else bf[-1]
+        lowLR = self.lowestLR if self.lowestLR.get() else lr[0]
+        highLR = self.highestLR if self.highestLR.get() else lr[-1]
+
+        plt.figure()
         plt.subplot(211)
         #The magic numbers of 0 and 40 define the size of the vertical bands, they provide good visualization aspect
-        plt.imshow(lr.reshape(1, len(lr)), cmap=plt.cm.viridis, extent=[r[0], r[-1], 0, 40])
+        plt.imshow(lr.reshape(1, len(lr)), vmin=lowLR, vmax=highLR, cmap=plt.cm.viridis, extent=[np.amin(r), np.amax(r), 0, 80])
+
         plt.xlabel('Residue')
         plt.title('Normalized Local Resolution')
         plt.colorbar()
 
         plt.subplot(212)
-        plt.imshow(bf.reshape(1, len(bf)), cmap=plt.cm.viridis, extent=[r[0], r[-1], 0, 40])
+
+        plt.imshow(bf.reshape(1, len(bf)), vmin=lowBF, vmax=highBF, cmap=plt.cm.viridis, extent=[np.amin(r), np.amax(r), 0, 80])
         plt.xlabel('Residue')
         plt.title('B-factor')
         plt.colorbar()
 
         plt.show()
-
diff --git a/xmipp3/viewers/viewer_resolution_fso.py b/xmipp3/viewers/viewer_resolution_fso.py
index 0d05b3bb1..94cf1dca8 100644
--- a/xmipp3/viewers/viewer_resolution_fso.py
+++ b/xmipp3/viewers/viewer_resolution_fso.py
@@ -195,8 +195,9 @@ def _showAnisotropyCurve(self, paramName=None):
         xTitle = 'Resolution (1/A)'
         yTitle = 'Anisotropy (a.u.)'
         mdLabelX = emlib.MDL_RESOLUTION_FREQ
-        mdLabelY = emlib.MDL_RESOLUTION_FSO
-        self._plotCurveAnisotropy(fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY)
+        mdLabelY1 = emlib.MDL_RESOLUTION_FSO
+        mdLabelY2 = emlib.MDL_RESOLUTION_ANISOTROPY
+        self._plotCurveAnisotropy(fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY1, mdLabelY2)
 
     def _showDirectionalFilter(self, param=None):
         """
@@ -292,7 +293,7 @@ def interpolRes(self, thr, x, y):
 
         return resInterp, okToPlot
 
-    def _plotCurveAnisotropy(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY):
+    def _plotCurveAnisotropy(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY1, mdLabelY2):
         """
         This function is called by _showAnisotropyCurve
         It shows the FSO curve in terms of the resolution
@@ -303,9 +304,11 @@ def _plotCurveAnisotropy(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY):
         xplotter.plot_title_fontsize = 11
 
         a = xplotter.createSubPlot(title, xTitle, yTitle, 1, 1)
-        xplotter.plotMdFile(md, mdLabelX, mdLabelY, 'g')
+        xplotter.plotMdFile(md, mdLabelX, mdLabelY1, 'g')
+
+        xx, yy = self._prepareDataForPlot(md, mdLabelX, mdLabelY1)
+        _, yyBingham = self._prepareDataForPlot(md, mdLabelX, mdLabelY2)
 
-        xx, yy = self._prepareDataForPlot(md, mdLabelX, mdLabelY)
         from matplotlib.ticker import FuncFormatter
         a.axes.xaxis.set_major_formatter(FuncFormatter(self._formatFreq))
         a.axes.set_ylim([-0.1, 1.1])
@@ -321,6 +324,8 @@ def _plotCurveAnisotropy(self, fnmd, title, xTitle, yTitle, mdLabelX, mdLabelY):
         res_05, okToPlot_05 = self.interpolRes(0.5, xx, yy)
         res_09, okToPlot_09 = self.interpolRes(0.9, xx, yy)
 
+        a.axes.plot(xx, yyBingham, 'r--')
+
         if (okToPlot_01 and okToPlot_05 and okToPlot_09):
             textstr = str(0.9) + ' --> ' + str("{:.2f}".format(res_09)) + 'A\n' + str(0.5) + ' --> ' + str(
                 "{:.2f}".format(res_05)) + 'A\n' + str(0.1) + ' --> ' + str("{:.2f}".format(res_01)) + 'A'
diff --git a/xmipp3/viewers/viewer_swarm.py b/xmipp3/viewers/viewer_swarm.py
index d25d797cb..e1dcaf97c 100644
--- a/xmipp3/viewers/viewer_swarm.py
+++ b/xmipp3/viewers/viewer_swarm.py
@@ -45,7 +45,7 @@ def __init__(self, **args):
 
     def _visualize(self, obj, **args):
         import os
-        fnVolume = self.protocol._getExtraPath("volumeAvg.vol")
+        fnVolume = self.protocol._getExtraPath("volumeAvg.mrc")
         if os.path.exists(fnVolume):
             fnDir = self.protocol._getExtraPath()
             samplingRate=self.protocol.readInfoField(fnDir,"sampling",emlib.MDL_SAMPLINGRATE)
diff --git a/xmipp3/wizards.py b/xmipp3/wizards.py
index 10217c725..7105e259b 100644
--- a/xmipp3/wizards.py
+++ b/xmipp3/wizards.py
@@ -528,8 +528,9 @@ def show(self, form):
 class ColorScaleWizard(ColorScaleWizardBase):
         _targets = ColorScaleWizardBase.defineTargets(XmippMonoResViewer, XmippResDeepResViewer, XmippProtFSOViewer)
 
+
 class XmippSelectPointinVolWizard(EmWizard):
-    _targets = [(XmippProtShiftParticles, ['x', 'y', 'z'])]
+    _targets = [(XmippProtShiftParticles, ['xin', 'yin', 'zin'])]
 
     def show(self, form):
         protocol = form.protocol
@@ -539,8 +540,7 @@ def show(self, form):
             return
         plt = MaskVolumeWizard(volume.getFileName())
         plt.initializePlot()
-        form.setVar('x', int(plt.origin[2]))
-        form.setVar('y', int(plt.origin[1]))
-        form.setVar('z', int(plt.origin[0]))
+        form.setVar('xin', int(plt.origin[2]))
+        form.setVar('yin', int(plt.origin[1]))
+        form.setVar('zin', int(plt.origin[0]))
         del plt
-


Test patch:
diff --git a/xmipp3/tests/test_protocols_continuousflex.py b/xmipp3/tests/test_protocols_continuousflex.py
new file mode 100644
index 000000000..5dc3bf639
--- /dev/null
+++ b/xmipp3/tests/test_protocols_continuousflex.py
@@ -0,0 +1,103 @@
+# **************************************************************************
+# * Authors:     Mohamad Harastani (mohamad.harastani@upmc.fr)
+# * IMPMC, Sorbonne University
+# *
+# * 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 continuousflex.protocols import FlexProtAlignmentNMAVol
+from pwem.protocols import ProtImportPdb
+from pwem.tests.workflows import TestWorkflow
+from pyworkflow.tests import setupTestProject, DataSet
+from continuousflex.protocols import (FlexProtNMA, FlexProtSynthesizeSubtomo, FlexProtSynthesizeImages, NMA_CUTOFF_ABS,
+                                      FlexProtAlignmentNMA, FlexProtSubtomogramAveraging)
+from continuousflex.protocols.protocol_subtomogrmas_synthesize import MODE_RELATION_3CLUSTERS
+
+
+class TestContinuousFlexBasics(TestWorkflow):
+    """ Test protocol for HEMNMA, HEMNMA-3D, StructMap, TomoFlow and STA backend functions. """
+    @classmethod
+    def setUpClass(cls):
+        # Create a new project
+        setupTestProject(cls)
+        cls.ds = DataSet.getDataSet('nma_V2.0')
+
+    def test_all(self):
+        # ------------------------------------------------
+        # Import a Pdb -> NMA
+        # ------------------------------------------------
+        protImportPdb = self.newProtocol(ProtImportPdb, inputPdbData=1,
+                                         pdbFile=self.ds.getFile('pdb'))
+        protImportPdb.setObjLabel('AK.pdb')
+        self.launchProtocol(protImportPdb)
+        # Launch NMA for PDB imported
+        protNMA = self.newProtocol(FlexProtNMA,
+                                   cutoffMode=NMA_CUTOFF_ABS)
+        protNMA.inputStructure.set(protImportPdb.outputPdb)
+        protNMA.setObjLabel('NMA')
+        self.launchProtocol(protNMA)
+        # ------------------------------------------------
+        # Create few volumes for testing
+        # ------------------------------------------------
+        protSynthesize = self.newProtocol(FlexProtSynthesizeSubtomo,
+                                          modeList='7-8',
+                                          numberOfVolumes=2,
+                                          modeRelationChoice=MODE_RELATION_3CLUSTERS,
+                                          noiseCTFChoice=1, # not adding noise
+                                          volumeSize=16)
+        protSynthesize.inputModes.set(protNMA.outputModes)
+        protSynthesize.setObjLabel('test volumes')
+        self.launchProtocol(protSynthesize)
+        # ------------------------------------------------
+        # Lunch HEMNMA-3D (this will test xmipp_nma_alignment_vol, also useful for StructMap)
+        # ------------------------------------------------
+        protAlignment = self.newProtocol(FlexProtAlignmentNMAVol,
+                                         modeList='7-8')
+        protAlignment.inputModes.set(protNMA.outputModes)
+        protAlignment.inputVolumes.set(protSynthesize.outputVolumes)
+        protAlignment.setObjLabel('HEMNMA-3D & StructMap backend')
+        self.launchProtocol(protAlignment)
+        # ------------------------------------------------
+        # Lunch STA on the synthesized subtomograms (this will test xmipp_volumeset_align, also useful for TomoFlow)
+        # ------------------------------------------------
+        protStA = self.newProtocol(FlexProtSubtomogramAveraging,
+                                   NumOfIters=1)
+        protStA.inputVolumes.set(protSynthesize.outputVolumes)
+        protStA.setObjLabel('STA and TomoFlow backend')
+        self.launchProtocol(protStA)
+        # ------------------------------------------------
+        # Create few images for testing
+        # ------------------------------------------------
+        protSynthesize = self.newProtocol(FlexProtSynthesizeImages,
+                                          modeList='7-8',
+                                          numberOfVolumes=2,
+                                          modeRelationChoice=MODE_RELATION_3CLUSTERS,
+                                          noiseCTFChoice=1, # not adding noise
+                                          volumeSize=32)
+        protSynthesize.inputModes.set(protNMA.outputModes)
+        protSynthesize.setObjLabel('test images')
+        self.launchProtocol(protSynthesize)
+        # ------------------------------------------------
+        # Lunch HEMNMA on the synthesized images (this will test xmipp_nma_alignment)
+        # ------------------------------------------------
+        protAlignment = self.newProtocol(FlexProtAlignmentNMA,
+                                         modeList='7-8')
+        protAlignment.inputModes.set(protNMA.outputModes)
+        protAlignment.inputParticles.set(protSynthesize.outputImages)
+        protAlignment.setObjLabel('HEMNMA backend')
+        self.launchProtocol(protAlignment)
\ No newline at end of file
diff --git a/xmipp3/tests/test_protocols_highres.py b/xmipp3/tests/test_protocols_highres.py
index ac62cd1d8..35d5d6a8e 100644
--- a/xmipp3/tests/test_protocols_highres.py
+++ b/xmipp3/tests/test_protocols_highres.py
@@ -88,6 +88,7 @@ def test(self):
                                    inputVolumes=self.protImportVol.outputVolume,
                                    particleRadius=180,
                                    symmetryGroup="i1",
+                                   nextResolutionCriterion=0.143,
                                    alignmentMethod=XmippProtReconstructHighRes.AUTOMATIC_ALIGNMENT,
                                    maximumTargetResolution="15 10 7",
                                    numberOfMpi=8)
diff --git a/xmipp3/tests/test_protocols_xmipp_3d.py b/xmipp3/tests/test_protocols_xmipp_3d.py
index b608cbcd3..04b515a77 100644
--- a/xmipp3/tests/test_protocols_xmipp_3d.py
+++ b/xmipp3/tests/test_protocols_xmipp_3d.py
@@ -1198,19 +1198,6 @@ def testXmippValidateNonTilt(self):
         self.assertIsNotNone(protValidate.outputVolumes, "There was a problem with Validate Non-Tilt")
 
 
-if __name__ == "__main__":
-    if len(sys.argv) > 1:
-        className = sys.argv[1]
-        cls = globals().get(className, None)
-        if cls:
-            suite = unittest.TestLoader().loadTestsFromTestCase(cls)
-            unittest.TextTestRunner(verbosity=2).run(suite)
-        else:
-            print("Test: '%s' not found." % className)
-    else:
-        unittest.main()
-
-
 class TestXmippVolSubtraction(TestXmippBase):
 
     @classmethod
@@ -1336,3 +1323,349 @@ def testXmippVolSub(self):
         self.launchProtocol(protVolConsensus)
         self.assertIsNotNone(protVolConsensus.outputVolume,
                              "There was a problem with Volumes consensus")
+
+
+class TestXmippVolPhantom(TestXmippBase):
+
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+
+    def testXmippPhantomVol(self):
+        protCreatePhantom = self.newProtocol(XmippProtPhantom)
+        self.launchProtocol(protCreatePhantom)
+        self.assertIsNotNone(protCreatePhantom.getFiles(),  "There was a problem with phantom creation")
+        self.assertEqual(protCreatePhantom.outputVolume.getSamplingRate(), 4,
+                         "There was a problem with the sampling rate value of output phantom")
+        self.assertEqual(protCreatePhantom.outputVolume.getDim(), (40, 40, 40),
+                         "There was a problem with the dimensions of output phantom")
+
+
+class TestXmippShiftParticlesAndVolume(TestXmippBase):
+
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+        cls.dataset = DataSet.getDataSet('xmipp_tutorial')
+        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,
+                                          samplingRate=7.08)
+        self.launchProtocol(protImportVol)
+        self.assertIsNotNone(protImportVol.getFiles(),
+                             "There was a problem with the volume import")
+
+        protCreateGallery = self.newProtocol(XmippProtCreateGallery,
+                                             inputVolume=protImportVol.outputVolume,
+                                             rotStep=15.0,
+                                             tiltStep=90.0)
+        self.launchProtocol(protCreateGallery)
+        self.assertIsNotNone(protCreateGallery.getFiles(),
+                             "There was a problem with create gallery")
+
+        protShiftParticles = self.newProtocol(XmippProtShiftParticles,
+                                              inputParticles=protCreateGallery.outputReprojections,
+                                              xin=2.0, yin=3.0, zin=4.0)
+        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.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"))
+
+        protCreateMask = self.newProtocol(XmippProtCreateMask3D,
+                                          inputVolume=protImportVol.outputVolume,
+                                          threshold=0.1)
+        self.launchProtocol(protCreateMask)
+        self.assertIsNotNone(protCreateMask.getFiles(),
+                             "There was a problem with the 3D mask ")
+
+        protShiftParticlesCenterOfMass = self.newProtocol(XmippProtShiftParticles,
+                                                          inputParticles=protCreateGallery.outputReprojections,
+                                                          inputMask=protCreateMask.outputMask,
+                                                          option=False)
+        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.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"))
+
+        protShiftVolPart = self.newProtocol(XmippProtShiftVolume,
+                                            inputVol=protImportVol.outputVolume,
+                                            inputProtocol=protShiftParticles)
+        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"))
+
+        protShiftVolCrop = self.newProtocol(XmippProtShiftVolume,
+                                            inputVol=protImportVol.outputVolume,
+                                            shiftBool=False,
+                                            x=5, y=5, z=5,
+                                            boxSizeBool=False,
+                                            boxSize=32)
+        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"))
+
+        protShiftVolPad = self.newProtocol(XmippProtShiftVolume,
+                                           inputVol=protImportVol.outputVolume,
+                                           shiftBool=False,
+                                           x=5, y=5, z=5,
+                                           boxSizeBool=False,
+                                           boxSize=80)
+        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"))
+
+
+class TestXmippProjSubtraction(TestXmippBase):
+
+    @classmethod
+    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 "
+
+        protCreatePhantom2items = self.newProtocol(XmippProtPhantom,
+                                               desc='80 80 80 0\nsph + 1 15 15 0 10\nsph + 5 -15 -15 0 10',
+                                               sampling=1.0)
+        self.launchProtocol(protCreatePhantom2items)
+        self.assertIsNotNone(protCreatePhantom2items.getFiles(),
+                             "There was a problem with phantom with 2 items creation")
+        protCreateGallery = self.newProtocol(XmippProtCreateGallery,
+                                             inputVolume=protCreatePhantom2items.outputVolume,
+                                             rotStep=15.0,
+                                             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")
+        protCreatePhantom1item = self.newProtocol(XmippProtPhantom,
+                                                  desc='80 80 80 0\nsph + 1 -15 -15 0 10',
+                                                  sampling=1.0)
+        self.launchProtocol(protCreatePhantom1item)
+        self.assertIsNotNone(protCreatePhantom1item.getFiles(),
+                             "There was a problem with phantom with 1 item creation")
+        protCreateMaskKeep = self.newProtocol(XmippProtCreateMask3D,
+                                              inputVolume=protCreatePhantom1item.outputVolume,
+                                              threshold=0.1)
+        self.launchProtocol(protCreateMaskKeep)
+        self.assertIsNotNone(protCreateMaskKeep.getFiles(),
+                             "There was a problem with the 3D mask of the 1 item phantom")
+        protSubtractProj = self.newProtocol(XmippProtSubtractProjection,
+                                            particles=protCreateGallery.outputReprojections,
+                                            vol=protCreatePhantom2items.outputVolume,
+                                            maskVol=protCreateMask.outputMask,
+                                            mask=protCreateMaskKeep.outputMask,
+                                            saveFiles=True)
+        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.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)
+
+        protSimulateCTF = self.newProtocol(XmippProtSimulateCTF,
+                                           inputParticles=protCreateGallery.outputReprojections)
+        self.launchProtocol(protSimulateCTF)
+        self.assertIsNotNone(protSimulateCTF.outputParticles,
+                             "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)
+        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.getFirstItem().getDim(), (80, 80, 1),
+                         (ERROR_DIM, "particles"))
+        self.assertEqual(protSubtractProjCTF.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+
+        protAddNoise = self.newProtocol(XmippProtAddNoiseParticles,
+                                        input=protCreateGallery.outputReprojections,
+                                        gaussianStd=50.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)
+        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.getFirstItem().getDim(), (80, 80, 1),
+                         (ERROR_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoise.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+        protAddNoiseCTF = self.newProtocol(XmippProtAddNoiseParticles,
+                                           input=protSimulateCTF.outputParticles,
+                                           gaussianStd=50.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)
+        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.getFirstItem().getDim(), (80, 80, 1),
+                         (ERROR_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTF.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+        protCreatePhantom2Over = self.newProtocol(XmippProtPhantom,
+                                                  desc='80 80 80 0\nsph + 1 5 5 0 10\nsph + 5 -5 -5 0 10',
+                                                  sampling=1.0)
+        self.launchProtocol(protCreatePhantom2Over)
+        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)
+        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")
+        protCreatePhantom1Over = self.newProtocol(XmippProtPhantom,
+                                                  desc='80 80 80 0\nsph + 1 -5 -5 0 10',
+                                                  sampling=1.0)
+        self.launchProtocol(protCreatePhantom1Over)
+        self.assertIsNotNone(protCreatePhantom1Over.getFiles(),
+                             "There was a problem with phantom with 1 item overlap creation")
+        protCreateMaskKeepOver = self.newProtocol(XmippProtCreateMask3D,
+                                                  inputVolume=protCreatePhantom1Over.outputVolume,
+                                                  threshold=0.1)
+        self.launchProtocol(protCreateMaskKeepOver)
+        self.assertIsNotNone(protCreateMaskKeepOver.getFiles(),
+                             "There was a problem with the 3D mask of the 1 item overlap phantom")
+        protSubtractProjOver = self.newProtocol(XmippProtSubtractProjection,
+                                                particles=protCreateGalleryOver.outputReprojections,
+                                                vol=protCreatePhantom2Over.outputVolume,
+                                                maskVol=protCreateMaskOver.outputMask,
+                                                mask=protCreateMaskKeepOver.outputMask,
+                                                saveFiles=True)
+        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.getFirstItem().getDim(), (80, 80, 1),
+                         (ERROR_DIM, "particles"))
+        self.assertEqual(protSubtractProjOver.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+        protSimulateCTFOver = self.newProtocol(XmippProtSimulateCTF,
+                                               inputParticles=protCreateGalleryOver.outputReprojections)
+        self.launchProtocol(protSimulateCTFOver)
+        self.assertIsNotNone(protSimulateCTFOver.outputParticles,
+                             "There was a problem with overlap CTF simulation")
+        protAddNoiseCTFOver = self.newProtocol(XmippProtAddNoiseParticles,
+                                               input=protSimulateCTFOver.outputParticles,
+                                               gaussianStd=50.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)
+        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.getFirstItem().getDim(), (80, 80, 1),
+                         (ERROR_DIM, "particles"))
+        self.assertEqual(protSubtractProjNoiseCTFOver.outputParticles.getSize(), 181, ERROR_SIZE_PART)
+
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        className = sys.argv[1]
+        cls = globals().get(className, None)
+        if cls:
+            suite = unittest.TestLoader().loadTestsFromTestCase(cls)
+            unittest.TextTestRunner(verbosity=2).run(suite)
+        else:
+            print("Test: '%s' not found." % className)
+    else:
+        unittest.main()
diff --git a/xmipp3/tests/test_protocols_xmipp_mics.py b/xmipp3/tests/test_protocols_xmipp_mics.py
index 8bd4eaad8..7f9e95b84 100644
--- a/xmipp3/tests/test_protocols_xmipp_mics.py
+++ b/xmipp3/tests/test_protocols_xmipp_mics.py
@@ -990,7 +990,7 @@ def testStreamingAndNonStreaming(self):
         self.launchProtocol(protCons1)
 
         self.assertTrue(protCons1.isFinished(), "Consensus failed")
-        self.assertSetSize(protCons1.consensusCoordinates,383,
+        self.assertSetSize(protCons1.consensusCoordinates,382,
                         "Output coordinates size for AND consensus is wrong.")
 
         protConsOr = self.newProtocol(XmippProtConsensusPicking,
@@ -1039,7 +1039,7 @@ def testStreamingAndNonStreaming(self):
 
         time.sleep(3)  # protDupl2 should be as long as protCons2, but just in case
         protCons2 = self._updateProtocol(protCons2)
-        self.assertSetSize(protCons2.consensusCoordinates, 383,
+        self.assertSetSize(protCons2.consensusCoordinates, 382,
                         "Output coordinates size does not is wrong.")
         protDupl2 = self._updateProtocol(protDupl2)
         self.assertSetSize(protDupl2.outputCoordinates, 245,
diff --git a/xmipp3/tests/test_protocols_zernike3d.py b/xmipp3/tests/test_protocols_zernike3d.py
index 2ffca6eba..c3ba43a66 100644
--- a/xmipp3/tests/test_protocols_zernike3d.py
+++ b/xmipp3/tests/test_protocols_zernike3d.py
@@ -145,8 +145,8 @@ def testPairAlignmentZernike3DCPU(self):
         self.assertEqual(clnm_gold[0], clnm_test[0], "There is a problem with Zernike degree")
         self.assertEqual(clnm_gold[1], clnm_test[1], "There is a problem with Spherical Harmonics degree")
         self.assertEqual(clnm_gold[2], clnm_test[2], "There is a problem with sphere radius")
-        diff_clnm = np.abs(np.sum(clnm_gold[3] - clnm_test[3]))
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+        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):
@@ -169,8 +169,8 @@ def testPairAlignmentZernike3DGPU(self):
         self.assertEqual(clnm_gold[0], clnm_test[0], "There is a problem with Zernike degree")
         self.assertEqual(clnm_gold[1], clnm_test[1], "There is a problem with Spherical Harmonics degree")
         self.assertEqual(clnm_gold[2], clnm_test[2], "There is a problem with sphere radius")
-        diff_clnm = np.abs(np.sum(clnm_gold[3] - clnm_test[3]))
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+        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):
@@ -190,22 +190,22 @@ def testStructMapZernike3DCPU(self):
                                      delimiter=' ')
         sm_test_def = np.loadtxt(structureMap._getExtraPath('CoordinateMatrix3.txt'),
                                      delimiter=' ')
-        diff_sm_def = np.abs(np.sum(sm_gold_def - sm_test_def))
-        self.assertAlmostEqual(diff_sm_def, 0, delta=0.2,
+        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.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+        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.abs(np.sum(sm_cpu_gold_cons - sm_cpu_test_cons))
-        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.2,
+        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):
@@ -222,8 +222,8 @@ def testStructMapCPU(self):
                                      delimiter=' ')
         sm_cpu_test_corr = np.loadtxt(structureMap._getExtraPath(self.CORR_FILE),
                                      delimiter=' ')
-        diff_sm_corr = np.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+        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)
 
 
@@ -244,22 +244,22 @@ def testStructMapZernike3DGPU(self):
                                      delimiter=' ')
         sm_test_def = np.loadtxt(structureMap._getExtraPath('CoordinateMatrix3.txt'),
                                      delimiter=' ')
-        diff_sm_def = np.abs(np.sum(sm_gold_def - sm_test_def))
-        self.assertAlmostEqual(diff_sm_def, 0, delta=0.2,
+        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.abs(np.sum(sm_cpu_gold_corr - sm_cpu_test_corr))
-        self.assertAlmostEqual(diff_sm_corr, 0, delta=0.2,
+        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.abs(np.sum(sm_cpu_gold_cons - sm_cpu_test_cons))
-        self.assertAlmostEqual(diff_sm_cons, 0, delta=0.2,
+        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):
@@ -279,8 +279,8 @@ def testAngularAlignZernike3DCPU(self):
         # Check if coefficients are fine
         aa_gold_clnm = self.readMetadata(self.dataset.getFile('gold_standard_aa/CPU/output_particles.xmd'))
         aa_test_clnm = self.readMetadata(alignZernike._getExtraPath('output_particles.xmd'))
-        diff_clnm = np.abs(np.sum(aa_gold_clnm - aa_test_clnm))
-        self.assertAlmostEqual(diff_clnm, 0, delta=0.2,
+        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):
@@ -300,6 +300,6 @@ def testAngularAlignZernike3DGPU(self):
         # Check if coefficients are fine
         aa_gold_clnm = self.readMetadata(self.dataset.getFile('gold_standard_aa/GPU/output_particles.xmd'))
         aa_test_clnm = self.readMetadata(alignZernike._getExtraPath('output_particles.xmd'))
-        diff_clnm = np.abs(np.sum(aa_gold_clnm - aa_test_clnm))
-        self.assertAlmostEqual(diff_clnm, 0, delta=4,
+        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