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

Issue:
Center particles: protocol simplified
Before the protocol to center particles had as output the centered classes, the centered particle and the coordinates.

With this PR the protocol is simplified and only centers the particle returning the classes and the particles, if the set of coordinates is needed,  you have to manually add the protocol to extract coordinates.


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 573543fb4..e57cc8df7 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,5 @@
+[![Build](https://github.com/I2PC/scipion-em-xmipp/actions/workflows/build.yml/badge.svg)](https://github.com/I2PC/scipion-em-xmipp/actions/workflows/build.yml)
+
 # scipion-em-xmipp
 
 Plugin to use Xmipp programs within the Scipion framework.
diff --git a/setup.py b/setup.py
index d06a58c61..028b96d3a 100644
--- a/setup.py
+++ b/setup.py
@@ -46,7 +46,7 @@
     # For a discussion on single-sourcing the version across setup.py and the
     # project code, see
     # https://packaging.python.org/en/latest/single_source_version.html
-    version='22.01.0',  # Required
+    version='22.04.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 cf45b64e4..d9d368c5e 100644
--- a/xmipp3/__init__.py
+++ b/xmipp3/__init__.py
@@ -39,7 +39,7 @@
 
 _logo = "xmipp_logo.png"
 _references = ['delaRosaTrevin2013', 'Sorzano2013']
-_currentVersion = '3.22.01.0'
+_currentVersion = '3.22.04.0'
 
 
 class Plugin(pwem.Plugin):
diff --git a/xmipp3/protocols.conf b/xmipp3/protocols.conf
index 07946b2aa..3239c5b24 100644
--- a/xmipp3/protocols.conf
+++ b/xmipp3/protocols.conf
@@ -63,6 +63,10 @@ Protocols SPA = [
 		]},
 		{"tag": "protocol_group", "text": "Classify", "openItem": "False", "children": [
 			{"tag": "protocol", "value": "XmippProtCL2D",   "text": "default"},			
+			{"tag": "protocol", "value": "XmippProtCoreAnalysis",   "text": "default"},
+			{"tag": "protocol", "value": "XmippProtGpuCrrCL2D",   "text": "default"},
+			{"tag": "protocol", "value": "XmippProtStrGpuCrrCL2D",   "text": "default"},
+			{"tag": "protocol", "value": "XmippProtStrGpuCrrSimple",   "text": "default"},
 			{"tag": "section", "text": "more", "openItem": "False", "children": [
 				{"tag": "protocol", "value": "XmippProtML2D", "text": "default"},
 				{"tag": "protocol", "value": "XmippProtKerdensom", "text": "default"},
@@ -121,6 +125,7 @@ Protocols SPA = [
 		{"tag": "protocol", "value": "XmippProtProjectZ", "text": "default"},
 		{"tag": "protocol", "value": "XmippProtValidateNonTilt",    "text": "default"},
 		{"tag": "protocol", "value": "XmippProtMultiRefAlignability", "text": "default"},
+		{"tag": "protocol", "value": "XmippProtAngularGraphConsistency", "text": "default"},
 		{"tag": "protocol", "value": "XmippProtValidateOverfitting",    "text": "default"}
 	]},
 	{"tag": "section", "text": "Resolution", "openItem": "False", "children": [
diff --git a/xmipp3/protocols/__init__.py b/xmipp3/protocols/__init__.py
index 25fa08b60..947db4f2d 100644
--- a/xmipp3/protocols/__init__.py
+++ b/xmipp3/protocols/__init__.py
@@ -30,6 +30,7 @@
 
 from .protocol_assignment_tilt_pair import XmippProtAssignmentTiltPair
 from .protocol_align_volume import XmippProtAlignVolume, XmippProtAlignVolumeForWeb
+from .protocol_angular_graph_consistency import XmippProtAngularGraphConsistency
 from .protocol_preprocess.protocol_add_noise import (XmippProtAddNoiseVolumes,
                                                      XmippProtAddNoiseParticles)
 from .protocol_apply_alignment import XmippProtApplyAlignment
@@ -45,6 +46,7 @@
 from .protocol_compare_reprojections import XmippProtCompareReprojections
 from .protocol_compare_angles import XmippProtCompareAngles
 from .protocol_convert_pdb import XmippProtConvertPdb
+from .protocol_core_analysis import XmippProtCoreAnalysis
 from .protocol_create_gallery import XmippProtCreateGallery
 from .protocol_ctf_consensus import XmippProtCTFConsensus
 from .protocol_ctf_micrographs import XmippProtCTFMicrographs
diff --git a/xmipp3/protocols/protocol_angular_graph_consistency.py b/xmipp3/protocols/protocol_angular_graph_consistency.py
new file mode 100644
index 000000000..9d2020961
--- /dev/null
+++ b/xmipp3/protocols/protocol_angular_graph_consistency.py
@@ -0,0 +1,568 @@
+# **************************************************************************
+# *
+# * Authors:         Jeison Mendez (jmendez@utp.edu.co) (2020)
+# *
+# * Universidad Nacional Autónoma de México, UNAM
+# *
+# * 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 glob import glob
+import math
+from os.path import join, isfile, exists
+from shutil import copyfile, copy
+import os
+
+from pyworkflow.protocol.params import (PointerParam, FloatParam,
+                                        STEPS_PARALLEL,
+                                        StringParam, BooleanParam, IntParam,
+                                        LEVEL_ADVANCED, USE_GPU, GPU_LIST)
+
+from pyworkflow.protocol import String, Float
+
+from pyworkflow.protocol import params
+
+from pyworkflow.utils.path import (moveFile, makePath, cleanPattern, copyFile, cleanPath)
+from pyworkflow.gui.plotter import Plotter
+
+from pwem.objects import Volume, SetOfVolumes
+from pwem import emlib
+from pwem.emlib.image import ImageHandler
+from pwem.emlib.metadata import getFirstRow, getSize, iterRows
+import pwem.emlib.metadata as md
+from pwem.protocols import ProtAnalysis3D
+from pwem.constants import ALIGN_PROJ
+
+from xmipp3.convert import writeSetOfParticles, writeSetOfVolumes, \
+    getImageLocation, setXmippAttributes, createItemMatrix, readSetOfParticles, readSetOfImages
+from _ast import Try
+
+
+class XmippProtAngularGraphConsistency(ProtAnalysis3D):
+    """
+    Performs soft alignment validation of a set of particles previously aligned
+    confronting them using Graph filtered correlations representation. This
+    protocol produces an histogram with two groups of particles.
+    """
+    _label = 'angular graph consistency'
+    
+    def __init__(self, *args, **kwargs):
+        ProtAnalysis3D.__init__(self, *args, **kwargs)
+
+    # --------------------------- DEFINE param functions --------------------------------------------
+    def _defineParams(self, form):
+
+        form.addSection(label='Input')
+        form.addParam('inputVolume', PointerParam, pointerClass='Volume',
+                      label="Input volume", important=True,
+                      help='Select the input volume(s).')
+        form.addParam('inputParticles', PointerParam,
+                      pointerClass='SetOfParticles',
+                      pointerCondition='hasAlignment',
+                      label="Input particles", important=True,
+                      help='Select the input projection images.')
+        form.addParam('symmetryGroup', StringParam, default='c1',
+                      label="Symmetry group",
+                      help='See [[Xmipp Symmetry][http://www2.mrc-lmb.cam.ac.uk/Xmipp/index.php/Conventions_%26_File_formats#Symmetry]] page '
+                           'for a description of the symmetry format accepted by Xmipp')
+        form.addParam('angularSampling', FloatParam, default=3,
+                      expertLevel=LEVEL_ADVANCED,
+                      label="Angular Sampling (degrees)",
+                      help='Angular distance (in degrees) between neighboring projection points ')
+        form.addParam('ccLevel', FloatParam, default=0.95,
+                      expertLevel=LEVEL_ADVANCED,
+                      label="correlation level for validation",
+                      help='threshold correlation to be used in validation. Keep this value in the range [0.90 -- 0.97]')
+        form.addParam('minTilt', FloatParam, default=0,
+                      expertLevel=LEVEL_ADVANCED,
+                      label="Minimum allowed tilt angle",
+                      help='Tilts below this value will not be considered for the alignment')
+        form.addParam('maxTilt', FloatParam, default=180,
+                      expertLevel=LEVEL_ADVANCED,
+                      label="Maximum allowed tilt angle without mirror check",
+                      help='Tilts above this value will not be considered for the alignment without mirror check')
+        form.addParam('maximumTargetResolution', FloatParam, default=8,
+                      label='Target resolution (A)',
+                      help='Low pass filter the particles to this resolution. This usually helps a lot obtaining good alignment. You should have a good'
+                           ' reason to modify this value outside the range  [8-10] A')
+
+        form.addParallelSection(threads=1, mpi=1)
+
+    # --------------------------- INSERT steps functions --------------------------------------------
+
+    def _insertAllSteps(self):
+        self.imgsFn = self._getExtraPath('images.xmd')
+        self._insertFunctionStep('convertInputStep',
+                                             self.inputParticles.get().getObjId())
+        self.TsOrig = self.inputParticles.get().getSamplingRate()
+        self._insertFunctionStep('doIteration000')
+        self._maximumTargetResolution = self.maximumTargetResolution.get()
+        self._insertFunctionStep('globalAssignment')  
+        self._insertFunctionStep('cleanDirectory') 
+        self._insertFunctionStep("createOutput")        
+
+    def convertInputStep(self, inputParticlesId):
+        writeSetOfParticles(self.inputParticles.get(), self.imgsFn)
+        self.runJob('xmipp_metadata_utilities', '-i %s --fill image1 constant noImage' % self.imgsFn, numberOfMpi=1)
+        self.runJob('xmipp_metadata_utilities', '-i %s --operate modify_values "image1=image"' % self.imgsFn, numberOfMpi=1)
+        self.runJob('xmipp_metadata_utilities', '-i %s --fill particleId constant 1' % self.imgsFn, numberOfMpi=1)
+        self.runJob('xmipp_metadata_utilities', '-i %s --operate modify_values "particleId=itemId"' % self.imgsFn, numberOfMpi=1)
+        imgsFnId = self._getExtraPath('imagesId.xmd')
+        self.runJob('xmipp_metadata_utilities', '-i %s --operate keep_column particleId -o %s' % (self.imgsFn, imgsFnId), numberOfMpi=1)
+
+    def createOutput(self):
+        fnLastDir = self._getExtraPath("Iter1") 
+        fnDirGlobal = fnLastDir
+        Ts = self.readInfoField(fnDirGlobal, "sampling", emlib.MDL_SAMPLINGRATE)
+
+        fnAnglesDisc=join(fnDirGlobal,"anglesDisc.xmd")
+        if exists(fnAnglesDisc):
+            fnAngles = self._getPath("angles.xmd")
+            self.runJob('xmipp_metadata_utilities','-i %s -o %s --operate modify_values "image=image1"'%(fnAnglesDisc,fnAngles), numberOfMpi=1)
+            self.runJob('xmipp_metadata_utilities','-i %s --operate sort particleId'%fnAngles,numberOfMpi=1)
+            self.runJob('xmipp_metadata_utilities','-i %s --operate drop_column image1'%fnAngles,numberOfMpi=1)
+            self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "itemId=particleId"'%fnAngles,numberOfMpi=1)
+            imgSet=self.inputParticles.get()
+            self.scaleFactor = Ts/imgSet.getSamplingRate()
+            outImgSet= self._createSetOfParticles()
+            outImgSet.copyInfo(imgSet)
+            outImgSet.setAlignmentProj()
+            outImgSet.setIsPhaseFlipped(imgSet.isPhaseFlipped())
+            self.iterMd = md.iterRows(fnAngles, md.MDL_PARTICLE_ID)
+            self.lastRow = next(self.iterMd)
+            outImgSet.copyItems(imgSet,
+                                updateItemCallback=self._updateItem)
+            self._defineOutputs(outputParticles=outImgSet)
+            self._defineSourceRelation(self.inputParticles, outImgSet)
+            
+            mdParticles = emlib.MetaData(fnAngles)
+            ccGraphPrevList = mdParticles.getColumnValues(emlib.MDL_GRAPH_CC_PREVIOUS)
+            ccAsigDirRefList = mdParticles.getColumnValues(emlib.MDL_ASSIGNED_DIR_REF_CC)
+            angDistPrevList = mdParticles.getColumnValues(emlib.MDL_GRAPH_DISTANCE2MAX_PREVIOUS)
+            
+            # threshold            
+            th_ccGraph = self.otsu(ccGraphPrevList)
+            th_ccAsigDir = self.otsu(ccAsigDirRefList)
+            th_angDist = self.otsu(angDistPrevList)
+            if (th_ccGraph < 0.50) and (th_ccAsigDir < 0.50):
+                print("applying default thresholds")
+                th_ccGraph = 0.95
+                th_ccAsigDir = 0.95
+                th_angDist = 3.0 * self.angularSampling.get() 
+                print(th_ccGraph, th_ccAsigDir, th_angDist)
+            
+            print("\033[1;33Thresholds:\033[0;33 \n")
+            print('correlation with projection in Graph max direction:',th_ccGraph,
+                  '\ncorrelation with assigned projection:',th_ccAsigDir,
+                  '\nangular distance to maxGraph:',th_angDist)
+            
+            # a second output set
+            fnOutParticles = self._getPath('anglesAux.xmd')
+            copy(fnAngles, fnOutParticles)
+            mdParticles = emlib.MetaData(fnOutParticles)
+            n_false = 0;
+            
+            nParticles = 0
+            for row in iterRows(mdParticles):
+                nParticles += 1
+                objId = row.getObjId()
+                assigDirRefCC = row.getValue(emlib.MDL_ASSIGNED_DIR_REF_CC)
+                graphCCPrev = row.getValue(emlib.MDL_GRAPH_CC_PREVIOUS)
+                graphDistMaxGraphPrev = row.getValue(emlib.MDL_GRAPH_DISTANCE2MAX_PREVIOUS)
+
+                if (assigDirRefCC < th_ccAsigDir) and ( graphCCPrev < th_ccGraph ) and (graphDistMaxGraphPrev > th_angDist):
+                    n_false += 1
+                    mdParticles.setValue(emlib.MDL_ENABLED, -1, objId)
+            mdParticles.write(fnOutParticles) 
+            print('to be disabled:',n_false)          
+            self.subsets = []
+            i=0
+            self.subsets.append(self._createSetOfParticles(str(i)))
+            self.subsets[i].copyInfo(self.inputParticles.get())
+            readSetOfParticles(fnOutParticles, self.subsets[i])
+            result = {'outputParticlesAux' : self.subsets[i]}
+            self._defineOutputs(**result)
+            self._store(self.subsets[i])
+
+            # extra output info
+            p_gsp = (1 - n_false/nParticles)* 100
+            self.percentage = Float(p_gsp)
+            self._store()
+            
+    def _updateItem(self, particle, row):
+        count = 0
+        while self.lastRow and particle.getObjId() == self.lastRow.getValue(md.MDL_PARTICLE_ID):
+            count += 1
+            if count:
+                self._createItemMatrix(particle, self.lastRow)
+            try:
+                self.lastRow = next(self.iterMd)
+            except StopIteration:
+                self.lastRow = None
+
+        particle._appendItem = count > 0
+        
+    def _createItemMatrix(self, particle, row):
+        row.setValue(emlib.MDL_SHIFT_X, row.getValue(emlib.MDL_SHIFT_X)*self.scaleFactor)
+        row.setValue(emlib.MDL_SHIFT_Y, row.getValue(emlib.MDL_SHIFT_Y)*self.scaleFactor)
+        
+        setXmippAttributes(particle, row, md.MDL_ANGLE_ROT, md.MDL_ANGLE_TILT,
+                           md.MDL_ANGLE_PSI, md.MDL_SHIFT_X, md.MDL_SHIFT_Y,
+                           md.MDL_MAXCC, md.MDL_WEIGHT,
+                           md.MDL_MAXCC_PREVIOUS, md.MDL_GRAPH_CC_PREVIOUS,
+                           md.MDL_ASSIGNED_DIR_REF_CC, md.MDL_GRAPH_DISTANCE2MAX_PREVIOUS,
+                           md.MDL_GRAPH_CC, md.MDL_GRAPH_DISTANCE2MAX)
+        
+        createItemMatrix(particle, row, align=ALIGN_PROJ)
+
+    def doIteration000(self):
+        fnDirCurrent = self._getExtraPath('Iter0')
+        makePath(fnDirCurrent)
+
+        # Get volume sampling rate
+        if self.inputVolume.get() is None:
+            TsCurrent = self.inputParticles.get().getSamplingRate()
+        else:
+            TsCurrent = self.inputVolume.get().getSamplingRate()
+        self.writeInfoField(fnDirCurrent, "sampling", emlib.MDL_SAMPLINGRATE, TsCurrent)
+
+        # Copy reference volume and window if necessary
+        XDim = self.inputParticles.get().getDimensions()[0]
+        newXDim = int(round(XDim * self.TsOrig / TsCurrent))
+        self.writeInfoField(fnDirCurrent, "size", emlib.MDL_XSIZE, newXDim)
+        
+        img = ImageHandler()
+
+        fnVol1 = join(fnDirCurrent, "volume.vol")  
+
+        vol = self.inputVolume.get()
+        img.convert(vol, fnVol1)
+        volXdim = vol.getDim()[0]
+        if newXDim != volXdim:
+            self.runJob('xmipp_transform_window', "-i %s --size %d" % (fnVol1, newXDim), numberOfMpi=1)
+
+    def globalAssignment(self):  
+        fnDirIter0 = self._getExtraPath("Iter0") 
+        fnGlobal = self._getExtraPath("Iter1")
+        makePath(fnGlobal)
+        previousResolution = self._maximumTargetResolution
+        targetResolution = self._maximumTargetResolution
+        TsCurrent = max(self.TsOrig,targetResolution/3)
+        getShiftsFrom = ''
+        self.prepareImages(fnDirIter0, fnGlobal, TsCurrent, getShiftsFrom)
+        TsCurrent = self.readInfoField(fnGlobal, "sampling", emlib.MDL_SAMPLINGRATE)  # Prepare images may have changed it
+        self.prepareReferences(fnDirIter0, fnGlobal, TsCurrent, targetResolution)
+        
+        # Calculate angular step at this resolution
+        ResolutionAlignment = previousResolution
+        self.nextLowPass = True
+        self.nextResolutionOffset = 2
+        if self.nextLowPass:
+            ResolutionAlignment += self.nextResolutionOffset
+        newXDim = self.readInfoField(fnGlobal, "size", emlib.MDL_XSIZE)
+        angleStep = self.calculateAngStep(newXDim, TsCurrent, ResolutionAlignment)
+        angleStep = max(angleStep, self.angularSampling.get())
+        self.writeInfoField(fnGlobal, "angleStep", emlib.MDL_ANGLE_DIFF, float(angleStep))
+
+        # Global alignment
+        self.numberOfPerturbations = 1
+        perturbationList = [chr(x) for x in range(ord('a'), ord('a') + self.numberOfPerturbations)]
+
+        fnDirAssignment = join(fnGlobal, "assignment")
+        fnImgs = join(fnGlobal, "images.xmd")
+        makePath(fnDirAssignment)
+
+        # Create defocus groups
+        row = getFirstRow(fnImgs)
+        if row.containsLabel(emlib.MDL_CTF_MODEL) or row.containsLabel(emlib.MDL_CTF_DEFOCUSU):
+            self.runJob("xmipp_ctf_group", "--ctfdat %s -o %s/ctf:stk --pad 1.0 --sampling_rate %f --phase_flipped  --error 0.1 --resol %f" % \
+                        (fnImgs, fnDirAssignment, TsCurrent, targetResolution), numberOfMpi=1)
+            moveFile("%s/ctf_images.sel" % fnDirAssignment, "%s/ctf_groups.xmd" % fnDirAssignment)
+            cleanPath("%s/ctf_split.doc" % fnDirAssignment)
+            mdInfo = emlib.MetaData("numberGroups@%s" % join(fnDirAssignment, "ctfInfo.xmd"))
+            numberGroups = mdInfo.getValue(emlib.MDL_COUNT, mdInfo.firstObject())
+            ctfPresent = True
+        else:
+            numberGroups = 1
+            ctfPresent = False
+            
+        # Generate projections
+        fnReferenceVol = join(fnGlobal, "volumeRef.vol") 
+        for subset in perturbationList:
+            fnGallery = join(fnDirAssignment, "gallery%s.stk" % subset)
+            fnGalleryMd = join(fnDirAssignment, "gallery%s.xmd" % subset)
+            args = "-i %s -o %s --sampling_rate %f --sym %s --min_tilt_angle %f --max_tilt_angle %f --perturb %f " % \
+                   (fnReferenceVol, fnGallery, angleStep, self.symmetryGroup, self.minTilt.get(), self.maxTilt.get(), math.sin(angleStep * math.pi / 180.0) / 4)
+            args += " --compute_neighbors --angular_distance -1 --experimental_images %s" % self._getExtraPath("images.xmd")
+            self.runJob("xmipp_angular_project_library", args, numberOfMpi = self.numberOfMpi.get() * self.numberOfThreads.get())
+            cleanPath(join(fnDirAssignment, "gallery_angles%s.doc" % subset))
+            moveFile(join(fnDirAssignment, "gallery%s.doc" % subset), fnGalleryMd)
+            fnAngles = join(fnGlobal, "anglesDisc%s.xmd" % subset)
+            for j in range(1, numberGroups + 1):
+                fnAnglesGroup = join(fnDirAssignment, "angles_group%03d%s.xmd" % (j, subset))
+                if not exists(fnAnglesGroup):
+                    if ctfPresent:
+                        fnGroup = "ctfGroup%06d@%s/ctf_groups.xmd" % (j, fnDirAssignment)
+                        fnCTF = "%d@%s/ctf_ctf.stk" % (j, fnDirAssignment)
+                        fnGalleryGroup = join(fnDirAssignment, "gallery%s_%06d.stk" % (subset, j))
+                        fnGalleryGroupMd = join(fnDirAssignment, "gallery%s_%06d.xmd" % (subset, j))
+                        self.runJob("xmipp_transform_filter", "-i %s -o %s --fourier binary_file %s --save_metadata_stack %s --keep_input_columns" % \
+                                    (fnGalleryMd, fnGalleryGroup, fnCTF, fnGalleryGroupMd),
+                                    numberOfMpi=min(self.numberOfMpi.get(), 24))
+                    else:
+                        fnGroup = fnImgs
+                        fnGalleryGroupMd = fnGalleryMd
+                    if getSize(fnGroup) == 0:  # If the group is empty
+                        continue
+                    self.angularMaxShift = 10
+                    maxShift = round(self.angularMaxShift * newXDim / 100)
+                    R = self.inputParticles.get().getDimensions()[0] / 2
+                    R = R * self.TsOrig / TsCurrent    
+                    args = '-i %s -o %s -ref %s -sampling %f -odir %s --Nsimultaneous %d -angleStep %f --maxShift %f --sym %s --useForValidation --refVol %s' % \
+                        (fnGroup, fnAnglesGroup, fnGalleryGroupMd, TsCurrent, fnDirAssignment, self.numberOfMpi.get() * self.numberOfThreads.get(), angleStep, maxShift, self.symmetryGroup, fnReferenceVol)
+                    self.runJob('xmipp_angular_assignment_mag', args, numberOfMpi=self.numberOfMpi.get() * self.numberOfThreads.get())
+                    if exists(fnAnglesGroup):
+                        if not exists(fnAngles) and exists(fnAnglesGroup):
+                            copyFile(fnAnglesGroup, fnAngles)
+                        else:
+                            if exists(fnAngles) and exists(fnAnglesGroup):
+                                self.runJob("xmipp_metadata_utilities", "-i %s --set union_all %s" % (fnAngles, fnAnglesGroup), numberOfMpi=1)
+            if exists(fnAngles) and exists(fnImgs):
+                self.runJob("xmipp_metadata_utilities", "-i %s --set join %s image" % (fnAngles, fnImgs), numberOfMpi=1)
+            self.saveSpace = True
+            if self.saveSpace and ctfPresent:
+                self.runJob("rm -f", fnDirAssignment + "/gallery*", numberOfMpi=1)  
+                
+        # Evaluate the stability of the alignment
+        fnOut = join(fnGlobal, "anglesDisc")
+        for set1 in perturbationList:
+            fnOut1 = join(fnGlobal, "anglesDisc%s" % set1)
+            fnAngles1 = fnOut1 + ".xmd"
+            counter2 = 0
+            for set2 in perturbationList:
+                if set1 == set2:
+                    continue
+                fnAngles2 = join(fnGlobal, "anglesDisc%s.xmd" % set2)
+                fnOut12 = join(fnGlobal, "anglesDisc%s%s" % (set1, set2))
+                self.runJob("xmipp_angular_distance", "--ang1 %s --ang2 %s --oroot %s --sym %s --compute_weights 1 particleId 0.5 --check_mirrors --set 0" % (fnAngles2, fnAngles1, fnOut12, self.symmetryGroup), numberOfMpi=1)
+                self.runJob("xmipp_metadata_utilities", '-i %s --operate keep_column "angleDiff0 shiftDiff0 weightJumper0"' % (fnOut12 + "_weights.xmd"), numberOfMpi=1)
+                if counter2 == 0:
+                    mdAllWeights = emlib.MetaData(fnOut12 + "_weights.xmd")
+                    counter2 = 1
+                else:
+                    mdWeight = emlib.MetaData(fnOut12 + "_weights.xmd")
+                    if mdWeight.size() == mdAllWeights.size():
+                        counter2 += 1
+                        for id1, id2 in izip(mdWeight, mdAllWeights):
+                            angleDiff0 = mdWeight.getValue(emlib.MDL_ANGLE_DIFF0, id1)
+                            shiftDiff0 = mdWeight.getValue(emlib.MDL_SHIFT_DIFF0, id1)
+                            weightJumper0 = mdWeight.getValue(emlib.MDL_WEIGHT_JUMPER0, id1)
+
+                            angleDiff0All = mdAllWeights.getValue(emlib.MDL_ANGLE_DIFF0, id2)
+                            shiftDiff0All = mdAllWeights.getValue(emlib.MDL_SHIFT_DIFF0, id2)
+                            weightJumper0All = mdAllWeights.getValue(emlib.MDL_WEIGHT_JUMPER0, id2)
+
+                            mdAllWeights.setValue(emlib.MDL_ANGLE_DIFF0, angleDiff0 + angleDiff0All, id2)
+                            mdAllWeights.setValue(emlib.MDL_SHIFT_DIFF0, shiftDiff0 + shiftDiff0All, id2)
+                            mdAllWeights.setValue(emlib.MDL_WEIGHT_JUMPER0, weightJumper0 + weightJumper0All, id2)
+            if counter2 > 1:
+                iCounter2 = 1.0 / counter2
+                for id in mdAllWeights:
+                    angleDiff0All = mdAllWeights.getValue(emlib.MDL_ANGLE_DIFF0, id)
+                    shiftDiff0All = mdAllWeights.getValue(emlib.MDL_SHIFT_DIFF0, id)
+                    weightJumper0All = mdAllWeights.getValue(emlib.MDL_WEIGHT_JUMPER0, id)
+
+                    mdAllWeights.setValue(emlib.MDL_ANGLE_DIFF0, angleDiff0All * iCounter2, id)
+                    mdAllWeights.setValue(emlib.MDL_SHIFT_DIFF0, shiftDiff0All * iCounter2, id)
+                    mdAllWeights.setValue(emlib.MDL_WEIGHT_JUMPER0, weightJumper0All * iCounter2, id)
+            if counter2 > 0:
+                mdAllWeights.write(fnOut1 + "_weights.xmd")
+                self.runJob("xmipp_metadata_utilities", '-i %s --set merge %s' % (fnAngles1, fnOut1 + "_weights.xmd"), numberOfMpi=1)
+            if not exists(fnOut + ".xmd") and exists(fnAngles1):
+                copyFile(fnAngles1, fnOut + ".xmd")
+            else:
+                if exists(fnAngles1) and exists(fnOut + ".xmd"):
+                    self.runJob("xmipp_metadata_utilities", '-i %s --set union_all %s' % (fnOut + ".xmd", fnAngles1), numberOfMpi=1)
+        cleanPath(join(fnGlobal, "anglesDisc*_weights.xmd"))    
+
+    def calculateAngStep(self, newXDim, TsCurrent, ResolutionAlignment):
+        k = newXDim * TsCurrent / ResolutionAlignment  # Freq. index
+        return math.atan2(1, k) * 180.0 / math.pi  # Corresponding angular step
+
+    def checkInfoField(self, fnDir, block):
+        fnInfo = join(fnDir, "iterInfo.xmd")
+        if not exists(fnInfo):
+            return False
+        blocks = emlib.getBlocksInMetaDataFile(fnInfo)
+        return block in blocks
+            
+    def writeInfoField(self, fnDir, block, label, value):
+        mdInfo = emlib.MetaData()
+        objId = mdInfo.addObject()
+        mdInfo.setValue(label, value, objId)
+        mdInfo.write("%s@%s" % (block, join(fnDir, "iterInfo.xmd")), emlib.MD_APPEND)
+        
+    def readInfoField(self, fnDir, block, label):
+        mdInfo = emlib.MetaData("%s@%s" % (block, join(fnDir, "iterInfo.xmd")))
+        return mdInfo.getValue(label, mdInfo.firstObject())    
+    
+    def prepareImages(self, fnDirPrevious, fnDir, TsCurrent, getShiftsFrom=''):
+        if self.checkInfoField(fnDir, "count"):
+            state = self.readInfoField(fnDir, "count", emlib.MDL_COUNT)
+            if state >= 1:
+                return        
+        print("Preparing images to sampling rate=", TsCurrent)
+        XDim = self.inputParticles.get().getDimensions()[0]
+        newXDim = int(round(XDim * self.TsOrig / TsCurrent))
+        if newXDim < 40:
+            newXDim = int(40)
+            TsCurrent = XDim * (self.TsOrig / newXDim)
+        elif newXDim % 2 == 1:
+            newXDim += 1
+            TsCurrent = XDim * (self.TsOrig / newXDim)
+        self.writeInfoField(fnDir, "sampling", emlib.MDL_SAMPLINGRATE, TsCurrent)
+        self.writeInfoField(fnDir, "size", emlib.MDL_XSIZE, newXDim)
+        
+        # Prepare particles
+        fnNewParticles = join(fnDir, "images.stk")
+        if newXDim != XDim:
+            self.runJob("xmipp_image_resize", "-i %s -o %s --fourier %d" % (self.imgsFn, fnNewParticles, newXDim), numberOfMpi=min(self.numberOfMpi.get(), 24))
+        else:
+            self.runJob("xmipp_image_convert", "-i %s -o %s --save_metadata_stack %s" % (self.imgsFn, fnNewParticles, join(fnDir, "images.xmd")),
+                        numberOfMpi=1)
+        R = self.inputParticles.get().getDimensions()[0] / 2
+        self.angularMaxShift = 10
+        R = min(round(R * self.TsOrig / TsCurrent * (1 + self.angularMaxShift * 0.01)), newXDim / 2)
+        self.runJob("xmipp_transform_mask", "-i %s --mask circular -%d" % (fnNewParticles, R), numberOfMpi=min(self.numberOfMpi.get(), 24))
+        fnSource = join(fnDir, "images.xmd")
+
+        if not self.inputParticles.get().isPhaseFlipped():
+            self.runJob("xmipp_ctf_correct_phase", "-i %s --sampling_rate %f" % (fnSource, TsCurrent),
+                        numberOfMpi=min(self.numberOfMpi.get(), 24))
+
+        self.writeInfoField(fnDir, "count", emlib.MDL_COUNT, int(1))   
+        
+    def prepareReferences(self, fnDirPrevious, fnDir, TsCurrent, targetResolution):
+        if self.checkInfoField(fnDir, "count"):
+            state = self.readInfoField(fnDir, "count", emlib.MDL_COUNT)
+            if state >= 2:
+                return
+
+        print("Preparing references to sampling rate=", TsCurrent)
+        fnMask = ''
+        newXDim = self.readInfoField(fnDir, "size", emlib.MDL_XSIZE)
+        oldXdim = self.readInfoField(fnDirPrevious, "size", emlib.MDL_XSIZE)
+        fnPreviousVol = join(fnDirPrevious, "volume.vol") 
+        fnReferenceVol = join(fnDir, "volumeRef.vol") 
+        if oldXdim != newXDim:
+            self.runJob("xmipp_image_resize", "-i %s -o %s --dim %d" % (fnPreviousVol, fnReferenceVol, newXDim), numberOfMpi=1)
+        else:
+            copyFile(fnPreviousVol, fnReferenceVol)
+        self.nextLowPass = True
+        if self.nextLowPass:
+            self.nextResolutionOffset = 2  # Resolution offset (A)
+            self.runJob('xmipp_transform_filter', '-i %s --fourier low_pass %f --sampling %f' % \
+                        (fnReferenceVol, targetResolution + self.nextResolutionOffset, TsCurrent), numberOfMpi=1)
+        self.nextSpherical = True
+        if self.nextSpherical:
+            R = self.inputParticles.get().getDimensions()[0] / 2
+            self.runJob('xmipp_transform_mask', '-i %s --mask circular -%d' % \
+                        (fnReferenceVol, round(R * self.TsOrig / TsCurrent)), numberOfMpi=1)
+        self.nextPositivity = True
+        if self.nextPositivity:
+            self.runJob('xmipp_transform_threshold', '-i %s --select below 0 --substitute value 0' % fnReferenceVol, numberOfMpi=1)
+
+        if fnMask != '':
+            cleanPath(fnMask)
+        self.writeInfoField(fnDir, "count", emlib.MDL_COUNT, int(2)) 
+
+    def cleanDirectory(self):
+        fnDirCurrent = self._getExtraPath("Iter1") 
+        if self.saveSpace:
+            fnGlobal = fnDirCurrent
+            if exists(fnGlobal):
+                cleanPath(join(fnGlobal, "images.stk"))
+            if exists(fnGlobal):
+                cleanPath(join(fnGlobal, "images.xmd"))
+                cleanPath(join(fnGlobal, "volumeRef.vol"))
+
+    # --------------------------- INFO functions --------------------------------------------
+    def _validate(self):
+        validateMsgs = []
+        # if there are Volume references, it cannot be empty.
+        if self.inputVolume.get() and not self.inputVolume.hasValue():
+            validateMsgs.append('Please provide an input reference volume.')
+        if self.inputParticles.get() and not self.inputParticles.hasValue():
+            validateMsgs.append('Please provide input particles.')
+        return validateMsgs
+    
+    def _summary(self):
+        summary = [
+            "Input particles:  %s" % self.inputParticles.get().getNameId()]
+        summary.append("-----------------")
+        
+        if (not hasattr(self, 'outputParticles')):
+            summary.append("Output not ready yet.")
+        else:
+            if self.percentage == 100.0:
+                text = 'After validation, all of the images are likely to be within the realiable assignment zone'
+            else:
+                text = 'After validation, a %.2f' % self.percentage
+                text += r'% of the images are likely to be within the reliable assignment zone'
+            summary.append(text)
+        return summary    
+
+    def _methods(self):
+        messages = []
+        if (hasattr(self, 'outputParticles')):
+            messages.append('correlation values of previous alignment process'
+                            ' are modified according to the spherical distance from'
+                            ' the assigned direction to a soft and high-valued correlation zone'
+                            ' of neighboring projection directions')
+        return messages
+    
+    def otsu(self, ccList):
+        # from golden-highres
+        import numpy as np
+
+        ccNumber = len(ccList)
+        meanWeight = 1.0 / ccNumber
+        his, bins = np.histogram(ccList, int((max(ccList)-min(ccList))/0.01))
+        finalThres = -1
+        finalVal = -1
+        ccArange = np.arange(min(ccList),max(ccList),0.01)
+        if len(ccArange)==(len(his)+1):
+            ccArange=ccArange[:-1]
+        for t, j in enumerate(bins[1:-1]):
+            idx = t+1
+            pcb = np.sum(his[:idx])
+            pcf = np.sum(his[idx:])
+            Wb = pcb * meanWeight
+            Wf = pcf * meanWeight
+
+            mub = np.sum(ccArange[:idx] * his[:idx]) / float(pcb)
+            muf = np.sum(ccArange[idx:] * his[idx:]) / float(pcf)
+            value = Wb * Wf * (mub - muf) ** 2
+
+            if value > finalVal:
+                finalThres = bins[idx]
+                finalVal = value
+
+        return finalThres
\ No newline at end of file
diff --git a/xmipp3/protocols/protocol_center_particles.py b/xmipp3/protocols/protocol_center_particles.py
index 11982cc2b..bd6a5ef22 100644
--- a/xmipp3/protocols/protocol_center_particles.py
+++ b/xmipp3/protocols/protocol_center_particles.py
@@ -1,6 +1,7 @@
 # ******************************************************************************
 # *
 # * Authors:     Amaya Jimenez Moreno (ajimenez@cnb.csic.es)
+# *              Alberto García Mena (alberto.garcia@cnb.csic.es)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -57,8 +58,8 @@ def _defineParams(self, form):
                       label="Set of micrographs",
                       help='Set of micrographs related to the selected input '
                            'classes')
-        form.addParallelSection(threads=0, mpi=0)
 
+        form.addParallelSection(threads=0, mpi=0)
     # --------------------------- INSERT steps functions -----------------------
     def _insertAllSteps(self):
         self._insertFunctionStep('realignStep')
@@ -136,6 +137,7 @@ def realignStep(self):
                     inPoint = np.array([[0.],[0.],[0.],[1.]])
                     invResultMat = np.linalg.inv(resultMatrix)
                     centerPoint = np.dot(invResultMat,inPoint)
+
                     rowOut.setValue(md.MDL_XCOOR, rowOut.getValue(
                         md.MDL_XCOOR)+int(centerPoint[0]))
                     rowOut.setValue(md.MDL_YCOOR, rowOut.getValue(
@@ -149,6 +151,7 @@ def realignStep(self):
 
     def createOutputStep(self):
         inputParticles = self.inputClasses.get().getImages()
+
         outputClasses = self._createSetOfClasses2D(inputParticles)
         self._fillClasses(outputClasses)
         result = {'outputClasses': outputClasses}
@@ -157,20 +160,12 @@ def createOutputStep(self):
 
         outputParticles = self._createSetOfParticles()
         outputParticles.copyInfo(inputParticles)
-        outputCoords = self._createSetOfCoordinates(self.inputMics.get())
-
-        self._fillParticles(outputParticles, outputCoords, inputParticles)
+        self._fillParticles(outputParticles)
         result = {'outputParticles': outputParticles}
         self._defineOutputs(**result)
         self._defineSourceRelation(self.inputClasses, outputParticles)
 
-        result = {'outputCoordinates': outputCoords}
-        self._defineOutputs(**result)
-        self._defineSourceRelation(self.inputClasses, outputCoords)
-
-
     # --------------------------- UTILS functions ------------------------------
-
     def _fillClasses(self, outputClasses):
         """ Create the SetOfClasses2D """
         inputSet = self.inputClasses.get().getImages()
@@ -189,80 +184,46 @@ def _fillClasses(self, outputClasses):
             newClass.setRepresentative(rep)
             outputClasses.append(newClass)
 
-        i=1
+        i = 1
         mdBlocks = md.getBlocksInMetaDataFile(self._getExtraPath(
             'final_classes.xmd'))
         for block in mdBlocks:
             if block.startswith('class00'):
                 mdClass = md.MetaData(block + "@" + self._getExtraPath(
-                                      'final_classes.xmd'))
+                    'final_classes.xmd'))
                 imgClassId = i
                 newClass = outputClasses[imgClassId]
                 newClass.enableAppend()
                 for row in md.iterRows(mdClass):
                     part = rowToParticle(row)
                     newClass.append(part)
-                i+=1
+                i += 1
                 newClass.setAlignment2D()
                 outputClasses.update(newClass)
 
-    def _fillParticles(self, outputParticles, outputCoords, inputParticles):
+    def _fillParticles(self, outputParticles):
         """ Create the SetOfParticles and SetOfCoordinates"""
         myParticles = md.MetaData(self._getExtraPath('final_images.xmd'))
         outputParticles.enableAppend()
-        outputCoords.enableAppend()
-
-        #Calculating the scale that relates the coordinates with the actual
-        # position in the mic
-        scale = inputParticles.getSamplingRate() / \
-                self.inputMics.get().getSamplingRate()
-        #Dictionary with the name and id of the inpt mics
-        micDictname = {}
-        micDictId = {}
-        for mic in self.inputMics.get():
-            micKey = mic.getMicName()
-            micDictname[micKey] = mic.clone()
-            micKey2 = mic.getObjId()
-            micDictId[micKey2] = mic.clone()
-
         for row in md.iterRows(myParticles):
             #To create the new particle
             p = rowToParticle(row)
-
-            #To create the new coordinate
-            newCoord = Coordinate()
-            coord = p.getCoordinate()
-            if coord.getMicName() is not None:
-                micKey = coord.getMicName()
-                micDict = micDictname
-            else:
-                micKey = coord.getMicId()
-                micDict = micDictId
-            mic = micDict.get(micKey, None)
-            if mic is None:
-                print("Skipping particle, key %s not found" % micKey)
-            else:
-                newCoord.copyObjId(p)
-                x, y = coord.getPosition()
-                newCoord.setPosition(x * scale, y * scale)
-                newCoord.setMicrograph(mic)
-                outputCoords.append(newCoord)
-                p.setCoordinate(newCoord)
-            #Storing the new particle
             outputParticles.append(p)
 
-        boxSize = inputParticles.getXDim() * scale
-        outputCoords.setBoxSize(boxSize)
-
-
 
     # --------------------------- INFO functions -------------------------------
-    def _validate(self):
-        pass
-
     def _summary(self):
         summary = []
         summary.append("Realignment of %s classes."
                        % self.inputClasses.get().getSize())
         return summary
 
+    def _validate(self):
+        errors = []
+        try:
+            self.inputClasses.get().getImages().getAcquisition()
+        except AttributeError:
+            errors.append('InputClasses has not clases')
+
+        return errors
+
diff --git a/xmipp3/protocols/protocol_classification_gpuCorr_semi.py b/xmipp3/protocols/protocol_classification_gpuCorr_semi.py
index 4dc7bf191..3c3c4ba99 100644
--- a/xmipp3/protocols/protocol_classification_gpuCorr_semi.py
+++ b/xmipp3/protocols/protocol_classification_gpuCorr_semi.py
@@ -251,7 +251,8 @@ def _checkNewOutput(self):
                 outputStep.setStatus(STATUS_NEW)
 
     def createOutputStep(self):
-        pass
+        self._closeOutputSet()
+        #pass
 
     # --------------------------- INFO functions -------------------------------
     def _validate(self):
diff --git a/xmipp3/protocols/protocol_compare_angles.py b/xmipp3/protocols/protocol_compare_angles.py
index 30d19825e..3f7997d56 100644
--- a/xmipp3/protocols/protocol_compare_angles.py
+++ b/xmipp3/protocols/protocol_compare_angles.py
@@ -28,6 +28,7 @@
 from pyworkflow import VERSION_2_0
 from pwem.protocols import ProtAnalysis3D
 import pwem.emlib.metadata as md
+from pwem.constants import ALIGN_PROJ
 
 from pwem import emlib
 from xmipp3.convert import (setXmippAttributes, writeSetOfParticles)
@@ -92,8 +93,8 @@ def convertInputStep(self, particlesId1, particlesId2):
         particlesId: is only need to detect changes in
         input particles and cause restart from here.
         """
-        writeSetOfParticles(self.inputParticles1.get(), self._getExtraPath("angles1.xmd"))
-        writeSetOfParticles(self.inputParticles2.get(), self._getExtraPath("angles2.xmd"))
+        writeSetOfParticles(self.inputParticles1.get(), self._getExtraPath("angles1.xmd"), alignType=ALIGN_PROJ)
+        writeSetOfParticles(self.inputParticles2.get(), self._getExtraPath("angles2.xmd"), alignType=ALIGN_PROJ)
 
     def analyzeDistanceStep(self, particlesId1, particlesId2, symmetryGroup):
         self.runJob("xmipp_metadata_utilities","-i %s -o %s --operate keep_column itemId"%\
diff --git a/xmipp3/protocols/protocol_compare_reprojections.py b/xmipp3/protocols/protocol_compare_reprojections.py
index 26ff60a99..cb4e4e611 100644
--- a/xmipp3/protocols/protocol_compare_reprojections.py
+++ b/xmipp3/protocols/protocol_compare_reprojections.py
@@ -1,6 +1,8 @@
 # **************************************************************************
 # *
 # * Authors:     Josue Gomez Blanco (josue.gomez-blanco@mcgill.ca)
+# *              Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *              (produce residuals improved with projection subtraction)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -28,8 +30,7 @@
 import os
 
 from pyworkflow import VERSION_1_1
-from pyworkflow.protocol.params import (PointerParam, StringParam, FloatParam,
-                                        BooleanParam)
+from pyworkflow.protocol.params import (PointerParam, StringParam, FloatParam, BooleanParam, IntParam)
 from pyworkflow.utils.path import cleanPath
 from pwem.protocols import ProtAnalysis3D
 from pwem.objects import (SetOfClasses2D, Image, SetOfAverages, SetOfParticles)
@@ -49,6 +50,7 @@ class XmippProtCompareReprojections(ProtAnalysis3D, ProjMatcher):
     The protocol also analyze the covariance matrix of the residual and computes the logarithm of
     its determinant [Cherian2013]. The extremes of this score (called zScoreResCov), that is
     values particularly low or high, may indicate outliers."""
+
     _label = 'compare reprojections'
     _lastUpdateVersion = VERSION_1_1
     
@@ -56,57 +58,69 @@ def __init__(self, **args):
         ProtAnalysis3D.__init__(self, **args)
         self._classesInfo = dict()
 
-    #--------------------------- DEFINE param functions --------------------------------------------
+    # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputSet', PointerParam, label="Input images", important=True, 
                       pointerClass='SetOfClasses2D, SetOfAverages, SetOfParticles')
         form.addParam('inputVolume', PointerParam, label="Volume to compare images to", important=True,
-                      pointerClass='Volume',
-                      help='Volume to be used for class comparison')
-        form.addParam('useAssignment', BooleanParam, default=True,
-                      label='Use input angular assignment (if available)')
-        form.addParam('optimizeGray', BooleanParam, default=False,
-                      label='Optimize gray scale')
-        form.addParam('ignoreCTF', BooleanParam, default=True,
-                      label='Ignore CTF',
+                      pointerClass='Volume', help='Volume to be used for class comparison')
+        form.addParam('useAssignment', BooleanParam, default=True, label='Use input angular assignment (if available)')
+        form.addParam('optimizeGray', BooleanParam, default=False, label='Optimize gray scale')
+        form.addParam('ignoreCTF', BooleanParam, default=True, 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('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')
-        form.addParam('symmetryGroup', StringParam, default="c1",
-                      label='Symmetry group', 
-                      help='See http://xmipp.cnb.uam.es/twiki/bin/view/Xmipp/Symmetry for a description of the symmetry groups format'
-                        'If no symmetry is present, give c1')
+                      help='If this option is chosen, then the residual covariance matrix is calculated and '
+                           'characterized. But this option takes time and disk space')
+        form.addParam('maskVol', PointerParam, pointerClass='VolumeMask', label='Volume mask', allowsNull=True,
+                      condition='doEvaluateResiduals', expertLevel=LEVEL_ADVANCED,
+                      help='3D mask for the input volume. This mask is not mandatory but advisable.')
+        form.addParam('symmetryGroup', StringParam, default="c1", label='Symmetry group',
+                      help='See http://xmipp.cnb.uam.es/twiki/bin/view/Xmipp/Symmetry for a description of the symmetry'
+                           ' groups format. If no symmetry is present, give c1')
         form.addParam('angularSampling', FloatParam, default=5, expertLevel=LEVEL_ADVANCED,
                       label='Angular sampling rate',
                       help='In degrees.'
                       ' This sampling defines how fine the projection gallery from the volume is explored.')
+        form.addParam('resol', FloatParam, label="Filter at resolution: ", default=3, allowsNull=True,
+                      expertLevel=LEVEL_ADVANCED,
+                      help='Resolution (A) at which subtraction will be performed, filtering the volume projections.'
+                           'Value 0 implies no filtering.')
+        form.addParam('sigma', FloatParam, label="Decay of the filter (sigma): ", default=3, condition='resol',
+                      help='Decay of the filter (sigma parameter) to smooth the mask transition',
+                      expertLevel=LEVEL_ADVANCED)
+        form.addParam('iter', IntParam, label="Number of iterations: ", default=5, expertLevel=LEVEL_ADVANCED,
+                      help='Number of iterations for the adjustment process of the images before the subtraction '
+                           'itself. Several iterations are recommended to improve the adjustment.')
+        form.addParam('rfactor', FloatParam, label="Relaxation factor (lambda): ",  expertLevel=LEVEL_ADVANCED,
+                      default=1, help='Relaxation factor for Fourier amplitude projector (POCS), it should be between '
+                                      '0 (no modification of the amplitudes of the volume 2) and 1 (no relaxation).')
         form.addParallelSection(threads=0, mpi=8)
     
-    #--------------------------- INSERT steps functions --------------------------------------------
+    # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
-        # Convert input images if necessary
-        self.imgsFn = self._getExtraPath('input_imgs.xmd')
+        """ Convert input images if necessary """
+        self.imgsFn = self._getExtraPath('residuals.xmd')
         vol = self.inputVolume.get()
-        
         self._insertFunctionStep("convertStep", self.imgsFn)
-        imgSet = self.inputSet.get()
-        if not self.useAssignment or isinstance(imgSet, SetOfClasses2D) or isinstance(imgSet, SetOfAverages) or (isinstance(imgSet, SetOfParticles) and not imgSet.hasAlignmentProj()):
+
+        if self._useProjMatching():
             anglesFn = self._getExtraPath('angles.xmd')
-            self._insertFunctionStep("projMatchStep", self.inputVolume.get().getFileName(), self.angularSampling.get(), self.symmetryGroup.get(), self.imgsFn,
-                                     anglesFn, self.inputVolume.get().getDim()[0])
+            self._insertFunctionStep("projMatchStep", self.inputVolume.get().getFileName(), self.angularSampling.get(),
+                                     self.symmetryGroup.get(), self.imgsFn, anglesFn, self.inputVolume.get().getDim()[0])
         else:
-            anglesFn=self.imgsFn
+            anglesFn = self.imgsFn
 
         self._insertFunctionStep("produceResiduals", vol.getFileName(), anglesFn, vol.getSamplingRate())
+
         if self.doEvaluateResiduals.get():
             self._insertFunctionStep("evaluateResiduals")
+
         self._insertFunctionStep("createOutputStep")
 
-    #--------------------------- STEPS functions ---------------------------------------------------
+    # --------------------------- STEPS functions ---------------------------------------------------
     def convertStep(self, imgsFn):
         from ..convert import writeSetOfClasses2D, writeSetOfParticles
         imgSet = self.inputSet.get()
@@ -118,38 +132,41 @@ def convertStep(self, imgsFn):
         img = ImageHandler()
         fnVol = self._getTmpPath("volume.vol")
         img.convert(self.inputVolume.get(), fnVol)
-        xdim=self.inputVolume.get().getDim()[0]
-        if xdim!=self._getDimensions():
-            self.runJob("xmipp_image_resize","-i %s --dim %d"%(fnVol,self._getDimensions()),numberOfMpi=1)
-    
-    def produceResiduals(self, fnVol, fnAngles, Ts):
-        fnVol = self._getTmpPath("volume.vol")
-        anglesOutFn=self._getExtraPath("anglesCont.stk")
+        xdim = self.inputVolume.get().getDim()[0]
+        if xdim != self._getDimensions():
+            self.runJob("xmipp_image_resize", "-i %s --dim %d" % (fnVol, self._getDimensions()), numberOfMpi=1)
 
-        projectionsOutFn=self._getExtraPath("projections.stk")
-        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)
+    def produceResiduals(self, fnVol, fnAngles, Ts):
+        anglesOutFn = self._getExtraPath("anglesCont.stk")
+        projectionsOutFn = self._getExtraPath("projections.stk")
+        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)
+        self.fnResiduals = self._getExtraPath("residuals.stk")
         if self.doEvaluateResiduals:
-            args+=" --oresiduals %s"%self._getExtraPath("residuals.stk")
+            args += " --oresiduals %s" % self.fnResiduals
+
         if self.ignoreCTF:
-            args+=" --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):
-            cleanPath(fnNewParticles)
-    
+
+        if self.doEvaluateResiduals:
+            self._computeResiduals(fnVol)
+
     def evaluateResiduals(self):
         # Evaluate each image
         fnAutoCorrelations = self._getExtraPath("autocorrelations.xmd")
         stkAutoCorrelations = self._getExtraPath("autocorrelations.stk")
-        stkResiduals = self._getExtraPath("residuals.stk")
-        anglesOutFn=self._getExtraPath("anglesCont.xmd")
-        self.runJob("xmipp_image_residuals", " -i %s -o %s --save_metadata_stack %s" % (stkResiduals, stkAutoCorrelations, fnAutoCorrelations), numberOfMpi=1)
-        self.runJob("xmipp_metadata_utilities", '-i %s --operate rename_column "image imageResidual"' % fnAutoCorrelations, numberOfMpi=1)
-        self.runJob("xmipp_metadata_utilities", '-i %s --set join %s imageResidual' % (anglesOutFn, fnAutoCorrelations), numberOfMpi=1)
+        stkResiduals = self.fnResiduals
+        anglesOutFn = self._getExtraPath("anglesCont.xmd")
+        self.runJob("xmipp_image_residuals", " -i %s -o %s --save_metadata_stack %s"
+                    % (stkResiduals, stkAutoCorrelations, fnAutoCorrelations), numberOfMpi=1)
+        self.runJob("xmipp_metadata_utilities", '-i %s --operate rename_column "image imageResidual"'
+                    % fnAutoCorrelations, numberOfMpi=1)
+        self.runJob("xmipp_metadata_utilities", '-i %s --set join %s imageResidual'
+                    % (anglesOutFn, fnAutoCorrelations), numberOfMpi=1)
         cleanPath(fnAutoCorrelations)
     
     def createOutputStep(self):
@@ -180,8 +197,7 @@ def createOutputStep(self):
                     outputSet.setAlignmentProj()
             self.iterMd = md.iterRows(imgFn, md.MDL_ITEM_ID)
             self.lastRow = next(self.iterMd)
-            outputSet.copyItems(imgSet,
-                                updateItemCallback=self._processRow)
+            outputSet.copyItems(imgSet, updateItemCallback=self._processRow)
 
         self._defineOutputs(reprojections=outputSet)
         self._defineSourceRelation(self.inputSet, outputSet)
@@ -203,10 +219,10 @@ def _getMdRow(self, mdFile, id):
             if row.getValue(md.MDL_ITEM_ID) == id:
                 return row
 
-        raise Exception("Missing row %s at %s" % (id,mdFile))
+        raise Exception("Missing row %s at %s" % (id, mdFile))
+
     def _processRow(self, particle, row):
         count = 0
-        
         while self.lastRow and particle.getObjId() == self.lastRow.getValue(md.MDL_ITEM_ID):
             count += 1
             if count:
@@ -215,7 +231,6 @@ def _processRow(self, particle, row):
                 self.lastRow = next(self.iterMd)
             except StopIteration:
                 self.lastRow = None
-                    
         particle._appendItem = count > 0
 
     def _createItemMatrix(self, particle, row):
@@ -245,7 +260,66 @@ def __setXmippImage(label):
             __setXmippImage(emlib.MDL_IMAGE_RESIDUAL)
             __setXmippImage(emlib.MDL_IMAGE_COVARIANCE)
 
-    #--------------------------- INFO functions --------------------------------------------
+    def _useProjMatching(self):
+        """ Determine if it is necessary to perform projection matching step (because there is not input alignment)"""
+        imgSet = self.inputSet.get()
+        if not self.useAssignment or isinstance(imgSet, SetOfClasses2D) \
+                or (isinstance(imgSet, SetOfAverages) and not imgSet.hasAlignmentProj()) or \
+                (isinstance(imgSet, SetOfParticles) and not imgSet.hasAlignmentProj()):
+            return True
+        else:
+            return False
+
+    def _computeResiduals(self, fnVol):
+        vol = self.inputVolume.get().clone()
+        involdim = vol.getDim()
+        if fnVol.endswith('.mrc'):
+            fnVol += ':mrc'
+        program = "xmipp_subtract_projection"
+        args = '-i %s --ref %s -o %s --iter %s --lambda %s --subAll' % \
+               (self.imgsFn, fnVol, self._getExtraPath("residuals"), self.iter.get(), self.rfactor.get())
+        args += ' --saveProj %s' % self._getExtraPath('')
+
+        if self.maskVol.get() is not None:
+            mskVol = self.maskVol.get().getFileName()
+        else:
+            fnDescr = self._getExtraPath("mask.descr")
+            fhDescr = open(fnDescr, 'w')
+            fhDescr.write("%d %d %d 0 \n sph + 1 0 0 0 %d" % (involdim[0], involdim[1], involdim[2], involdim[0] / 2))
+            fhDescr.close()
+            mskVol = self._getExtraPath("mask.mrc")
+            args_mask = "-i %s -o %s" % (fnDescr, mskVol)
+            self.runJob("xmipp_phantom_create", args_mask, numberOfMpi=1)
+            args_imageheader = "-i %s --sampling_rate %f" % (mskVol, vol.getSamplingRate())
+            self.runJob("xmipp_image_header", args_imageheader, numberOfMpi=1)
+        args += ' --maskVol %s' % mskVol
+
+        fnDescr2 = self._getExtraPath("mask.descr2")
+        fhDescr2 = open(fnDescr2, 'w')
+        fhDescr2.write("%d %d %d 0" % (involdim[0], involdim[1], involdim[2]))
+        fhDescr2.close()
+        mskKeep = self._getExtraPath("mask_keep.mrc")
+        args_mask_keep = "-i %s -o %s" % (fnDescr2, mskKeep)
+        self.runJob("xmipp_phantom_create", args_mask_keep, numberOfMpi=1)
+        args_imageheader2 = "-i %s --sampling_rate %f" % (mskKeep, vol.getSamplingRate())
+        self.runJob("xmipp_image_header", args_imageheader2, numberOfMpi=1)
+        args += ' --mask %s' % mskKeep
+
+        resol = self.resol.get()
+        if resol:
+            fc = vol.getSamplingRate() / resol
+            args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())
+        self.runJob(program, args, numberOfMpi=1)
+        mrcsresiduals = self._getExtraPath("residuals.mrcs")
+        args2 = " -i %s -o %s" % (mrcsresiduals, self.fnResiduals)
+        self.runJob("xmipp_image_convert", args2, numberOfMpi=1)
+        fnNewParticles = self._getExtraPath("images.stk")
+        if os.path.exists(fnNewParticles):
+            cleanPath(fnNewParticles)
+        if os.path.exists(mrcsresiduals):
+            cleanPath(mrcsresiduals)
+
+    # --------------------------- INFO functions --------------------------------------------
     def _summary(self):
         summary = []
         summary.append("Images evaluated: %i" % self.inputSet.get().getSize())
@@ -255,12 +329,14 @@ def _summary(self):
     def _methods(self):
         methods = []
         if hasattr(self, 'outputParticles'):
-            methods.append("We evaluated %i input images %s regarding to volume %s."\
-                           %(self.inputSet.get().getSize(), self.getObjectTag('inputSet'), self.getObjectTag('inputVolume')) )
-            methods.append("The residuals were evaluated according to their mean, variance and covariance structure [Cherian2013].")
+            methods.append("We evaluated %i input images %s regarding to volume %s." %
+                           (self.inputSet.get().getSize(), self.getObjectTag('inputSet'),
+                            self.getObjectTag('inputVolume')))
+            methods.append("The residuals were evaluated according to their mean, variance and covariance structure "
+                           "[Cherian2013].")
         return methods
     
-    #--------------------------- UTILS functions --------------------------------------------
+    # --------------------------- UTILS functions --------------------------------------------
     def _getDimensions(self):
         imgSet = self.inputSet.get()
         if isinstance(imgSet, SetOfClasses2D):
diff --git a/xmipp3/protocols/protocol_convert_pdb.py b/xmipp3/protocols/protocol_convert_pdb.py
index c58e9800c..f5bb7dde9 100644
--- a/xmipp3/protocols/protocol_convert_pdb.py
+++ b/xmipp3/protocols/protocol_convert_pdb.py
@@ -134,10 +134,10 @@ def convertPdbStep(self):
         if self.setSize:
             args += ' --size'
 
-        if self.size_x.hasValue():
+        if self.size_x.hasValue() and self.setSize:
             args += ' %d' % self.size_x.get()
 
-        if self.size_y.hasValue() and self.size_z.hasValue():
+        if self.size_y.hasValue() and self.size_z.hasValue() and self.setSize:
             args += ' %d %d' % (self.size_y.get(), self.size_z.get())
 
         self.info("Input file: " + pdbFn)
diff --git a/xmipp3/protocols/protocol_core_analysis.py b/xmipp3/protocols/protocol_core_analysis.py
new file mode 100644
index 000000000..caeed6ba6
--- /dev/null
+++ b/xmipp3/protocols/protocol_core_analysis.py
@@ -0,0 +1,226 @@
+# ******************************************************************************
+# *
+# * Authors:     J.M. De la Rosa Trevin (jmdelarosa@cnb.csic.es)
+# *              Carlos Oscar Sánchez Sorzano (coss@cnb.csic.es)
+# *              Daniel Marchán Torres (da.marchan@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 join, exists
+from glob import glob
+
+
+import pyworkflow.protocol.params as param
+import pyworkflow.protocol.constants as const
+from pyworkflow.utils.path import cleanPath, makePath
+
+import pwem.emlib.metadata as md
+from pwem.protocols import ProtClassify2D
+from pwem.constants import ALIGN_2D
+
+from xmipp3.convert import (writeSetOfClasses2D, xmippToLocation,
+                            rowToAlignment)
+
+
+CLASSES_CORE = '_core'
+
+
+class XmippProtCoreAnalysis(ProtClassify2D):
+    """ Analyzes the core of a 2D classification. The core is calculated through
+    the Mahalanobis distance from each image to the center of the class. """
+    
+    _label = 'core analysis'
+    
+    def __init__(self, **args):
+        ProtClassify2D.__init__(self, **args)
+        if self.numberOfMpi.get() < 2:
+            self.numberOfMpi.set(2)
+
+    #--------------------------- DEFINE param functions ------------------------
+    def _defineParams(self, form):
+        form.addSection(label='Input')
+        form.addParam('inputClasses', param.PointerParam,
+                      label="Input classes",
+                      pointerClass='SetOfClasses2D',
+                      help='Set of input classes to be analyzed')
+        form.addParam('thZscore', param.FloatParam, default=3,
+                      label='Junk Zscore',
+                      help='Which is the average Z-score to be considered as '
+                           'junk. Typical values go from 1.5 to 3. For the '
+                           'Gaussian distribution 99.5% of the data is '
+                           'within a Z-score of 3. Lower Z-scores reject more '
+                           'images. Higher Z-scores accept more images.')
+        form.addParam('thPCAZscore', param.FloatParam, default=3,
+                      label='PCA Zscore',
+                      help='Which is the PCA Z-score to be considered as junk. '
+                           'Typical values go from 1.5 to 3. For the Gaussian '
+                           'distribution 99.5% of the data is within a '
+                           'Z-score of 3. Lower Z-scores reject more images. '
+                           'Higher Z-scores accept more images.')
+        form.addParallelSection(threads=0, mpi=4)
+
+    #--------------------------- INSERT steps functions ------------------------
+    def _insertAllSteps(self):
+        self._defineFileNames()
+        self._insertFunctionStep('analyzeCore')
+        self._insertFunctionStep('createOutputStep')
+
+    def analyzeCore(self):
+        # Put in a function convertInputStep
+        fnLevel = self._getExtraPath('level_00')
+        makePath(fnLevel)
+        inputMdName = join(fnLevel, 'level_classes.xmd')
+        writeSetOfClasses2D(self.inputClasses.get(), inputMdName, writeParticles=True)
+
+        args = " --dir %s --root level --computeCore %f %f"%(self._getExtraPath(),
+                                                            self.thZscore, self.thPCAZscore)
+        self.runJob('xmipp_classify_CL2D_core_analysis', args)
+        self.runJob("xmipp_classify_evaluate_classes", "-i %s"%\
+                    self._getExtraPath(join("level_00","level_classes_core.xmd")), numberOfMpi=1)
+
+    #--------------------------- STEPS functions -------------------------------
+    def _defineFileNames(self):
+        """ Centralize how files are called within the protocol. """
+        self.levelPath = self._getExtraPath('level_%(level)02d/')
+        myDict = {
+                  'final_classes': self._getPath('classes2D%(sub)s.sqlite'),
+                  'output_particles': self._getExtraPath('images.xmd'),
+                  'level_classes' : self.levelPath + 'level_classes%(sub)s.xmd',
+                  'level_images' : self.levelPath + 'level_images%(sub)s.xmd',
+                  'classes_scipion': (self.levelPath + 'classes_scipion_level_'
+                                                   '%(level)02d%(sub)s.sqlite'),
+                  }
+        self._updateFilenamesDict(myDict)
+
+    def createOutputStep(self):
+        """ Store the SetOfClasses2D object
+        resulting from the protocol execution.
+        """
+        inputParticles = self.inputClasses.get().getImages()
+        level = self._getLastLevel()
+        subset = CLASSES_CORE
+
+        subsetFn = self._getFileName("level_classes", level=level, sub=subset)
+
+        if exists(subsetFn):
+            classes2DSet = self._createSetOfClasses2D(inputParticles, subset)
+            self._fillClassesFromLevel(classes2DSet, 'last', subset)
+            result = {'outputClasses' + subset: classes2DSet}
+            self._defineOutputs(**result)
+            self._defineSourceRelation(inputParticles, classes2DSet)
+
+    #--------------------------- INFO functions --------------------------------
+    def _validate(self):
+        validateMsgs = []
+        if self.numberOfMpi <= 1:
+            validateMsgs.append('Mpi needs to be greater than 1.')
+        return validateMsgs
+    
+    def _citations(self):
+        citations=['Sorzano2014']
+        return citations
+
+    def _methods(self):
+        strline ='We calculated the class cores %s. [Sorzano2014]' % self.getObjectTag('outputClasses_core')
+        return [strline]
+    
+    #--------------------------- UTILS functions -------------------------------
+    def _updateParticle(self, item, row):
+        item.setClassId(row.getValue(md.MDL_REF))
+        item.setTransform(rowToAlignment(row, ALIGN_2D))
+
+    def _updateClass(self, item):
+        classId = item.getObjId()
+
+        if classId in self._classesInfo:
+            index, fn, _ = self._classesInfo[classId]
+            item.setAlignment2D()
+            rep = item.getRepresentative()
+            rep.setLocation(index, fn)
+            rep.setSamplingRate(self.inputClasses.get().getImages().getSamplingRate())
+
+    def _loadClassesInfo(self, filename):
+        """ Read some information about the produced 2D classes
+        from the metadata file.
+        """
+        self._classesInfo = {}  # store classes info, indexed by class id
+
+        mdClasses = md.MetaData(filename)
+
+        for classNumber, row in enumerate(md.iterRows(mdClasses)):
+            index, fn = xmippToLocation(row.getValue(md.MDL_IMAGE))
+            # Store info indexed by id, we need to store the row.clone() since
+            # the same reference is used for iteration
+            self._classesInfo[index] = (index, fn, row.clone())
+
+    def _fillClassesFromLevel(self, clsSet, level, subset):
+        """ Create the SetOfClasses2D from a given iteration. """
+        classRf = ''
+        self._loadClassesInfo(self._getLevelMdClasses(lev=level, subset=classRf))
+        if subset == '' and level == "last":
+            xmpMd = self._getFileName('output_particles')
+            if not exists(xmpMd):
+                xmpMd = self._getLevelMdImages(level, subset)
+        else:
+            xmpMd = self._getLevelMdImages(level, subset)
+        iterator = md.SetMdIterator(xmpMd, sortByLabel=md.MDL_ITEM_ID,
+                                    updateItemCallback=self._updateParticle,
+                                    skipDisabled=True)
+        # itemDataIterator is not neccesary because, the class SetMdIterator
+        # contain all the information about the metadata
+        clsSet.classifyItems(updateItemCallback=iterator.updateItem,
+                             updateClassCallback=self._updateClass)
+
+    def _getLevelMdClasses(self, lev=0, block="classes", subset=""):
+        """ Return the classes metadata for this iteration.
+        block parameter can be 'info' or 'classes'."""
+        if lev == "last":
+            lev = self._getLastLevel()
+        mdFile = self._getFileName('level_classes', level=lev, sub=subset)
+        if block:
+            mdFile = block + '@' + mdFile
+
+        return mdFile
+
+    def _getLevelMdImages(self, level, subset):
+        if level == "last":
+            level = self._getLastLevel()
+
+        xmpMd = self._getFileName('level_images', level=level, sub=subset)
+        if not exists(xmpMd):
+            self._createLevelMdImages(level, subset)
+
+        return xmpMd
+
+    def _createLevelMdImages(self, level, sub):
+        if level == "last":
+            level = self._getLastLevel()
+        mdClassesFn = self._getLevelMdClasses(lev=level, block="", subset=sub)
+        mdImgs = md.joinBlocks(mdClassesFn, "class0")
+        mdImgs.write(self._getFileName('level_images', level=level, sub=sub))
+
+    def _getLastLevel(self):
+        """ Find the last Level number """
+        clsFn = self._getFileName('level_classes', level=0, sub="")
+        levelTemplate = clsFn.replace('level_00', 'level_??')
+        lev = len(glob(levelTemplate)) - 1
+        return lev
diff --git a/xmipp3/protocols/protocol_ctf_consensus.py b/xmipp3/protocols/protocol_ctf_consensus.py
index fcb29bd7d..e4bae3a74 100644
--- a/xmipp3/protocols/protocol_ctf_consensus.py
+++ b/xmipp3/protocols/protocol_ctf_consensus.py
@@ -48,6 +48,8 @@
 
 ACCEPTED = 'Accepted'
 DISCARDED = 'Discarded'
+INPUT1 = 1
+INPUT2 = 2
 
 class XmippProtCTFConsensus(ProtCTFMicrographs):
     """
@@ -654,7 +656,14 @@ def compareValue(ctf, label, comp, crit):
             minNonAstigmatic, maxNonAstigmatic = \
                 self._getCritNonAstigmaticValidity()
 
-            ctfX = self.xmippCTF[ctfId]
+            if self.xmippCTF == INPUT1:
+                ctfX = ctf
+            elif self.xmippCTF == INPUT2:
+                ctfX = self.allCtf2.get(ctfId)
+            else:
+                raise Exception("No xmipp ctf was found")
+
+
             secondCondition = (
                 compareValue(ctfX, '_xmipp_ctfCritFirstZero', 'lt', firstZero) or
                 compareValue(ctfX, '_xmipp_ctfCritfirstZeroRatio', 'lt', minFirstZero) or
@@ -815,16 +824,16 @@ def _validate(self):
         errors = []
         if self.useCritXmipp.get() and not self.calculateConsensus.get():
             if self.usingXmipp(self.inputCTF.get().getFirstItem()):
-                self.xmippCTF = self.inputCTF.get()
+                self.xmippCTF = INPUT1
             else:
                 errors.append("The primary CTF input ( _Input CTF_ ) must be "
                               "estimated using the _Xmipp - CTF estimation_ "
                               "protocol.")
         if self.useCritXmipp.get() and self.calculateConsensus.get():
             if self.usingXmipp(self.inputCTF.get().getFirstItem()):
-                self.xmippCTF = self.inputCTF.get()
+                self.xmippCTF = INPUT1
             elif self.usingXmipp(self.inputCTF2.get().getFirstItem()):
-                self.xmippCTF = self.inputCTF2.get()
+                self.xmippCTF = INPUT2
             else:
                 errors.append("One of the CTF inputs ( _Input CTF_ or "
                               "_Secundary CTF_) must be estimated using the "
diff --git a/xmipp3/protocols/protocol_denoise_particles.py b/xmipp3/protocols/protocol_denoise_particles.py
index 73e786072..3745bd120 100644
--- a/xmipp3/protocols/protocol_denoise_particles.py
+++ b/xmipp3/protocols/protocol_denoise_particles.py
@@ -28,6 +28,7 @@
 from pyworkflow.object import String
 from pyworkflow.protocol.params import IntParam, PointerParam, LEVEL_ADVANCED
 from pwem.protocols import ProtProcessParticles
+from pwem.objects import SetOfVolumes
 
 from xmipp3.convert import (writeSetOfParticles, writeSetOfClasses2D,
                             xmippToLocation)
@@ -52,7 +53,7 @@ def _defineProcessParams(self, form):
                                                       'respect to the chosen set of classes. This is the standard situation '
                                                       'after CL2D or ML2D.')
         form.addParam('inputClasses', PointerParam, label='Input Classes', important=True,
-                      pointerClass='SetOfClasses', 
+                      pointerClass='SetOfClasses, SetOfAverages',
                       help='Select the input classes for the basis construction against images will be projected to.')
         
         form.addSection(label='Basis construction')
@@ -88,7 +89,10 @@ def denoiseImages(self, inputId, inputClassesId):
         imagesMd = self._getPath('images.xmd')
         writeSetOfParticles(self.inputParticles.get(), imagesMd)
         classesMd = self._getPath('classes.xmd')
-        writeSetOfClasses2D(self.inputClasses.get(), classesMd)
+        if isinstance(self.inputClasses.get(), SetOfVolumes):
+            writeSetOfClasses2D(self.inputClasses.get(), classesMd)
+        else:
+            writeSetOfParticles(self.inputParticles.get(), classesMd)
 
         fnRoot = self._getExtraPath('pca')
         fnRootDenoised = self._getExtraPath('imagesDenoised')
diff --git a/xmipp3/protocols/protocol_extract_particles.py b/xmipp3/protocols/protocol_extract_particles.py
index a11717982..4c9a8853d 100644
--- a/xmipp3/protocols/protocol_extract_particles.py
+++ b/xmipp3/protocols/protocol_extract_particles.py
@@ -68,7 +68,7 @@ def _definePreprocessParams(self, form):
                            'Non-integer downsample factors are possible. ')
 
         form.addParam('boxSize', params.IntParam,
-                      label='Particle box size (px)',
+                      label='Particle box size (px)', allowsPointers=True, 
                       # validators=[params.Positive],
                       help='This is size of the boxed particles (in pixels). '
                            'Note that if you use downsample option, the '
diff --git a/xmipp3/protocols/protocol_particle_pick_consensus.py b/xmipp3/protocols/protocol_particle_pick_consensus.py
index b5e4f967e..342beffd9 100644
--- a/xmipp3/protocols/protocol_particle_pick_consensus.py
+++ b/xmipp3/protocols/protocol_particle_pick_consensus.py
@@ -79,7 +79,7 @@ def _defineParams(self, form):
                       label="Input coordinates", important=True,
                       help='Select the set of coordinates to compare')
         form.addParam('consensusRadius', params.IntParam, default=10,
-                      label="Radius",
+                      label="Radius",  allowsPointers=True, allowsNull=True,
                       help="All coordinates within this radius (in pixels) "
                            "are presumed to correspond to the same particle")
         form.addParam('consensus', params.IntParam, default=-1,
diff --git a/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py b/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
index b2f0202f2..3a1d319bc 100644
--- a/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
+++ b/xmipp3/protocols/protocol_preprocess/protocol_add_noise.py
@@ -48,7 +48,8 @@ class XmippProtAddNoise(ProtRefine3D):
     STUDENT_NOISE = 1
     UNIFORM_NOISE = 2
     _lastUpdateVersion = VERSION_1_1
-    #--------------------------- DEFINE param functions ------------------------
+    # --------------------------- DEFINE param functions ------------------------
+
     def _defineParams(self, form):
         
         form.addParam('noiseType', EnumParam,
@@ -95,13 +96,9 @@ def _defineParams(self, form):
                       condition='noiseType == %d' % self.UNIFORM_NOISE,
                       label="Maximum Value", 
                       help='Please, introduce the maximum value (default = 1).')
-        
-        
         form.addParallelSection(threads=1, mpi=1)
 
-    #--------------------------- INSERT steps functions ------------------------
-
-
+    # --------------------------- INSERT steps functions ------------------------
     def _insertAllSteps(self):        
         self.micsFn = self._getPath()
         # Convert input into xmipp Metadata format
@@ -122,9 +119,8 @@ def _getTypeOfNoise(self):
             noiseParams = '%f %f' % (self.uniformMin, self.uniformMax)
         return kindNoise, noiseParams
 
-    #--------------------------- INFO functions -------------------------------
+    # --------------------------- INFO functions -------------------------------
     def _validate(self):
-        
         validateMsgs = []
         if self.input and not self.input.hasValue():
             validateMsgs.append('Please provide input volume.')  
@@ -132,19 +128,20 @@ def _validate(self):
 
     def _summary(self):
         summary = []
-
-        if  (not hasattr(self,'outputVolume')):
-            summary.append("Output volume not ready yet.")
-        else:
-            summary.append("Volume with %s noise has been obtained" 
-                           % (self.getEnumText("noiseType")) )
-            #
+        if hasattr(self, 'outputVolume'):
+            summary.append("Volume with %s noise has been obtained" % (self.getEnumText("noiseType")))
+        elif hasattr(self, 'outputParticles'):
+            summary.append("Particles with %s noise has been obtained" % (self.getEnumText("noiseType")))
+        elif not hasattr(self, 'outputVolume') or not hasattr(self, 'outputParticles'):
+                summary.append("Output not ready yet.")
         return summary
     
     def _methods(self):
         messages = []
-        if (hasattr(self,'outputVolume')):
+        if hasattr(self, 'outputVolume'):
             messages.append('Noisy volume has been obtained')
+        elif hasattr(self, 'outputParticles'):
+            messages.append('Noisy particles have been obtained')
         return messages
     
     def _citations(self):
@@ -164,7 +161,7 @@ class XmippProtAddNoiseVolumes(XmippProtAddNoise):
     """
     _label = 'add noise volume/s'
     
-    #--------------------------- DEFINE param functions ------------------------
+    # --------------------------- DEFINE param functions ------------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
        
@@ -232,6 +229,7 @@ def _updateNoisyPath(self, vol, row):
     def _isSingleVolume(self):
         return isinstance(self.input.get(), Volume)
 
+
 class XmippProtAddNoiseParticles(XmippProtAddNoise):
     """    
     Given a set of particles, the protocol will add noise to them 
@@ -239,7 +237,7 @@ class XmippProtAddNoiseParticles(XmippProtAddNoise):
     """
     _label = 'add noise particles'
     
-    #--------------------------- DEFINE param functions --------------------
+    # --------------------------- DEFINE param functions --------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
        
diff --git a/xmipp3/protocols/protocol_resolution_deepres.py b/xmipp3/protocols/protocol_resolution_deepres.py
index 60eff8b9a..9f4e0a5dc 100644
--- a/xmipp3/protocols/protocol_resolution_deepres.py
+++ b/xmipp3/protocols/protocol_resolution_deepres.py
@@ -111,7 +111,7 @@ def _createFilenameTemplates(self):
         myDict = {
                  MASK_DILATE: self._getTmpPath('Mask_dilate.vol'),  
                  OPERATE_VOL: self._getTmpPath('operateVolume.vol'),                             
-                 RESIZE_MASK: self._getExtraPath('binaryMask.vol'),                
+                 RESIZE_MASK: self._getTmpPath('binaryMask.vol'),                
                  RESIZE_VOL: self._getExtraPath('originalVolume.vol'),
                  OUTPUT_RESOLUTION_FILE_CHIMERA: self._getExtraPath('chimera_resolution.vol'),
 #                 OUTPUT_RESOLUTION_FILE_CHIMERA: self._getExtraPath(CHIMERA_RESOLUTION_VOL),                                 
@@ -230,6 +230,7 @@ def resolutionStep(self):
         params += ' -s  %f' % self.inputVolume.get().getSamplingRate()            
         params += ' -o  %s' % self._getFileName(OUTPUT_RESOLUTION_FILE)
 
+        os.environ['TF_FORCE_GPU_ALLOW_GROWTH'] = 'true'
         self.runJob("xmipp_deepRes_resolution", params, numberOfMpi=1,
                     env=self.getCondaEnv())
         
diff --git a/xmipp3/protocols/protocol_screen_deepConsensus.py b/xmipp3/protocols/protocol_screen_deepConsensus.py
index c95cafb2e..f7d983453 100644
--- a/xmipp3/protocols/protocol_screen_deepConsensus.py
+++ b/xmipp3/protocols/protocol_screen_deepConsensus.py
@@ -3,6 +3,7 @@
 # * Authors:    Ruben Sanchez Garcia (rsanchez@cnb.csic.es)
 # *             David Maluenda (dmaluenda@cnb.csic.es)
 # *             Daniel Del Hoyo (daniel.delhoyo.gomez@alumnos.upm.es)
+# *             Daniel Marchán (da.marchan@cnb.csic.es)
 # *
 # * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 # *
@@ -44,7 +45,6 @@
 import pyworkflow.utils as pwutils
 from pyworkflow.protocol import params, STATUS_NEW
 import pwem.emlib.metadata as md
-
 from pwem import emlib
 import xmipp3
 from xmipp3 import XmippProtocol
@@ -54,6 +54,9 @@
                             writeSetOfCoordinates, readSetOfCoordsFromPosFnames, readSetOfCoordinates)
 
 MIN_NUM_CONSENSUS_COORDS = 256
+AND = 'by_all'
+OR = 'by_at_least_one'
+UNION_INTERSECTIONS = 'by_at_least_two'
 
 try:
   from xmippPyModules.deepConsensusWorkers.deepConsensus_networkDef import DEEP_PARTICLE_SIZE
@@ -79,6 +82,7 @@ class XmippProtScreenDeepConsensus(ProtParticlePicking, XmippProtocol):
     _stepsCheckSecs = 5              # time in seconds to check the steps
 
     USING_INPUT_COORDS = False
+    USING_INPUT_MICS = False
     CONSENSUS_COOR_PATH_TEMPLATE="consensus_coords_%s"
     CONSENSUS_PARTS_PATH_TEMPLATE="consensus_parts_%s"
     PRE_PROC_MICs_PATH="preProcMics"
@@ -112,6 +116,8 @@ class XmippProtScreenDeepConsensus(ProtParticlePicking, XmippProtocol):
                   'AND': False}
 
     TRAIN_BATCH_MAX = 20
+    PREPROCESS_BATCH_MAX = 200
+    PREDICT_BATCH_MAX = 20
     TO_TRAIN_MICFNS = []
     TRAINED_PARAMS_PATH = 'trainedParams.pickle'
     TRAINING = False
@@ -120,9 +126,11 @@ class XmippProtScreenDeepConsensus(ProtParticlePicking, XmippProtocol):
     LAST_ROUND = False
     ENDED = False
 
+    counter = 0
+
     def __init__(self, **args):
         ProtParticlePicking.__init__(self, **args)
-        self.stepsExecutionMode = params.STEPS_PARALLEL
+        #self.stepsExecutionMode = params.STEPS_PARALLEL
 
     def _defineParams(self, form):
         # GPU settings
@@ -137,7 +145,7 @@ def _defineParams(self, form):
                             " if you do not know what we are talking about."
                             " First core index is 0, second 1 and so on.")
                             
-        form.addParallelSection(threads=1, mpi=1)
+        form.addParallelSection(threads=2, mpi=1)
         
         form.addSection(label='Input')
         
@@ -243,6 +251,10 @@ def _defineParams(self, form):
                            'Warning: Sometimes convergency seems to be reached, '
                            'but after time, improvement can still happen. '
                            'Not recommended for very small data sets (<100 true particles)')
+        form.addParam('maxValAcc', params.FloatParam,
+                      label="Training mean val_acc threshold", default=0.95,
+                      condition="modelInitialization==%s or not skipTraining" % self.ADD_MODEL_TRAIN_NEW,
+                      help='Stop training if at any training batch the selected threshold is achieved')
         form.addParam('l2RegStrength', params.FloatParam,
                       label="Regularization strength",
                       default=1e-5, expertLevel=params.LEVEL_ADVANCED,
@@ -414,7 +426,7 @@ def _defineParams(self, form):
                       label="Training batch size",
                       help='Size of the training batches (in number of micrographs).'
                            'The CNN needs a minimum number of particles to train for each batch, if there are not'
-                           'enough particles, the batch size must be increased')
+                           ' enough particles, the batch size must be increased')
 
 
     def _validate(self):
@@ -476,12 +488,14 @@ def _stepsCheck(self):
         if self.ENDED:
           return
         # Preprocessing
-        if len(self.readyToPreprocessMics(shared = False)) > 0 and not self.PREPROCESSING:
+        if len(self.readyToPreprocessMics(shared=False)) > 0 and not self.PREPROCESSING and (self.counter%2 == 0):
+            #print('----------------------------------ENTERING PREPROCESSING STEP-----------------------------')
             self.PREPROCESSING = True
             self.lastDeps = [self._insertFunctionStep("preprocessMicsStep", prerequisites=self.initDeps)]
 
         #Particle extraction OR (for predictions)
         if len(self.readyToExtractMicFns('OR')) >= self.extractingBatch.get() and not self.EXTRACTING['OR']:
+          #print('----------------------------ENTERING OR EXTRACTION FOR PREDICTION-----------------------------')
           self.EXTRACTING['OR'] = True
           self.newSteps += self.insertCaculateConsensusSteps('OR', prerequisites=self.initDeps)
           self.newSteps += self.insertExtractPartSteps('OR', prerequisites=self.newSteps)
@@ -489,17 +503,20 @@ def _stepsCheck(self):
         #Particle extraction for training and training
         trainedParams = self.loadTrainedParams()
         toTrainSize = self.toTrainDataSize.get() if self.toTrainDataSize.get() != -1 else 1e10
-        if not self.skipTraining.get() and trainedParams['posParticlesTrained'] < toTrainSize:
+        if self.trainingOn() and trainedParams['posParticlesTrained'] < toTrainSize and trainedParams['trainingPass'] != '':
+          #print('----------------------------------ENTERING AND EXTRACTION AND/OR TRAINING--------------------')
           self.doTraining()
 
         elif trainedParams['posParticlesTrained'] >= toTrainSize and trainedParams['trainingPass'] != '':
+          #print('-------------------------ENTERING IN THE CHANGING THE TRAINING_PASS LOGIC-------------------------')
           lastTrainingPass = trainedParams['trainingPass']
           self.retrievePreviousPassModel('', lastTrainingPass)
           trainedParams['trainingPass'] = ''
           self.saveTrainedParams(trainedParams)
 
         #Prediction
-        if self.networkReadyToPredict() and self.cnnFree() and self.predictionsOn():
+        if self.networkReadyToPredict() and self.cnnFree() and self.predictionsOn() and len(self.readyToPredictMicFns()) > 0:
+            #print('---------------------------------------------ENTERING PREDICTION---------------------------------')
             self.PREDICTING = True
             depPredict = self._insertFunctionStep('predictCNN', prerequisites= self.newSteps)
             self.newSteps += [self._insertFunctionStep('endPredictingStep', prerequisites=[depPredict])]
@@ -507,6 +524,7 @@ def _stepsCheck(self):
 
         #Last round with batch size == 1 to include all input
         if self.allFree() and not self.LAST_ROUND and self.checkIfParentsFinished():
+          #print('----------------------------------NOT LAST ROUND BUT ACTIVATES LAST_STEPS---------------------------')
           protLast = self._steps[self.lastStep - 1]
           protLast.addPrerequisites(*self.newSteps)
           protLast.setStatus(STATUS_NEW)
@@ -516,9 +534,11 @@ def _stepsCheck(self):
         protEnd.addPrerequisites(*self.newSteps)
         #Ending the protocol when everything is done
         if self.LAST_ROUND and self.allFree():
+          #print('----------------------------------LAST ROUND ACTIVATES END_STEP---------------------------')
           protEnd.setStatus(STATUS_NEW)
-        self.updateSteps()
 
+        self.updateSteps()
+        self.counter += 1
 
     def endPredictingStep(self):
       self.PREDICTING = False
@@ -526,6 +546,22 @@ def endPredictingStep(self):
     def endTrainingStep(self):
       self.saveTrainedParams(self.curTrainedParams)
       self.TRAINING = False
+      mean_acc = self.loadMeanAccuracy()
+      threshold = self.maxValAcc.get()
+      trainedParams = self.loadTrainedParams()
+      toTrainSize = self.toTrainDataSize.get() if self.toTrainDataSize.get() != -1 else 1e10
+
+      if (mean_acc != None and mean_acc > threshold) or \
+              (trainedParams['posParticlesTrained'] >= toTrainSize and trainedParams['trainingPass'] != ''):
+        #print('-------------------------ENTERING IN THE CHANGING THE TRAININGPASS LOGIC-------------------------')
+        lastTrainingPass = trainedParams['trainingPass']
+        self.retrievePreviousPassModel('', lastTrainingPass)
+        trainedParams['trainingPass'] = ''
+        self.saveTrainedParams(trainedParams)
+        if mean_acc > threshold:
+            print('Mean accuracy %f surpass training accuracy threshold %f -> end training'
+                %(mean_acc, threshold))
+
 
     def lastRoundStep(self):
       '''Starts the last round of training and predictions with the remainign microgrpahs
@@ -545,6 +581,7 @@ def endProtocolStep(self):
         self._insertFunctionStep('createOutputStep', True, prerequisites=[self.depLastPredict])
       else:
         self.updateOutput(closeStream=True)
+        self.ENDED = True
 
 
     def initializeStep(self):
@@ -599,6 +636,7 @@ def initializeStep(self):
 
         #Initializing outputs
         self.inSamplingRate = self._getInputMicrographs().getSamplingRate()
+        self.USING_INPUT_MICS = False
         self.preCorrectedParSet, self.preCoordSet = [], []
 
 
@@ -609,26 +647,29 @@ def preprocessMicsStep(self):
           mics_ = self._getInputMicrographs()
           micsFnameSet = {mics_[micId].getMicName(): mics_[micId].getFileName() for micId in micIds
                           if mics_[micId] is not None}  # to skip failed mics
+          self.USING_INPUT_MICS = False
+          dirName = os.path.dirname(list(micsFnameSet.values())[0])
           toPreprocessMicFns = self.readyToPreprocessMics(shared=False)
-          print('New mics being preprocessed: ', toPreprocessMicFns)
+          print('New mics to be preprocessed: %d' % len(toPreprocessMicFns))
+          if len(toPreprocessMicFns) > self.PREPROCESS_BATCH_MAX:
+              toPreprocessMicFns = toPreprocessMicFns[:self.PREPROCESS_BATCH_MAX]
+
           if self.ignoreCTF.get():
             preproMicsContent="#mics\n"
-            for micName in micsFnameSet:
-              preproMicsContent+= "%s\n"%(micsFnameSet[micName])
+            for micFileName in toPreprocessMicFns:
+              preproMicsContent+= "%s\n"%os.path.join(dirName, micFileName)
           else:
             preproMicsContent="#mics ctfs\n"
             setOfMicCtf= self.ctfRelations.get()
             assert setOfMicCtf is not None, "Error, CTFs must be provided to compute phase flip"
             for ctf in setOfMicCtf:
               ctf_mic = ctf.getMicrograph()
-              ctfMicName = ctf_mic.getMicName()
-              if ctfMicName in micsFnameSet:
+              ctfFnName = ctf_mic.getFileName()
+              if os.path.basename(ctfFnName) in toPreprocessMicFns:
                 ctf_mic.setCTF(ctf)
-                ctfMicName= micsFnameSet[ctfMicName]
-                fnCTF = self._getTmpPath("%s.ctfParam" % os.path.basename(ctfMicName))
+                fnCTF = self._getTmpPath("%s.ctfParam" % os.path.basename(ctfFnName))
                 micrographToCTFParam(ctf_mic, fnCTF)
-                preproMicsContent+= "%s %s\n"%(ctfMicName, fnCTF)
-
+                preproMicsContent+= "%s %s\n"%(ctfFnName, fnCTF)
 
           inputsFname= self._getTmpPath("preprocMic_inputs.txt")
           ouputDir= self._getExtraPath(self.PRE_PROC_MICs_PATH)
@@ -655,122 +696,154 @@ def insertCaculateConsensusSteps(self, mode, prerequisites):
         newDep = self._insertFunctionStep('calculateCoorConsensusStep', outCoordsDataPath, mode, prerequisites=prerequisites)
         newDep = self._insertFunctionStep('loadCoords', outCoordsDataPath, mode, prerequisites=[newDep])
         return [newDep]
-    
+
+    def waitFreeInputMics(self):
+      sl = 0
+      while self.USING_INPUT_MICS:
+        time.sleep(1)
+        sl+=1
+        print('Waiting mics: ', sl)
+        if sl > 60:
+            print('WARNING: a thread has been waiting for more than a minute to access the input micrographs')
+
+      self.USING_INPUT_MICS = True
+
     def waitFreeInputCoords(self):
+      sl = 0
       while self.USING_INPUT_COORDS:
         time.sleep(1)
+        sl += 1
+        print('Waiting coords: ', sl)
+        if sl > 60:
+          print('WARNING: a thread has been waiting for more than a minute to access the input coordinates')
+
       self.USING_INPUT_COORDS = True
 
     def calculateCoorConsensusStep(self, outCoordsDataPath, mode):
-      '''Calculates the consensus coordinates from micrographs whose particles haven't been extracted yet in "mode"'''
-      #Only calculate consensus for coordinates that has not been extracted yet
-      if self.checkIfNewMics(mode):
-        extractedSetOfCoordsFns = []
-        for micFn in self.getExtractedMicFns(mode):
-          extractedSetOfCoordsFns.append(pwutils.path.replaceBaseExt(micFn,'pos'))
+        '''Calculates the consensus coordinates from micrographs whose particles haven't been extracted yet in "mode"'''
+        #Only calculate consensus for coordinates that has not been extracted yet
+        trainedParams = self.loadTrainedParams()
+        if trainedParams['trainingPass'] != '' or mode != 'AND':
+            if self.checkIfNewMics(mode):
+                extractedSetOfCoordsFns = []
+                for micFn in self.getExtractedMicFns(mode):
+                    extractedSetOfCoordsFns.append(pwutils.path.replaceBaseExt(micFn,'pos'))
+
+                self.waitFreeInputCoords()
+                inputCoordsFnames = self.getInpCoordsFns(mode, extractedSetOfCoordsFns)
+                self.USING_INPUT_COORDS = False
+
+                inputFileHeader="#pos_i\n"
+                inputFileStr=inputFileHeader
+                for baseName in inputCoordsFnames:
+                    fnames= inputCoordsFnames[baseName]
+                    inputFileStr+=" ".join(fnames)+"\n"
+
+                assert len(inputFileStr)>len(inputFileHeader), "Error, no consensus can be computed as there " \
+                                                               "is a mismatch in coordinate sets filenames"
+                #Concensus decision
+                consensus = UNION_INTERSECTIONS if mode=="AND" else OR
+                # THE MODE CAN BE A NEW PARAMETER TO BE ASKED TO THE USER
+                configFname= self._getTmpPath("consensus_%s_inputs.txt"%(mode) )
+                with open(configFname, "w") as f:
+                    f.write(inputFileStr)
+
+                args = "-i %s -s %d -c %s -d %f -o %s -t %d" % (configFname, self._getBoxSize(), consensus, self.consensusRadius.get(),
+                                                                outCoordsDataPath, self.numberOfThreads.get())
+                self.runJob('xmipp_coordinates_consensus', args, numberOfMpi=1)
+                self.TO_EXTRACT_MICFNS[mode] = self.readyToExtractMicFns(mode)
+
+                return
+
 
-        self.waitFreeInputCoords()
-        inputCoordsFnames = self.getInpCoordsFns(mode, extractedSetOfCoordsFns)
-        self.USING_INPUT_COORDS = False
-          
-        inputFileHeader="#pos_i\n"
-        inputFileStr=inputFileHeader
-        for baseName in inputCoordsFnames:
-          fnames= inputCoordsFnames[baseName]
-          inputFileStr+=" ".join(fnames)+"\n"
-
-        assert len(inputFileStr)>len(inputFileHeader), "Error, no consensus can be computed as there " \
-                                                       "are mismatch in coordinate sets filenames"
-        consensus = -1 if mode=="AND" else 1
-        configFname= self._getTmpPath("consensus_%s_inputs.txt"%(mode) )
-        with open(configFname, "w") as f:
-            f.write(inputFileStr)
-
-        args="-i %s -s %d -c %d -d %f -o %s -t %d"%(configFname, self._getBoxSize(), consensus, self.consensusRadius.get(),
-                                                   outCoordsDataPath, self.numberOfThreads.get())
-        self.runJob('xmipp_coordinates_consensus', args, numberOfMpi=1)
-        self.TO_EXTRACT_MICFNS[mode] = self.readyToExtractMicFns(mode)
 
     def pickNoise(self):
         '''Find noise coordinates from micrographs in order to use them as negatives in the training process'''
-        orPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "OR")
-        outputPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "NOISE")
-        if not "OR" in self.coordinatesDict:  # fill self.coordinatesDict['OR']
-          self.loadCoords(orPosDir, 'OR')
-        # Getting the extracted mics where noise have not been picked yet
-        micsDir = self._getExtraPath(self.PRE_PROC_MICs_PATH)
-        toPickNoiseFns = list(set(self.getExtractedMicFns('OR')) -
-                              (set(self.getExtractedMicFns('NOISE')) | set(self.TO_EXTRACT_MICFNS['NOISE'])))
-        if len(toPickNoiseFns) > 0:
-          micNoiseSet = self.loadMicSetFromFns(micsDir, toPickNoiseFns)
-          toNoiseSetOfCoords = self._createSetOfCoordinates(micNoiseSet)
-          toNoiseSetOfCoords.setBoxSize(self._getBoxSize())
-          readSetOfCoordinates(orPosDir, micSet=micNoiseSet, coordSet=toNoiseSetOfCoords)
-
-          # Write the tonoise mic files into a tmp directory
-          coordsDir = self._getTmpPath(IN_COORDS_POS_DIR_BASENAME)
-          toPickMicsDir = micsDir + '_toPickNoise'
-          toPickCoorsDir = coordsDir + '_toPickNoise'
-          if not os.path.exists(coordsDir):
-            makePath(coordsDir)
-          if os.path.exists(toPickMicsDir):
-            cleanPath(toPickMicsDir)
-          if os.path.exists(toPickCoorsDir):
-            cleanPath(toPickCoorsDir)
-          makePath(toPickMicsDir)
-          makePath(toPickCoorsDir)
-          writeSetOfCoordinates(coordsDir, toNoiseSetOfCoords)
-
-          for micFn in toPickNoiseFns:
-            posFn = pwutils.path.replaceBaseExt(micFn, 'pos')
-            shutil.copyfile(micsDir + '/' + micFn, toPickMicsDir + '/' + micFn)
-            shutil.copyfile(coordsDir + '/' + posFn, toPickCoorsDir + '/' + posFn)
-
-          argsDict = pickNoise_prepareInput(toNoiseSetOfCoords, self._getTmpPath())
-          argsDict['toPickMicsDir'] = toPickMicsDir
-          argsDict['toPickCoorsDir'] = toPickCoorsDir
-          argsDict["outputPosDir"] = outputPosDir
-          argsDict["nThrs"] = self.numberOfThreads.get()
-          argsDict["nToPick"] = -1
-          args = (" -i %(toPickMicsDir)s -c %(toPickCoorsDir)s -o %(outputPosDir)s -s %(boxSize)s " +
-                  "-n %(nToPick)s -t %(nThrs)s") % argsDict
-
-          self.runJob('xmipp_pick_noise', args, numberOfMpi=1)
-
-          self.loadCoords(outputPosDir, 'NOISE', micSet=micNoiseSet)
-          self.TO_EXTRACT_MICFNS['NOISE'] = toPickNoiseFns
-          print('Adding to extract {} {} micrographs'.format(len(toPickNoiseFns), 'NOISE'))
+        trainedParams = self.loadTrainedParams()
+        if trainedParams['trainingPass'] != '':
+            orPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "OR")
+            outputPosDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % "NOISE")
+            if not "OR" in self.coordinatesDict:  # fill self.coordinatesDict['OR']
+              self.loadCoords(orPosDir, 'OR')
+            # Getting the extracted mics where noise have not been picked yet
+            micsDir = self._getExtraPath(self.PRE_PROC_MICs_PATH)
+            toPickNoiseFns = list(set(self.getExtractedMicFns('OR')) -
+                                  (set(self.getExtractedMicFns('NOISE')) | set(self.TO_EXTRACT_MICFNS['NOISE'])))
+            if len(toPickNoiseFns) > 0:
+              micNoiseSet = self.loadMicSetFromFns(micsDir, toPickNoiseFns)
+              toNoiseSetOfCoords = self._createSetOfCoordinates(micNoiseSet)
+              toNoiseSetOfCoords.setBoxSize(self._getBoxSize())
+              readSetOfCoordinates(orPosDir, micSet=micNoiseSet, coordSet=toNoiseSetOfCoords)
+
+              # Write the tonoise mic files into a tmp directory
+              coordsDir = self._getTmpPath(IN_COORDS_POS_DIR_BASENAME)
+              toPickMicsDir = micsDir + '_toPickNoise'
+              toPickCoorsDir = coordsDir + '_toPickNoise'
+              if not os.path.exists(coordsDir):
+                makePath(coordsDir)
+              if os.path.exists(toPickMicsDir):
+                cleanPath(toPickMicsDir)
+              if os.path.exists(toPickCoorsDir):
+                cleanPath(toPickCoorsDir)
+              makePath(toPickMicsDir)
+              makePath(toPickCoorsDir)
+              writeSetOfCoordinates(coordsDir, toNoiseSetOfCoords)
+
+              for micFn in toPickNoiseFns:
+                posFn = pwutils.path.replaceBaseExt(micFn, 'pos')
+                shutil.copyfile(micsDir + '/' + micFn, toPickMicsDir + '/' + micFn)
+                shutil.copyfile(coordsDir + '/' + posFn, toPickCoorsDir + '/' + posFn)
+
+              argsDict = pickNoise_prepareInput(toNoiseSetOfCoords, self._getTmpPath())
+              argsDict['toPickMicsDir'] = toPickMicsDir
+              argsDict['toPickCoorsDir'] = toPickCoorsDir
+              argsDict["outputPosDir"] = outputPosDir
+              argsDict["nThrs"] = self.numberOfThreads.get()
+              argsDict["nToPick"] = -1
+              args = (" -i %(toPickMicsDir)s -c %(toPickCoorsDir)s -o %(outputPosDir)s -s %(boxSize)s " +
+                      "-n %(nToPick)s -t %(nThrs)s") % argsDict
+
+              self.runJob('xmipp_pick_noise', args, numberOfMpi=1)
+
+              self.loadCoords(outputPosDir, 'NOISE', micSet=micNoiseSet)
+              self.TO_EXTRACT_MICFNS['NOISE'] = toPickNoiseFns
+              print('Adding to extract {} {} micrographs'.format(len(toPickNoiseFns), 'NOISE'))
+
+              return
 
     def loadCoords(self, posCoorsPath, mode, micSet=[]):
         #Upload coords sqlite
-        if not micSet == []:
-          #Load coordinates from an specific set of mics
-          batchSetOfCoordinates = self._createSetOfCoordinates(micSet)
-          batchSetOfCoordinates.setBoxSize(self._getBoxSize())
-          readSetOfCoordinates(posCoorsPath, micSet=micSet, coordSet = batchSetOfCoordinates)
-          if mode in self.coordinatesDict:
-            for newCoord in batchSetOfCoordinates:
-              apCoord = Coordinate()
-              apCoord.copy(newCoord, copyId=False)
-              self.coordinatesDict[mode].append(apCoord)
-          else:
-            self.coordinatesDict[mode] = batchSetOfCoordinates
-        else:
-          sqliteName = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode) + ".sqlite"
-          if os.path.isfile(self._getExtraPath(sqliteName)):
-            cleanPath(self._getExtraPath(sqliteName))
-          
-          self.waitFreeInputCoords()
-          totalSetOfCoordinates = readSetOfCoordsFromPosFnames(posCoorsPath,
-                                                               setOfInputCoords=self.inputCoordinates[0].get(),
-                                                               sqliteOutName=sqliteName, write=True)
-          self.USING_INPUT_COORDS=False
-          print("Coordinates %s size: %d" % (mode, totalSetOfCoordinates.getSize()))
-          assert totalSetOfCoordinates.getSize() > MIN_NUM_CONSENSUS_COORDS, \
-            ("Error, the consensus (%s) of your input coordinates was too small (%s). " +
-             "It must be > %s. Try a different input..."
-             ) % (mode, str(totalSetOfCoordinates.getSize()), str(MIN_NUM_CONSENSUS_COORDS))
-          self.coordinatesDict[mode] = totalSetOfCoordinates
+        trainedParams = self.loadTrainedParams()
+        if trainedParams['trainingPass'] != '' or mode != 'AND':
+            if len(micSet):
+                  #Load coordinates from an specific set of mics
+                  batchSetOfCoordinates = self._createSetOfCoordinates(micSet)
+                  batchSetOfCoordinates.setBoxSize(self._getBoxSize())
+                  readSetOfCoordinates(posCoorsPath, micSet=micSet, coordSet = batchSetOfCoordinates)
+                  if mode in self.coordinatesDict:
+                    for newCoord in batchSetOfCoordinates:
+                      apCoord = Coordinate()
+                      apCoord.copy(newCoord, copyId=False)
+                      self.coordinatesDict[mode].append(apCoord)
+                  else:
+                       self.coordinatesDict[mode] = batchSetOfCoordinates
+            else:
+                  sqliteName = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode) + ".sqlite"
+                  if os.path.isfile(self._getExtraPath(sqliteName)):
+                    cleanPath(self._getExtraPath(sqliteName))
+
+                  self.waitFreeInputCoords()
+                  totalSetOfCoordinates = readSetOfCoordsFromPosFnames(posCoorsPath,
+                                                                       setOfInputCoords=self.inputCoordinates[0].get(),
+                                                                       sqliteOutName=sqliteName, write=True)
+                  print("Coordinates %s size: %d" % (mode, totalSetOfCoordinates.getSize()))
+                  assert totalSetOfCoordinates.getSize() > MIN_NUM_CONSENSUS_COORDS, \
+                    ("Error, the consensus (%s) of your input coordinates was too small (%s). " +
+                     "It must be > %s. Try a different input..."
+                     ) % (mode, str(totalSetOfCoordinates.getSize()), str(MIN_NUM_CONSENSUS_COORDS))
+                  self.coordinatesDict[mode] = totalSetOfCoordinates
+
+            self.USING_INPUT_COORDS = False
 
 
     def insertExtractPartSteps(self, mode, prerequisites):
@@ -782,78 +855,89 @@ def insertExtractPartSteps(self, mode, prerequisites):
 
     def extractParticles(self, mode):
         '''Extract the particles from a set of micrographs with their corresponding coordinates'''
-        micsFnameSet = {}
-        posDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode)
-        preprocMicsPath = self._getExtraPath(self.PRE_PROC_MICs_PATH)
-        toExtractMicFns = self.TO_EXTRACT_MICFNS[mode]
-        print('To extract in mode {}: {}'.format(mode, toExtractMicFns))
-        if len(toExtractMicFns) <= 0:
-          return
+        trainedParams = self.loadTrainedParams()
+        if trainedParams['trainingPass'] == '' and (mode.startswith("ADDITIONAL_COORDS") or mode == 'AND' or mode == 'NOISE'):
+            print('Training already finished')
 
-        for micFname in toExtractMicFns:
-          micFnameBase= pwutils.removeExt(micFname)
-          micFname= os.path.join(preprocMicsPath, micFname)
-          micsFnameSet[micFnameBase]= micFname
-        extractCoordsContent="#mics coords\n"
-
-        if mode.startswith("ADDITIONAL_COORDS"):
-          if not os.path.exists(posDir):
-            os.mkdir(posDir)
-          if mode.endswith("TRUE"):
-            coordSet= self.trainTrueSetOfCoords.get()
-          elif mode.endswith("FALSE"):
-            coordSet= self.trainFalseSetOfCoords.get()
-          writeSetOfCoordinates(posDir, coordSet)
-
-        for posFname in os.listdir(posDir):
-          posNameBase=  pwutils.removeExt(posFname)
-          posFname= os.path.join(posDir, posFname)
-          if posNameBase in micsFnameSet:
-            extractCoordsContent+= "%s particles@%s\n"%(micsFnameSet[posNameBase], posFname)
-
-        inputsFname= self._getTmpPath("extractParticles_inputs_%s.txt"%mode)
-        outputDir= self._getConsensusParticlesDir(mode)
-        if not os.path.exists(outputDir):
-          makePath(outputDir)
-        nThrs= self.numberOfThreads.get()
-        with open(inputsFname, "w") as f:
-          f.write(extractCoordsContent)
-        downFactor= self._getDownFactor()
-        args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, DEEP_PARTICLE_SIZE, downFactor, outputDir, nThrs)
-
-        self.runJob('xmipp_extract_particles', args, numberOfMpi=1)
+        else:
+            micsFnameSet = {}
+            posDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode)
+            preprocMicsPath = self._getExtraPath(self.PRE_PROC_MICs_PATH)
+            toExtractMicFns = self.TO_EXTRACT_MICFNS[mode]
+            print('To extract {} in mode {}: {}'.format(len(toExtractMicFns), mode, toExtractMicFns))
+            if len(toExtractMicFns) <= 0:
+              return
+
+            for micFname in toExtractMicFns:
+              micFnameBase= pwutils.removeExt(micFname)
+              micFname= os.path.join(preprocMicsPath, micFname)
+              micsFnameSet[micFnameBase]= micFname
+            extractCoordsContent="#mics coords\n"
+
+            if mode.startswith("ADDITIONAL_COORDS"):
+              if not os.path.exists(posDir):
+                os.mkdir(posDir)
+              if mode.endswith("TRUE"):
+                coordSet= self.trainTrueSetOfCoords.get()
+              elif mode.endswith("FALSE"):
+                coordSet= self.trainFalseSetOfCoords.get()
+              writeSetOfCoordinates(posDir, coordSet)
+
+            for posFname in os.listdir(posDir):
+              posNameBase=  pwutils.removeExt(posFname)
+              posFname= os.path.join(posDir, posFname)
+              if posNameBase in micsFnameSet:
+                extractCoordsContent+= "%s particles@%s\n"%(micsFnameSet[posNameBase], posFname)
+
+            inputsFname= self._getTmpPath("extractParticles_inputs_%s.txt"%mode)
+            outputDir= self._getConsensusParticlesDir(mode)
+            if not os.path.exists(outputDir):
+              makePath(outputDir)
+            nThrs= self.numberOfThreads.get()
+            with open(inputsFname, "w") as f:
+              f.write(extractCoordsContent)
+            downFactor= self._getDownFactor()
+            args= "-i %s -s %s -d %s -o %s -t %d"%(inputsFname, DEEP_PARTICLE_SIZE, downFactor, outputDir, nThrs)
+
+            self.runJob('xmipp_extract_particles', args, numberOfMpi=1)
 
     def joinSetOfParticlesStep(self, mode, micFns='', trainingPass='', clean=False):
         '''Stores the particles extracted from a set of micrographs in a images.xmd metadata file'''
-        fnImages = self._getExtraPath("particles_{}{}.xmd".format(mode, trainingPass))
-        posDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode)
-        if micFns == '':
-          micFns = self.TO_EXTRACT_MICFNS[mode]
-          self.TO_EXTRACT_MICFNS[mode] = []
-          print('Updating set of particles written in {}'.format(fnImages))
-        else:
-          print('Creating set of particles in {}'.format(fnImages))
-
-        imgsXmd = md.MetaData()
-        for micFn in micFns:
-          posFn = os.path.join(posDir, pwutils.replaceBaseExt(micFn, "pos"))
-          xmdFn = os.path.join(self._getConsensusParticlesDir(mode),
-                               pwutils.replaceBaseExt(posFn, "xmd"))
-          if os.path.exists(xmdFn):
-            mdFn = md.MetaData(xmdFn)
-            mdPos = md.MetaData('particles@%s' % posFn)
-            mdPos.merge(mdFn)
-            imgsXmd.unionAll(mdPos)
-          else:
-            self.warning("The coord file %s wasn't used for extraction! "
-                         % os.path.basename(posFn))
-        if not os.path.exists(fnImages) or clean:
-          imgsXmd.write(fnImages)
+        trainedParams = self.loadTrainedParams()
+        if trainedParams['trainingPass'] == '' and (mode.startswith("ADDITIONAL_COORDS") or mode == 'AND' or mode == 'NOISE'):
+            print('Training already finished')
+
         else:
-          imgsXmd.append(fnImages)
+            fnImages = self._getExtraPath("particles_{}{}.xmd".format(mode, trainingPass))
+            posDir = self._getExtraPath(self.CONSENSUS_COOR_PATH_TEMPLATE % mode)
+            if micFns == '':
+              micFns = self.TO_EXTRACT_MICFNS[mode]
+              self.TO_EXTRACT_MICFNS[mode] = []
+              print('Updating set of particles written in {}'.format(fnImages))
+            else:
+              print('Creating set of particles in {}'.format(fnImages))
+
+            imgsXmd = md.MetaData()
+            for micFn in micFns:
+              posFn = os.path.join(posDir, pwutils.replaceBaseExt(micFn, "pos"))
+              xmdFn = os.path.join(self._getConsensusParticlesDir(mode),
+                                   pwutils.replaceBaseExt(posFn, "xmd"))
+              if os.path.exists(xmdFn):
+                mdFn = md.MetaData(xmdFn)
+                mdPos = md.MetaData('particles@%s' % posFn)
+                mdPos.merge(mdFn)
+                imgsXmd.unionAll(mdPos)
+              else:
+                self.warning("The coord file %s wasn't used for extraction! "
+                             % os.path.basename(posFn))
+            if not os.path.exists(fnImages) or clean:
+              imgsXmd.write(fnImages)
+            else:
+              imgsXmd.append(fnImages)
+
         self.EXTRACTING[mode] = False
 
-    def doTraining(self):
+    def doTraining(self): # PUT SOMETHING LIKE IN THE PCIKNOISE BUT MAKE SURE YOU ARE ONLY USING IT FOR AND SETS
       '''Prepares the positive (AND) and negative (NOISE) coordinates for the training and executes it'''
       if len(self.readyToExtractMicFns('NOISE')) >= self.extractingBatch.get() and not self.EXTRACTING['NOISE']:
         self.EXTRACTING['NOISE'] = True
@@ -981,6 +1065,10 @@ def predictCNN(self):
           mdORPath = self._getExtraPath("particles_OR.xmd")
         else:
           toPredictMicFns = self.readyToPredictMicFns()
+          print("Mics ready to predict {}".format(len(toPredictMicFns)))
+          if len(toPredictMicFns) > self.PREDICT_BATCH_MAX:
+              toPredictMicFns = toPredictMicFns[:self.PREDICT_BATCH_MAX]
+
           predExten = '_partial'
           self.joinSetOfParticlesStep(mode='OR', micFns=toPredictMicFns, trainingPass=predExten, clean=True)
           mdORPath = self._getExtraPath("particles_OR{}.xmd".format(predExten))
@@ -1001,6 +1089,7 @@ def predictCNN(self):
 
         mdObject = md.MetaData(mdORPath)
         print('Predicting on {} true particles'.format(mdObject.size()))
+              #'in {} micrographs'.format(mdObject.size(), len(toPredictMicFns)))
         predictDict = {mdORPath: 1}
 
         if self.doTesting.get() and self.testTrueSetOfParticles.get() and self.testFalseSetOfParticles.get() and not\
@@ -1053,7 +1142,6 @@ def createFinalOutput(self, closeStream=False):
       partSet.setSamplingRate(self._getDownFactor() * self.inSamplingRate)
 
       self.outputParticles, self.outputCoordinates = self.getParticlesOutput(partSet), self.getCoordinatesOutput()
-
       downFactor = self._getDownFactor()
       for part in partSet:
         coord = part.getCoordinate().clone()
@@ -1064,7 +1152,7 @@ def createFinalOutput(self, closeStream=False):
         newPart.copy(part, copyId=False)
         newPart.scaleCoordinate(downFactor)
         if (self.threshold.get() < 0 or
-                getattr(newPart, deepZscoreLabel) > self.threshold.get()):
+          getattr(newPart, deepZscoreLabel) > self.threshold.get()):
           self.outputCoordinates.append(coord)
           self.outputParticles.append(newPart)
 
@@ -1246,6 +1334,7 @@ def getAllCoordsInputMicrographs(self, shared=False):
       all microgrpah present in any set (Intersection vs Union)
       Do not create a set, because of concurrency in the database'''
       self.waitFreeInputCoords()
+      self.waitFreeInputMics()
       micDict, micFns = {}, set([])
       for inputCoord in self.inputCoordinates:
         newMics = inputCoord.get().getMicrographs()
@@ -1259,7 +1348,7 @@ def getAllCoordsInputMicrographs(self, shared=False):
           micFns = micFns | set(newMicFns)
         else:
           micFns = micFns & set(newMicFns)
-      self.USING_INPUT_COORDS = False
+      self.USING_INPUT_COORDS, self.USING_INPUT_MICS = False, False
       sharedMicDict = {}
       for micFn in micFns:
         sharedMicDict[micFn] = micDict[micFn]
@@ -1268,6 +1357,7 @@ def getAllCoordsInputMicrographs(self, shared=False):
 
     def _getInputMicrographs(self):
       '''Return a list with the micrographs corresponding the input coordinates'''
+      self.waitFreeInputMics()
       if not hasattr(self, "inputMicrographs") or not self.inputMicrographs:
         self.waitFreeInputCoords()
         if len(self.inputCoordinates) == 0:
@@ -1292,8 +1382,7 @@ def _getBoxSize(self):
 
     def _getDownFactor(self):
       if not hasattr(self, "downFactor") or not self.downFactor:
-        firstCoords = self._getInputMicrographs()
-        self.boxSize = firstCoords.getBoxSize()
+        self.boxSize = self._getBoxSize()
         self.downFactor = self.boxSize / float(DEEP_PARTICLE_SIZE)
         assert self.downFactor >= 1, \
           "Error, the particle box size must be greater or equal than 128."
@@ -1311,12 +1400,13 @@ def loadMicSetFromFns(self, inputDir, micFns):
         micPath = os.path.join(inputDir, micFn)
         micSet.append(Micrograph(micPath))
       micSet.copyInfo(self._getInputMicrographs())
+      self.USING_INPUT_MICS = False
       return micSet
 
     def getPreProcParamsFromForm(self):
         mics_= self._getInputMicrographs()
-        mic = mics_.getFirstItem()
-        fnMic = mic.getFileName()
+        fnMic = mics_.getFirstItem().getFileName()
+        self.USING_INPUT_MICS = False
         pathToMics= os.path.split(fnMic)[0]
         pathToCtfs= "None"
 
@@ -1353,7 +1443,26 @@ def _getEffectiveNumPartsTrain(self, dictTrueData):
           nParts+= mdObject.size()
         return nParts
 
+    def loadMeanAccuracy(self):
+        trainedParams = self.curTrainedParams
+        trPass = trainedParams['trainingPass']
+        netDataPath = self._getExtraPath(self.NET_TEMPLATE.format(trPass))
+        netMeanAccFname = os.path.join(netDataPath, "netsMeanValAcc.txt")
+        if os.path.exists(netDataPath):
+            with open(netMeanAccFname) as f:
+                lines = f.readlines()
+                mean_accuracy = float(lines[0].split()[1])
+            return mean_accuracy
+        else:
+            return None
+
     #STREAMING and state checks
+    def trainingOn(self):
+      '''Return a boolean for whether to perform training. True if the training must not be skipped and if the finish
+      training criteria has not been reached yet'''
+      trainedParams = self.loadTrainedParams()
+      return not self.skipTraining.get() and trainedParams['keepTraining']
+
     def predictionsOn(self):
       '''Return a boolean for whether to perform a prediction. True if there must be a preliminar prediction or
       if the training process has finished (trainingPass=='') '''
@@ -1386,6 +1495,7 @@ def allFree(self):
       for mode in ['OR', 'NOISE', 'AND']:
         if self.EXTRACTING[mode]:
           gExtracting = True
+
       return not self.PREDICTING and not self.TRAINING and not gExtracting and not self.PREPROCESSING
 
     def checkIfParentsFinished(self):
@@ -1416,7 +1526,9 @@ def checkIfNewMics(self, mode=''):
     def getMicsIds(self, filterOutNoCoords=False):
         '''Returns the input micrographs Ids'''
         if not filterOutNoCoords:
-          return self._getInputMicrographs().getIdSet()
+          idSet = self._getInputMicrographs().getIdSet()
+          self.USING_INPUT_MICS = False
+          return idSet
         self.waitFreeInputCoords()
         micFnames, micIds = set([]), set([])
         for coordinatesP in self.inputCoordinates:
@@ -1527,7 +1639,8 @@ def loadTrainedParams(self):
                   'trainingPass': 0,
                   'predictionPasses': [],
                   'doneExtraTesting': False,
-                  'firstTraining': True
+                  'firstTraining': True,
+                  'keepTraining': True,
                   }
       return params
 
diff --git a/xmipp3/protocols/protocol_shift_particles.py b/xmipp3/protocols/protocol_shift_particles.py
index 162834ddb..680c072d0 100644
--- a/xmipp3/protocols/protocol_shift_particles.py
+++ b/xmipp3/protocols/protocol_shift_particles.py
@@ -61,6 +61,12 @@ def _defineParams(self, form):
         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('applyShift', BooleanParam, label='Apply shift to particles?', default='True',
+                      help='Yes: The shift is applied to particle images and zero shift is stored in the metadata. No: '
+                           'The shift is stored in the transformation matrix in the metadata, but not applied to the '
+                           'particle image (i.e. the output images are the same of input images). This option takes '
+                           'less time and the shift could be applied later using protocol "xmipp3 - apply alignment 2d"'
+                           ' or by re-extracting the particles.')
         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',
@@ -109,7 +115,9 @@ def shiftStep(self):
             interp = 'linear'
         else:
             interp = 'spline'
-        args += ' --apply_transform --dont_wrap --interp %s' % interp
+        if self.applyShift.get():
+            args += ' --apply_transform'
+        args += ' --dont_wrap --interp %s' % interp
         if self.inv.get():
             args += ' --inverse'
         self.runJob(program, args)
diff --git a/xmipp3/protocols/protocol_subtract_projection.py b/xmipp3/protocols/protocol_subtract_projection.py
index 5843416a5..ce183d8af 100644
--- a/xmipp3/protocols/protocol_subtract_projection.py
+++ b/xmipp3/protocols/protocol_subtract_projection.py
@@ -51,7 +51,8 @@ def _defineParams(self, form):
         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.')
+                      help='Specify a 3D mask for the region of the input volume that you want to keep. If no mask is '
+                           'selected, the whole image will be subtracted')
         form.addParam('resol', FloatParam, label="Filter at resolution: ", default=3, allowsNull=True,
                       expertLevel=LEVEL_ADVANCED,
                       help='Resolution (A) at which subtraction will be performed, filtering the volume projections.'
@@ -76,7 +77,7 @@ def _insertAllSteps(self):
     # --------------------------- 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()
@@ -96,6 +97,18 @@ def subtractionStep(self):
             args += ' --maskVol %s' % self.maskVol.get().getFileName()
         if self.mask.get() is not None:
             args += ' --mask %s' % self.mask.get().getFileName()
+        else:
+            involdim = vol.getDim()
+            fnDescr = self._getExtraPath("mask.descr")
+            fhDescr = open(fnDescr, 'w')
+            fhDescr.write("%d %d %d 0" % (involdim[0], involdim[1], involdim[2]))
+            fhDescr.close()
+            mskKeep = self._getExtraPath("mask_keep.mrc")
+            args_mask_keep = "-i %s -o %s" % (fnDescr, mskKeep)
+            self.runJob("xmipp_phantom_create", args_mask_keep, numberOfMpi=1)
+            args_imageheader2 = "-i %s --sampling_rate %f" % (mskKeep, vol.getSamplingRate())
+            self.runJob("xmipp_image_header", args_imageheader2, numberOfMpi=1)
+            args += ' --mask %s --subAll' % mskKeep
         self.runJob(program, args)
 
     def createOutputStep(self):
diff --git a/xmipp3/protocols/protocol_validate_fscq.py b/xmipp3/protocols/protocol_validate_fscq.py
index 8dde147c4..504670922 100644
--- a/xmipp3/protocols/protocol_validate_fscq.py
+++ b/xmipp3/protocols/protocol_validate_fscq.py
@@ -25,7 +25,7 @@
 # *
 # **************************************************************************
 
-
+import os
 from pyworkflow import VERSION_3_0
 from pyworkflow.protocol.params import (PointerParam, BooleanParam,
                                         IntParam, FileParam, FloatParam)
@@ -34,8 +34,9 @@
 from pwem.protocols import ProtAnalysis3D
 from pyworkflow.utils import getExt
 from pyworkflow.object import (Float, Integer)
-from pwem.objects import Volume
+from pwem.objects import AtomStruct
 from pwem.convert import Ccp4Header
+from pwem.convert.atom_struct import toPdb, toCIF, AtomicStructHandler, addScipionAttribute
 import xmipp3
 
 
@@ -48,7 +49,7 @@
 RESTA_FILE_MRC = 'diferencia.map'
 PDB_VALUE_FILE = 'pdb_fsc-q.pdb'
 MASK_FILE_MRC = 'mask.map'
-MASK_FILE = 'mask.vol' 
+MASK_FILE = 'mask.vol'
 FN_VOL = 'vol.map'
 FN_HALF1 = 'half1.map'
 FN_HALF2 = 'half2.map'
@@ -56,65 +57,73 @@
 MD2_MEANS = 'params2.xmd'
 RESTA_FILE_NORM = 'diferencia_norm.map'
 PDB_NORM_FILE = 'pdb_fsc-q_norm.pdb'
+OUTPUT_CIF = 'fscq_struct.cif'
 
 
 class XmippProtValFit(ProtAnalysis3D):
-    """    
+    """
     The protocol assesses the quality of the fit.
     """
     _label = 'validate fsc-q'
     _lastUpdateVersion = VERSION_3_0
-    
+    _ATTRNAME = 'fscq_score'
+    _OUTNAME = 'outputAtomStruct'
+    _possibleOutputs = {_OUTNAME: AtomStruct}
+
     def __init__(self, **args):
         ProtAnalysis3D.__init__(self, **args)
         self.stepsExecutionMode = 1
-    
+
     # --------------------------- DEFINE param functions ----------------------
     def _defineParams(self, form):
         form.addSection(label='Input')
-
-        form.addParam('inputVolume', PointerParam, pointerClass='Volume',
+        group = form.addGroup('Input')
+        group.addParam('inputVolume', PointerParam, pointerClass='Volume',
                       label="Input Volume", important=True,
                       help='Select a volume.')
-        
-#         form.addParam('inputPDB', PointerParam, pointerClass='PdbFile',
-#                       label="Refined PDB", important=True, )
-        form.addParam('inputPDB', FileParam,
-                      label="PDB File path", important=True,
+
+        group.addParam('fromFile', BooleanParam, default=False,
+                      label='Input PDB from file: ')
+        group.addParam('inputPDBObj', PointerParam, pointerClass='AtomStruct', allowsNull=True,
+                      label="Refined PDB: ", important=True, condition='not fromFile',
+                      help='Specify the desired input structure.')
+        group.addParam('inputPDB', FileParam,  condition='fromFile',
+                      label="PDB File path: ", important=True,
                       help='Specify a path to desired PDB structure.')
-        
-        form.addParam('pdbMap', PointerParam, pointerClass='Volume',
-                      label="Volume from PDB", allowsNull=True,
+
+        group.addParam('pdbMap', PointerParam, pointerClass='Volume',
+                      label="Volume from PDB: ", allowsNull=True,
                       help='Volume created from the PDB.'
                            ' The volume should be aligned with the reconstruction map.'
                            ' If the volume is not entered,' 
-                           ' it is automatically created from the PDB.')        
+                           ' it is automatically created from the PDB.')
 
-        form.addParam('inputMask', PointerParam, pointerClass='VolumeMask', 
+        group.addParam('inputMask', PointerParam, pointerClass='VolumeMask',
                       allowsNull=True,
-                      label="Soft Mask", 
+                      label="Soft Mask",
                       help='The mask determines which points are specimen'
                       ' and which are not. If the mask is not passed,' 
                       ' the method creates an automatic mask from the PDB.')
-        
-        form.addParam('box', IntParam, default=20,
+
+        group = form.addGroup('Parameters')
+        group.addParam('box', IntParam, default=20,
                       label="window size",
                       help='Kernel size (slidding window) for determining'
                       ' local resolution (pixels/voxels).')
-        
-        form.addParam('setOrigCoord', BooleanParam,
+
+        group.addParam('setOrigCoord', BooleanParam,
                       label="Set origin of coordinates",
                       help="Option YES:\nA new volume will be created with "
                            "the given ORIGIN of coordinates. ",
-                      default=False)       
-        
-        form.addParam('xcoor', FloatParam, default=0, condition='setOrigCoord',
+                      default=False)
+
+        group.addParam('xcoor', FloatParam, default=0, condition='setOrigCoord',
                       label="x", help="offset along x axis")
-        form.addParam('ycoor', FloatParam, default=0, condition='setOrigCoord',
+        group.addParam('ycoor', FloatParam, default=0, condition='setOrigCoord',
                       label="y", help="offset along y axis")
-        form.addParam('zcoor', FloatParam, default=0, condition='setOrigCoord',
+        group.addParam('zcoor', FloatParam, default=0, condition='setOrigCoord',
                       label="z", help="offset along z axis")
-        
+
         form.addParallelSection(threads=8, mpi=1)
 
     # --------------------------- INSERT steps functions --------------------------------------------
@@ -130,33 +139,40 @@ def _createFilenameTemplates(self):
                  RESTA_FILE_MRC: self._getExtraPath('diferencia.map'),
                  RESTA_FILE_NORM: self._getExtraPath('diferencia_norm.map'),
                  PDB_VALUE_FILE:  self._getExtraPath('pdb_fsc-q.pdb'),
-                 PDB_NORM_FILE: self._getExtraPath('pdb_fsc-q_norm.pdb'), 
-                 MASK_FILE_MRC : self._getExtraPath('mask.map'),  
-                 MASK_FILE: self._getTmpPath('mask.vol'), 
+                 PDB_NORM_FILE: self._getExtraPath('pdb_fsc-q_norm.pdb'),
+                 MASK_FILE_MRC : self._getExtraPath('mask.map'),
+                 MASK_FILE: self._getTmpPath('mask.vol'),
                  FN_VOL: self._getTmpPath("vol.map"),
                  FN_HALF1: self._getTmpPath("half1.map"),
-                 FN_HALF2: self._getTmpPath("half2.map"),  
-                 MD_MEANS: self._getExtraPath('params.xmd'),  
-                 MD2_MEANS: self._getExtraPath('params2.xmd')      
+                 FN_HALF2: self._getTmpPath("half2.map"),
+                 MD_MEANS: self._getExtraPath('params.xmd'),
+                 MD2_MEANS: self._getExtraPath('params2.xmd')
                  }
         self._updateFilenamesDict(myDict)
 
     def _insertAllSteps(self):
- 
-        self._createFilenameTemplates() 
+
+        self._createFilenameTemplates()
         input = self._insertFunctionStep('convertInputStep')
         id = []
         for i in range(2):
-            id.append(self._insertFunctionStep('runBlocresStep', i, prerequisites=[input]))   
-        input1 = self._insertFunctionStep('substractBlocresStep',prerequisites=id)    
-        input2 = self._insertFunctionStep('assignPdbStep', prerequisites=[input1])  
-        self._insertFunctionStep('createOutputStep', prerequisites=[input2])         
-           
+            id.append(self._insertFunctionStep('runBlocresStep', i, prerequisites=[input]))
+        input1 = self._insertFunctionStep('substractBlocresStep',prerequisites=id)
+        input2 = self._insertFunctionStep('assignPdbStep', prerequisites=[input1])
+        self._insertFunctionStep('createOutputStep', prerequisites=[input2])
+
 
     def convertInputStep(self):
+        """ Convert inputs to desired format."""
+        #Convert Input to pdb
+        if self.getInputStructFile().endswith('.pdb') or self.getInputStructFile().endswith('.ent'):
+            os.symlink(os.path.abspath(self.getInputStructFile()), self.getPDBFile())
+        else:
+            toPdb(self.getInputStructFile(), self.getPDBFile())
+
         """ Read the input volume."""
         self.volume = self.inputVolume.get()
-        
+
         """Read the Origin."""
         if self.setOrigCoord.get():
             self.x = self.xcoor.get()
@@ -165,116 +181,116 @@ def convertInputStep(self):
         else:
             self.x = 0
             self.y = 0
-            self.z = 0  
-                                       
+            self.z = 0
+
         self.sampling = self.volume.getSamplingRate()
-        self.origin = (-self.x*self.sampling, -self.y*self.sampling, -self.z*self.sampling) 
-        
+        self.origin = (-self.x*self.sampling, -self.y*self.sampling, -self.z*self.sampling)
+
         self.vol = self.volume.getFileName()
         self.half1, self.half2 = self.volume.getHalfMaps().split(',')
 
-        extVol = getExt(self.vol)        
+        extVol = getExt(self.vol)
         if (extVol == '.mrc') or (extVol == '.map'):
-            self.vol_xmipp = self.vol + ':mrc' 
+            self.vol_xmipp = self.vol + ':mrc'
         else:
-            self.vol_xmipp = self.vol     
- 
-        self.fnvol = self._getFileName(FN_VOL)           
+            self.vol_xmipp = self.vol
+
+        self.fnvol = self._getFileName(FN_VOL)
         self.fnvol1 = self._getFileName(FN_HALF1)
         self.fnvol2 = self._getFileName(FN_HALF2)
 
         Ccp4Header.fixFile(self.vol, self.fnvol, self.origin, self.sampling,
                         Ccp4Header.START)
         Ccp4Header.fixFile(self.half1, self.fnvol1, self.origin, self.sampling,
-                        Ccp4Header.START) 
+                        Ccp4Header.START)
         Ccp4Header.fixFile(self.half2, self.fnvol2, self.origin, self.sampling,
-                        Ccp4Header.START)        
-        
-        """Create map from PDB """         
-        if self.pdbMap.hasValue():   
+                        Ccp4Header.START)
+
+        """Create map from PDB """
+        if self.pdbMap.hasValue():
             pdbvolume = self.pdbMap.get()
             self.pdbvol = pdbvolume.getFileName()
-            
+
             self.pdbmap = self._getFileName(OUTPUT_PDBMRC_FILE)
             Ccp4Header.fixFile(self.pdbvol, self.pdbmap, self.origin, self.sampling,
-                        Ccp4Header.START) 
-      
-        else:         
-            """ Convert PDB to Map """           
+                        Ccp4Header.START)
+
+        else:
+            """ Convert PDB to Map """
             params = ' --centerPDB '
-            params += ' -v 0 '    
-            params += ' --sampling %f' % self.inputVolume.get().getSamplingRate()        
+            params += ' -v 0 '
+            params += ' --sampling %f' % self.inputVolume.get().getSamplingRate()
             params += ' --size %d' % self.inputVolume.get().getXDim()
-            params += ' -i %s' % self.inputPDB.get()        
+            params += ' -i %s' % self.getPDBFile()
             params += ' -o %s' % self._getFileName(OUTPUT_PDBVOL_FILE)
-            self.runJob('xmipp_volume_from_pdb', params)            
-    
-            """ Align pdbMap to reconstruction Map """  
-              
-            params = ' --i1 %s' % self.vol_xmipp        
+            self.runJob('xmipp_volume_from_pdb', params)
+
+            """ Align pdbMap to reconstruction Map """
+
+            params = ' --i1 %s' % self.vol_xmipp
             params += ' --i2 %s' % self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'
-            params += ' --local --apply'  
-            params += ' %s' % self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'                    
-            self.runJob('xmipp_volume_align', params)  
-            
+            params += ' --local --apply'
+            params += ' %s' % self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'
+            self.runJob('xmipp_volume_align', params)
+
             """ convert align vol to mrc format """
-            
-            self.pdbvol = self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol' 
+
+            self.pdbvol = self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'
             self.pdbmap = self._getFileName(OUTPUT_PDBMRC_FILE)
             Ccp4Header.fixFile(self.pdbvol, self.pdbmap, self.origin, self.sampling,
-                        Ccp4Header.START) 
-        
-        
-        """ Create a mask"""               
+                        Ccp4Header.START)
+
+
+        """ Create a mask"""
         if self.inputMask.hasValue():
             self.maskIn = self.inputMask.get().getFileName()
-            self.maskFn = self._getFileName(MASK_FILE_MRC)           
+            self.maskFn = self._getFileName(MASK_FILE_MRC)
             Ccp4Header.fixFile(self.maskIn, self.maskFn, self.origin, self.sampling,
                         Ccp4Header.START)
-           
-            self.mask_xmipp = self.maskFn + ':mrc' 
-            
+
+            self.mask_xmipp = self.maskFn + ':mrc'
+
         else:
-            
+
             self.maskFn = self._getFileName(MASK_FILE_MRC)
-            self.mask_xmipp = self._getFileName(MASK_FILE) 
-            
+            self.mask_xmipp = self._getFileName(MASK_FILE)
+
             if (not self.pdbMap.hasValue()):
-                            
-                params = ' -i %s' % self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'          
+
+                params = ' -i %s' % self._getFileName(OUTPUT_PDBVOL_FILE)+'.vol'
                 params += ' -o %s' % self.mask_xmipp
-                params += ' --select below 0.02 --substitute binarize'                   
-                self.runJob('xmipp_transform_threshold', params) 
-                 
-                params = ' -i %s' % self.mask_xmipp        
+                params += ' --select below 0.02 --substitute binarize'
+                self.runJob('xmipp_transform_threshold', params)
+
+                params = ' -i %s' % self.mask_xmipp
                 params += ' -o %s' % self.mask_xmipp
-                params += ' --binaryOperation dilation --size 3'                    
-                self.runJob('xmipp_transform_morphology', params) 
-                
+                params += ' --binaryOperation dilation --size 3'
+                self.runJob('xmipp_transform_morphology', params)
+
                 """ convert mask.vol to mrc format """
-                self.maskFn = self._getFileName(MASK_FILE_MRC)           
+                self.maskFn = self._getFileName(MASK_FILE_MRC)
                 Ccp4Header.fixFile(self.mask_xmipp, self.maskFn, self.origin, self.sampling,
                         Ccp4Header.START)
-                
+
             else:
-                
-                """ create mask from pdbMap """                
-                params = ' -i %s' % self._getFileName(OUTPUT_PDBMRC_FILE)+':mrc'          
+
+                """ create mask from pdbMap """
+                params = ' -i %s' % self._getFileName(OUTPUT_PDBMRC_FILE)+':mrc'
                 params += ' -o %s' % self.mask_xmipp
-                params += ' --select below 0.02 --substitute binarize'                    
-                self.runJob('xmipp_transform_threshold', params) 
-                 
-                params = ' -i %s' % self.mask_xmipp        
+                params += ' --select below 0.02 --substitute binarize'
+                self.runJob('xmipp_transform_threshold', params)
+
+                params = ' -i %s' % self.mask_xmipp
                 params += ' -o %s' % self.mask_xmipp
-                params += ' --binaryOperation dilation --size 3'                    
-                self.runJob('xmipp_transform_morphology', params) 
-                 
+                params += ' --binaryOperation dilation --size 3'
+                self.runJob('xmipp_transform_morphology', params)
+
                 """ convert mask.vol to mrc format """
-                self.maskFn = self._getFileName(MASK_FILE_MRC)           
+                self.maskFn = self._getFileName(MASK_FILE_MRC)
                 Ccp4Header.fixFile(self.mask_xmipp, self.maskFn, self.origin, self.sampling,
                         Ccp4Header.START)
-                       
-                              
+
+
     def runBlocresStep(self, i):
         # Local import to prevent discovery errors
         import bsoft
@@ -283,7 +299,7 @@ def runBlocresStep(self, i):
 
             """ Calculate FSC map-PDB """
 
-            params = ' -criterio FSC -nofill -smooth -pad 1 '
+            params = '  -nofill -smooth -pad 1 '
             params += ' -cutoff 0.67'
             params += ' -maxresolution 0.5 '
             params += ' -step 1 '
@@ -291,86 +307,84 @@ def runBlocresStep(self, i):
             params += ' -sampling %f,%f,%f' % (self.inputVolume.get().getSamplingRate(),
                                                self.inputVolume.get().getSamplingRate(),
                                                self.inputVolume.get().getSamplingRate())
-#             params += ' -origin %f,%f,%f' % ((self.shifts[0], self.shifts[1], self.shifts[2]))          
+#             params += ' -origin %f,%f,%f' % ((self.shifts[0], self.shifts[1], self.shifts[2]))
             params += ' -Mask %s' % self.maskFn
             params += ' %s  %s' % (self.fnvol, self._getFileName(OUTPUT_PDBMRC_FILE))
             params += ' %s' % self._getFileName(BLOCRES_AVG_FILE)
 
-            self.runJob(bsoft.Plugin.getProgram('blocres', bsoftVersion=bsoft.V1_9_0), params,
-                        env=bsoft.Plugin.getEnviron(bsoftVersion=bsoft.V1_9_0))
+            self.runJob(bsoft.Plugin.getProgram('blocres'), params)
         else:
 
             """ Calculate FSC half1-half2 """
 
-            params = ' -criterio FSC -nofill -smooth -pad 1 '
+            params = '  -nofill -smooth -pad 1 '
             params += ' -cutoff 0.5'
             params += ' -maxresolution 0.5 '
             params += ' -step 1 '
             params += ' -box %d ' % self.box.get()
             params += ' -sampling %f,%f,%f' % (self.inputVolume.get().getSamplingRate(),
                                                self.inputVolume.get().getSamplingRate(),
-                                               self.inputVolume.get().getSamplingRate()) 
-#             params += ' -origin %f,%f,%f' % ((self.shifts[0], self.shifts[1], self.shifts[2])) 
+                                               self.inputVolume.get().getSamplingRate())
+#             params += ' -origin %f,%f,%f' % ((self.shifts[0], self.shifts[1], self.shifts[2]))
             params += ' -Mask %s' % self.maskFn
             params += ' %s  %s' % (self.fnvol1, self.fnvol2)
             params += ' %s' % self._getFileName(BLOCRES_HALF_FILE)
 
-            self.runJob(bsoft.Plugin.getProgram('blocres', bsoftVersion=bsoft.V1_9_0), params,
-                        env=bsoft.Plugin.getEnviron(bsoftVersion=bsoft.V1_9_0))
+            self.runJob(bsoft.Plugin.getProgram('blocres'), params)
 
     def substractBlocresStep(self):
-        
-        params = ' -i %s' % self._getFileName(BLOCRES_AVG_FILE)+':mrc'  
-        params += ' --minus %s' % self._getFileName(BLOCRES_HALF_FILE)+':mrc'     
-        params += ' -o %s ' % self._getFileName(RESTA_FILE)       
+
+        params = ' -i %s' % self._getFileName(BLOCRES_AVG_FILE)+':mrc'
+        params += ' --minus %s' % self._getFileName(BLOCRES_HALF_FILE)+':mrc'
+        params += ' -o %s ' % self._getFileName(RESTA_FILE)
         self.runJob('xmipp_image_operate', params)
-        
-        Ccp4Header.fixFile(self._getFileName(RESTA_FILE), self._getFileName(RESTA_FILE_MRC), 
+
+        Ccp4Header.fixFile(self._getFileName(RESTA_FILE), self._getFileName(RESTA_FILE_MRC),
                            self.origin, self.sampling, Ccp4Header.START)
-        
-        """Diveded by resolution"""       
+
+        """Diveded by resolution"""
         Vx = xmipp3.Image(self._getFileName(RESTA_FILE))
         V=Vx.getData()
         Vmask = xmipp3.Image(self._getFileName(MASK_FILE_MRC)+':mrc').getData()
         Vres = xmipp3.Image(self._getFileName(BLOCRES_HALF_FILE)+':mrc').getData()
         Vt = V
         Zdim, Ydim, Xdim = V.shape
-              
+
         for z in range(0,Zdim):
             for y in range(0,Ydim):
                 for x in range(0,Xdim):
-                    if (Vmask[z,y,x] > 0.001 and Vres[z,y,x]>0.001): 
+                    if (Vmask[z,y,x] > 0.001 and Vres[z,y,x]>0.001):
                         Vt[z,y,x] = (V[z,y,x]/Vres[z,y,x])
-        Vx.setData(Vt) 
+        Vx.setData(Vt)
         Vx.write(self._getFileName(RESTA_FILE_NORM))
-        Ccp4Header.fixFile(self._getFileName(RESTA_FILE_NORM), self._getFileName(RESTA_FILE_NORM), 
-                           self.origin, self.sampling, Ccp4Header.START) 
-                    
-        
+        Ccp4Header.fixFile(self._getFileName(RESTA_FILE_NORM), self._getFileName(RESTA_FILE_NORM),
+                           self.origin, self.sampling, Ccp4Header.START)
+
+
     def assignPdbStep(self):
-                 
-        params = ' --pdb %s ' % self.inputPDB.get()  
-        params += ' --vol %s ' % self._getFileName(RESTA_FILE_MRC) 
-        params += ' --mask %s ' % self.mask_xmipp         
-        params += ' -o %s ' % self._getFileName(PDB_VALUE_FILE)    
+
+        params = ' --pdb %s ' % self.getPDBFile()
+        params += ' --vol %s ' % self._getFileName(RESTA_FILE_MRC)
+        params += ' --mask %s ' % self.mask_xmipp
+        params += ' -o %s ' % self._getFileName(PDB_VALUE_FILE)
         params += ' --sampling %f' % self.inputVolume.get().getSamplingRate()
         params += ' --origin %f %f %f' %(self.x, self.y, self.z)
-        params += ' --radius 0.8'  
-        params += ' --md %s' % self._getFileName(MD_MEANS)     
-        self.runJob('xmipp_pdb_label_from_volume', params)   
-        
+        params += ' --radius 0.8'
+        params += ' --md %s' % self._getFileName(MD_MEANS)
+        self.runJob('xmipp_pdb_label_from_volume', params)
+
         """Diveded by resolution"""
-        params = ' --pdb %s ' % self.inputPDB.get()  
-        params += ' --vol %s ' % self._getFileName(RESTA_FILE_NORM) 
-        params += ' --mask %s ' % self.mask_xmipp         
-        params += ' -o %s ' % self._getFileName(PDB_NORM_FILE)    
+        params = ' --pdb %s ' % self.getPDBFile()
+        params += ' --vol %s ' % self._getFileName(RESTA_FILE_NORM)
+        params += ' --mask %s ' % self.mask_xmipp
+        params += ' -o %s ' % self._getFileName(PDB_NORM_FILE)
         params += ' --sampling %f' % self.inputVolume.get().getSamplingRate()
         params += ' --origin %f %f %f' %(self.x, self.y, self.z)
-        params += ' --radius 0.8' 
-        params += ' --md %s' % self._getFileName(MD2_MEANS) 
-        self.runJob('xmipp_pdb_label_from_volume', params)    
+        params += ' --radius 0.8'
+        params += ' --md %s' % self._getFileName(MD2_MEANS)
+        self.runJob('xmipp_pdb_label_from_volume', params)
 
-    def _getMetrics(self):
+    def _setMetrics(self):
         """ Internal method to compute some metrics. """
         # mean values of FSC-Q
 
@@ -410,28 +424,48 @@ def _getMetrics(self):
         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)
-        }
+        self.mean = Float(mean)
+        self.meanA = Float(meanA)
+        self.mean2 = Float(mean2)
+        self.meanA2 = Float(meanA2)
+        self.total_atom = Integer(total_atom)
+        self.fscq_greater = Integer(fscq_greater)
+        self.fscq_less = Integer(fscq_less)
+        self.porc_greater = Float(porc_greater)
+        self.porc_less = Float(porc_less)
+        self._store()
+
+    def getFscqAttrDic(self):
+        fscqDic = {}
+        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')):
+                    resNumber, resChain, atomName = lin[22:26].strip(), lin[21].strip(), lin[12:16].strip()
+                    resId = '{}:{}@{}'.format(resChain, resNumber, atomName)
+                    fscq_atom = lin[54:60].strip()
+                    fscqDic[resId] = fscq_atom
+        return fscqDic
 
     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._setMetrics()
+        fscqDic = self.getFscqAttrDic()
+
+        AS = self.getInputStruct()
+        ASH = AtomicStructHandler()
+        inpAS = toCIF(AS.getFileName(), self._getTmpPath('inputStruct.cif'))
+        cifDic = ASH.readLowLevel(inpAS)
+        cifDic = addScipionAttribute(cifDic, fscqDic, self._ATTRNAME, recipient='atoms')
+        ASH._writeLowLevel(self._getPath(OUTPUT_CIF), cifDic)
+
+        outAS = AS.clone()
+        outAS.setFileName(self._getPath(OUTPUT_CIF))
+        outAS.setVolume(self.inputVolume.get())
+
+        self._defineOutputs(outputAtomStruct=outAS)
+        if not self.fromFile:
+            self._defineTransformRelation(self.inputPDBObj, outAS)
 
-        self._defineOutputs(fscq_Volume=volume, **metrics)
-        self._defineTransformRelation(self.inputVolume, volume)
-        
 
     # --------------------------- INFO functions ------------------------------
 
@@ -441,30 +475,30 @@ def _methods(self):
             messages.append(
                 'Information about the method/article in ' + VALIDATE_METHOD_URL)
         return messages
-    
+
     def _summary(self):
         summary = []
         summary.append("Deviation from the signal of the Half Maps")
         if self.hasAttribute('mean'):
             summary.append("Mean FSC-Q: %.2f" % (self.mean.get()))
-        if self.hasAttribute('meanA'):    
+        if self.hasAttribute('meanA'):
             summary.append("Absotute Mean FSC-Q: %.2f" % (self.meanA.get()))
-             
-        summary.append(" ")            
+
+        summary.append(" ")
         summary.append("Deviation from the signal of the Half Maps divided by local resolution")
         if self.hasAttribute('mean2'):
             summary.append("Mean FSC-Qr: %.2f" % (self.mean2.get()))
-        if self.hasAttribute('meanA2'):    
-            summary.append("Absotute Mean FSC-Qr: %.2f" % (self.meanA2.get()))      
-             
-        summary.append("------------------------------------------")   
-        if self.hasAttribute('total_atom'):                    
+        if self.hasAttribute('meanA2'):
+            summary.append("Absotute Mean FSC-Qr: %.2f" % (self.meanA2.get()))
+
+        summary.append("------------------------------------------")
+        if self.hasAttribute('total_atom'):
             summary.append("Total number of atoms analyzed: %d" % (self.total_atom.get()))
-        if (self.hasAttribute('fscq_greater') and self.hasAttribute('porc_greater')):            
-            summary.append("Number of atoms with FSC-Q>0.5: %d.  Percentage of total: %.2f." 
+        if (self.hasAttribute('fscq_greater') and self.hasAttribute('porc_greater')):
+            summary.append("Number of atoms with FSC-Q>0.5: %d.  Percentage of total: %.2f."
                        % (self.fscq_greater.get(), self.porc_greater.get()))
-        if (self.hasAttribute('fscq_less') and self.hasAttribute('porc_less')):             
-            summary.append("Number of atoms with FSC-Q<-0.5: %d.  Percentage of total: %.2f." 
+        if (self.hasAttribute('fscq_less') and self.hasAttribute('porc_less')):
+            summary.append("Number of atoms with FSC-Q<-0.5: %d.  Percentage of total: %.2f."
                        % (self.fscq_less.get(), self.porc_less.get()))
         return summary
 
@@ -481,6 +515,11 @@ def _validate(self):
                 errors.append("Input Volume needs to have half maps. "
                 "If you have imported the volume, be sure to import the half maps.")
 
+        if self.fromFile and not self.inputPDB.get():
+            errors.append('You have to provide a PDB file as input')
+        elif not self.fromFile and not self.inputPDBObj.get():
+            errors.append('You have to provide a PDB object as input')
+
         try:
             import bsoft
             if bsoft.__version__ in ["3.0.0", "3.0.1", "3.0.4"]:
@@ -488,9 +527,24 @@ def _validate(self):
                               " You have %s. Update it using the plugin manager or command line" % bsoft.__version__)
         except Exception as e:
             errors.append("This protocol requires bsoft plugin 3.0.5 or above to run. Update it using the plugin manager or command line")
-            
-        return errors    
+
+        return errors
 
     def _citations(self):
         return ['Ramirez-Aportela 2020']
 
+    def getInputStructFile(self):
+        if self.fromFile:
+            return self.inputPDB.get()
+        else:
+            return self.inputPDBObj.get().getFileName()
+
+    def getPDBFile(self):
+        return self._getExtraPath('inputPDB.pdb')
+
+    def getInputStruct(self):
+        if self.fromFile:
+            return AtomStruct(filename=self.inputPDB.get())
+        else:
+            return self.inputPDBObj.get()
+
diff --git a/xmipp3/viewers/__init__.py b/xmipp3/viewers/__init__.py
index 367bc596e..f425ae317 100644
--- a/xmipp3/viewers/__init__.py
+++ b/xmipp3/viewers/__init__.py
@@ -66,4 +66,4 @@
 
 # TODO: Import from continuousflex the modules needed to create clusters (soft dependency)
 from .viewer_angular_alignment_sph import XmippAngularAlignmentSphViewer
-
+from .viewer_analyze_local_ctf import XmippAnalyzeLocalCTFViewer
diff --git a/xmipp3/viewers/viewer.py b/xmipp3/viewers/viewer.py
index 0013a9a8a..ca7ab84ba 100644
--- a/xmipp3/viewers/viewer.py
+++ b/xmipp3/viewers/viewer.py
@@ -46,6 +46,7 @@
 from xmipp3.protocols import XmippProtCTFMicrographs
 from xmipp3.protocols import XmippProtValidateNonTilt
 from xmipp3.protocols import XmippProtMultiRefAlignability
+from xmipp3.protocols import XmippProtAngularGraphConsistency
 from xmipp3.protocols import XmippProtAssignmentTiltPair
 from xmipp3.protocols import XmippProtMovieGain
 from xmipp3.protocols import XmippProtDeepDenoising
@@ -73,6 +74,7 @@ class XmippViewer(DataViewer):
                 XmippProtValidateNonTilt,
                 XmippProtAssignmentTiltPair,
                 XmippProtMultiRefAlignability,
+                XmippProtAngularGraphConsistency,
                 XmippProtMovieGain,
                 XmippProtDeepDenoising,
                 XmippProtParticleBoxsize
@@ -324,6 +326,19 @@ def _visualize(self, obj, **kwargs):
                                emlib.MDL_SCORE_BY_ALIGNABILITY_ACCURACY,
                                marker='.', markersize=.55, color='red', linestyle='')
             self._views.append(plotter)
+    
+        elif issubclass(cls, XmippProtAngularGraphConsistency):
+            fn = obj.outputParticles.getFileName()
+            labels = ('id enabled _index _filename _xmipp_assignedDirRefCC '
+                      '_xmipp_maxCCprevious _xmipp_graphCCPrevious _xmipp_distance2MaxGraphPrevious '
+                      '_xmipp_maxCC _xmipp_graphCC _xmipp_distance2MaxGraph')
+            labelRender = "_filename"
+            self._views.append(ObjectView(self._project, obj.outputParticles.strId(), fn,
+                                          viewParams={ORDER: labels,
+                                                      VISIBLE: labels,
+                                                      SORT_BY: '_xmipp_assignedDirRefCC desc',
+                                                      RENDER: labelRender,
+                                                      MODE: MODE_MD}))
 
         elif issubclass(cls, XmippProtExtractParticlesPairs):
             self._visualize(obj.outputParticlesTiltPair)
@@ -333,6 +348,7 @@ def _visualize(self, obj, **kwargs):
             DataViewer._visualize(self, obj, **kwargs)
 
         return self._views
+            
 
     def getCTFViews(self, ctfSet):
         # This could be used by any CTF viewer to show CTF plus, phaseShift plot
diff --git a/xmipp3/viewers/viewer_validate_fscq.py b/xmipp3/viewers/viewer_validate_fscq.py
index 4368b3e15..9ca459510 100644
--- a/xmipp3/viewers/viewer_validate_fscq.py
+++ b/xmipp3/viewers/viewer_validate_fscq.py
@@ -35,7 +35,7 @@
 from pyworkflow.utils import getExt, removeExt
 from os.path import abspath
 from pwem.viewers import (LocalResolutionViewer, EmPlotter, ChimeraView,
-                          DataView)
+                          DataView, ChimeraAttributeViewer)
 from pyworkflow.gui import *
 from pyworkflow.protocol.params import (LabelParam, EnumParam)
 from pyworkflow.viewer import ProtocolViewer, DESKTOP_TKINTER
@@ -47,7 +47,7 @@
                                                      PDB_NORM_FILE)
 
 
-class XmippProtValFitViewer(LocalResolutionViewer):
+class XmippProtValFitViewer(LocalResolutionViewer, ChimeraAttributeViewer):
     """
     Visualization tools for validation fsc-q.
     
@@ -62,12 +62,12 @@ class XmippProtValFitViewer(LocalResolutionViewer):
     
     def __init__(self, *args, **kwargs):
         ProtocolViewer.__init__(self, *args, **kwargs)
-    
+
     def _defineParams(self, form):
         self._env = os.environ.copy()
         form.addSection(label='FSC-Q results')
         
-        group = form.addGroup('Visualization in Chimera')
+        group = form.addGroup('Chimera visualization')
         
         group.addParam('displayVolume', LabelParam,
                       important=True,
@@ -105,9 +105,15 @@ def _defineParams(self, form):
                       label='Amino acids with low resolvability',   
                       help='Amino acids that have atoms with FSC-Q > 1'
                       ' are determined. It is suggested that these amino acids '
-                      ' be re-checked in order to improve the fit.')       
-            
-        
+                      ' be re-checked in order to improve the fit.')
+
+        super()._defineParams(form)
+        from pwem.wizards.wizard import ColorScaleWizardBase
+        group = form.addGroup('Color settings')
+        ColorScaleWizardBase.defineColorScaleParams(group, defaultLowest=-3, defaultHighest=3, defaultIntervals=21,
+                                                    defaultColorMap='RdBu_r')
+
+
     def _getVisualizeDict(self):
         self.protocol._createFilenameTemplates()
         visualizeDict = {'displayVolume': self._visualize_vol,
@@ -116,6 +122,7 @@ def _getVisualizeDict(self):
                          'displayNormPDB': self._visualize_norm_pdb,
                          'calculateFscqNeg': self._statistics,
                          'calculateFscqPos': self._statistics}
+        visualizeDict.update(ChimeraAttributeViewer._getVisualizeDict(self))
         return visualizeDict     
     
     def _create_legend(self, scale):
diff --git a/xmipp3/xmipp_logo.png b/xmipp3/xmipp_logo.png
index 084bfe97a..1ada98399 100644
Binary files a/xmipp3/xmipp_logo.png and b/xmipp3/xmipp_logo.png differ


Test patch:
diff --git a/xmipp3/tests/test_protocol_angular_graph_consistency.py b/xmipp3/tests/test_protocol_angular_graph_consistency.py
new file mode 100644
index 000000000..2cbf93609
--- /dev/null
+++ b/xmipp3/tests/test_protocol_angular_graph_consistency.py
@@ -0,0 +1,99 @@
+# ***************************************************************************
+# * Authors:     Jeison Méndez García (jmendez@utp.edu.co)
+# *                based on TestMultireferenceAlignability (Javier Vargas, J.M. De la Rosa Trevin)
+# *                 
+# *
+# * This program is free software; you can redistribute it and/or modify
+# * it under the terms of the GNU General Public License as published by
+# * the Free Software Foundation; either version 2 of the License, or
+# * (at your option) any later version.
+# *
+# * This program is distributed in the hope that it will be useful,
+# * but WITHOUT ANY WARRANTY; without even the implied warranty of
+# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# * GNU General Public License for more details.
+# *
+# * You should have received a copy of the GNU General Public License
+# * along with this program; if not, write to the Free Software
+# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+# * 02111-1307  USA
+# *
+# *  All comments concerning this program package may be sent to the
+# *  e-mail address 'scipion@cnb.csic.es'
+# ***************************************************************************/
+
+from pyworkflow.tests import BaseTest, setupTestProject, DataSet
+from pwem.protocols import ProtImportParticles, ProtImportVolumes, ProtSubSet
+
+from xmipp3.protocols import XmippProtAngularGraphConsistency
+
+
+class TestAngularGraphConsistency(BaseTest):
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+        cls.dsRelion = DataSet.getDataSet('relion_tutorial')
+
+    def checkOutput(self, prot, outName, conditions=[]):
+        """ Check that an ouput was generated and
+        the condition is valid.
+        """
+        o = getattr(prot, outName, None)
+        locals()[outName] = o
+        self.assertIsNotNone(o, "Output: %s is None" % outName)
+        for c in conditions:
+            self.assertTrue(eval(c), 'Condition failed: ' + c)
+
+    def importVolumes(self):
+        importProt = self.newProtocol(ProtImportVolumes,
+                                objLabel='import volume',
+                                filesPath=self.dsRelion.getFile('volumes/reference.mrc'),
+                                samplingRate=7.08)
+        self.launchProtocol(importProt)
+
+        return importProt
+
+    def importParticles(self, numberOfParticles, path):
+        """ Import an EMX file with Particles and defocus
+        """
+        importProt = self.newProtocol(ProtImportParticles,
+                                 objLabel='from relion (auto-refine 3d)',
+                                 importFrom=ProtImportParticles.IMPORT_FROM_RELION,
+                                 starFile=self.dsRelion.getFile(path),
+                                 magnification=10000,
+                                 samplingRate=7.08,
+                                 haveDataBeenPhaseFlipped=True
+                                 )
+        self.launchProtocol(importProt)
+        self.checkOutput(importProt, 'outputParticles',
+                         ['outputParticles.hasAlignmentProj()',
+                          'outputParticles.isPhaseFlipped()'])
+        # subset to speed-up following processes
+        subsetProt = self.newProtocol(ProtSubSet,
+                                      objLabel='subset of particles',
+                                      chooseAtRandom=True,
+                                      nElements=numberOfParticles)
+        subsetProt.inputFullSet.set(importProt.outputParticles)
+
+        self.launchProtocol(subsetProt)
+        self.checkOutput(subsetProt, 'outputParticles',
+                         ['outputParticles.getSize() == %d'
+                          % numberOfParticles])
+
+        return subsetProt
+
+    def test_validate(self):
+        protImportVols = self.importVolumes()
+        pathToFile = 'import/case2/relion_it015_data.star'
+        importParticlesProt = self.importParticles(2000,pathToFile)
+
+        validateProt = self.newProtocol(XmippProtAngularGraphConsistency,
+                                objLabel='angular graph consistency',
+                                angularSampling=3,
+                                maximumTargetResolution=10,
+                                numberOfMpi=4, numberOfThreads=1)
+        validateProt.inputParticles.set(importParticlesProt.outputParticles)
+        validateProt.inputVolume.set(protImportVols.outputVolume)
+
+        self.launchProtocol(validateProt)
+        self.checkOutput(validateProt, 'outputParticles')       
diff --git a/xmipp3/tests/test_protocol_validate_fscq.py b/xmipp3/tests/test_protocol_validate_fscq.py
new file mode 100644
index 000000000..224fb072f
--- /dev/null
+++ b/xmipp3/tests/test_protocol_validate_fscq.py
@@ -0,0 +1,82 @@
+# **************************************************************************
+# *
+# * Authors:    Daniel Del Hoyo Gomez (ddelhoyo@cnb.csic.es)
+# *
+# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
+# *
+# * This program is free software; you can redistribute it and/or modify
+# * it under the terms of the GNU General Public License as published by
+# * the Free Software Foundation; either version 2 of the License, or
+# * (at your option) any later version.
+# *
+# * This program is distributed in the hope that it will be useful,
+# * but WITHOUT ANY WARRANTY; without even the implied warranty of
+# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# * GNU General Public License for more details.
+# *
+# * You should have received a copy of the GNU General Public License
+# * along with this program; if not, write to the Free Software
+# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+# * 02111-1307  USA
+# *
+# *  All comments concerning this program package may be sent to the
+# *  e-mail address 'scipion@cnb.csic.es'
+# *
+# **************************************************************************
+
+from pyworkflow.tests import BaseTest, setupTestProject, DataSet
+from pwem.protocols import ProtImportPdb, ProtImportVolumes
+from xmipp3.protocols import XmippProtValFit
+
+class TestValidateFSCQ(BaseTest):
+  @classmethod
+  def setUpClass(cls):
+    cls.ds = DataSet.getDataSet('model_building_tutorial')
+
+    setupTestProject(cls)
+    cls._runImportPDB()
+    cls._runImportVolume()
+
+  @classmethod
+  def _runImportPDB(cls):
+    protImportPDB = cls.newProtocol(
+      ProtImportPdb,
+      inputPdbData=1,
+      pdbFile=cls.ds.getFile('PDBx_mmCIF/5ni1.pdb'))
+    cls.launchProtocol(protImportPDB)
+    cls.protImportPDB = protImportPDB
+
+  @classmethod
+  def _runImportVolume(cls):
+    args = {'filesPath': cls.ds.getFile(
+      'volumes/emd_3488.map'),
+      'samplingRate': 1.05,
+      'setHalfMaps': True,
+      'half1map': cls.ds.getFile('volumes/emd_3488_Noisy_half1.vol'),
+      'half2map': cls.ds.getFile('volumes/emd_3488_Noisy_half2.vol'),
+      'setOrigCoord': True,
+      'x': 0.0,
+      'y': 0.0,
+      'z': 0.0
+    }
+    protImportVolume = cls.newProtocol(ProtImportVolumes, **args)
+    cls.launchProtocol(protImportVolume)
+    cls.protImportVolume = protImportVolume
+
+  def _runFSCQ(self):
+    protFSCQ = self.newProtocol(
+      XmippProtValFit,
+      fromFile=False,
+      inputPDBObj=self.protImportPDB.outputPdb,
+      inputVolume=self.protImportVolume.outputVolume)
+
+    self.launchProtocol(protFSCQ)
+    pdbOut = getattr(protFSCQ, 'outputAtomStruct', None)
+    self.assertIsNotNone(pdbOut)
+
+  def testFSCQ(self):
+    self._runFSCQ()
+
+
+
+
diff --git a/xmipp3/tests/test_protocols_realignment_classes.py b/xmipp3/tests/test_protocols_realignment_classes.py
index b95d6c7f0..979b6d606 100644
--- a/xmipp3/tests/test_protocols_realignment_classes.py
+++ b/xmipp3/tests/test_protocols_realignment_classes.py
@@ -166,14 +166,11 @@ def test_pattern(self):
         if protRealing.isFailed():
             self.assertTrue(False)
         if not protRealing.hasAttribute('outputClasses') or not \
-                protRealing.hasAttribute('outputParticles') or not \
-                protRealing.hasAttribute('outputCoordinates'):
+                protRealing.hasAttribute('outputParticles'):
             self.assertTrue(False)
         if protRealing.outputClasses.getSize() != numClasses:
             self.assertTrue(False)
-        if protRealing.outputParticles.getSize() != \
-                protRealing.outputCoordinates.getSize():
-            self.assertTrue(False)
+
 
 
 
diff --git a/xmipp3/tests/test_protocols_xmipp_2d.py b/xmipp3/tests/test_protocols_xmipp_2d.py
index e9ef6756a..1afab6fc6 100644
--- a/xmipp3/tests/test_protocols_xmipp_2d.py
+++ b/xmipp3/tests/test_protocols_xmipp_2d.py
@@ -944,6 +944,39 @@ def test_cl2d(self):
         self.assertIsNotNone(protCL2DInitialCore.outputClasses, "There was a problem with CL2D with initial class and core analysis")
 
 
+class TestXmippCL2DCoreAnalysis(TestXmippBase):
+    """This class check if the protocol to core Analysis with CL2D in Xmipp works properly."""
+
+    @classmethod
+    def setUpClass(cls):
+        setupTestProject(cls)
+        TestXmippBase.setData('mda')
+        cls.protImport = cls.runImportParticles(cls.particlesFn, 3.5)
+
+    def test_cl2dCoreAnalysis(self):
+        print("Run CL2D core analysis")
+        # Run CL2D with random class and no core analysis
+        protCL2DRandomNoCore = self.newProtocol(XmippProtCL2D,
+                                                numberOfClasses=3, numberOfInitialClasses=1,
+                                                doCore=False, numberOfIterations=4, numberOfMpi=2)
+        protCL2DRandomNoCore.inputParticles.set(self.protImport.outputParticles)
+        protCL2DRandomNoCore.setObjLabel("CL2D with random class and no core analysis")
+        self.launchProtocol(protCL2DRandomNoCore)
+        self.assertIsNotNone(protCL2DRandomNoCore.outputClasses,
+                             "There was a problem with CL2D with random class and no core analysis")
+
+        # Run Core Analysis
+        protCL2DCoreAnalysis = self.newProtocol(XmippProtCoreAnalysis,
+                                                numberOfClasses=3, thZscore=1.5,
+                                                thPCAZscore=1.5, numberOfMpi=2)
+        protCL2DCoreAnalysis.inputClasses.set(protCL2DRandomNoCore.outputClasses)
+        protCL2DCoreAnalysis.setObjLabel("Core analysis")
+        self.launchProtocol(protCL2DCoreAnalysis)
+        self.assertIsNotNone(protCL2DCoreAnalysis.outputClasses_core,
+                             "There was a problem with CL2D Core Analysis")
+
+
+
 class TestXmippProtCL2DAlign(TestXmippBase):
     """This class check if the protocol to align particles in Xmipp works properly."""
     @classmethod
@@ -1165,6 +1198,10 @@ def setUpClass(cls):
         cls.protImportPart = cls.runImportParticles(cls.particlesFn, 3.5)
         cls.protImportAvgs = cls.runImportAverages(cls.particlesFn, 3.5)
         cls.protImportVol = cls.runImportVolume(cls.volumesFn, 3.5)
+        cls.protCreateMask = cls.newProtocol(XmippProtCreateMask3D, threshold=0.01)
+        cls.protCreateMask.inputVolume.set(cls.protImportVol.outputVolume)
+        cls.launchProtocol(cls.protCreateMask)
+        cls.assertIsNotNone(cls.protCreateMask.getFiles(), "There was a problem with the 3D mask")
         cls.protClassify = cls.runClassify(cls.protImportPart.outputParticles)
         cls.protProjMatch = cls.newProtocol(XmippProtProjMatch,
                                             doCTFCorrection=False,
@@ -1203,8 +1240,51 @@ def test_particles3(self):
         prot.inputSet.set(self.protProjMatch.outputParticles)
         prot.inputVolume.set(self.protImportVol.outputVolume)
         self.launchProtocol(prot)
-        self.assertIsNotNone(prot.reprojections, "There was a problem with Compare Reprojections from projections with angles")
+        self.assertIsNotNone(prot.reprojections,
+                             "There was a problem with Compare Reprojections from projections with angles")
 
+    def test_particles4(self):
+        print("Run Compare Reprojections from classes evaluating residuals")
+        prot = self.newProtocol(XmippProtCompareReprojections,
+                                symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
+        prot.inputSet.set(self.protClassify.outputClasses)
+        prot.inputVolume.set(self.protImportVol.outputVolume)
+        prot.maskVol.set(self.protCreateMask.outputMask)
+        self.launchProtocol(prot)
+        self.assertIsNotNone(prot.reprojections,
+                             "There was a problem with Compare Reprojections from classes evaluating residuals")
+
+    def test_particles5(self):
+        print("Run Compare Reprojections from averages evaluating residuals")
+        prot = self.newProtocol(XmippProtCompareReprojections,
+                                symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
+        prot.inputSet.set(self.protImportAvgs.outputAverages)
+        prot.inputVolume.set(self.protImportVol.outputVolume)
+        prot.maskVol.set(self.protCreateMask.outputMask)
+        self.launchProtocol(prot)
+        self.assertIsNotNone(prot.reprojections,
+                             "There was a problem with Compare Reprojections from averages evaluating residuals")
+
+    def test_particles6(self):
+        print("Run Compare Reprojections from projections with angles evaluating residuals")
+        prot = self.newProtocol(XmippProtCompareReprojections,
+                                symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
+        prot.inputSet.set(self.protProjMatch.outputParticles)
+        prot.inputVolume.set(self.protImportVol.outputVolume)
+        prot.maskVol.set(self.protCreateMask.outputMask)
+        self.launchProtocol(prot)
+        self.assertIsNotNone(prot.reprojections, "There was a problem with Compare Reprojections from projections"
+                                                 " with angles evaluating residuals")
+
+    def test_particles7(self):
+        print("Run Compare Reprojections from projections with angles evaluating residuals without user mask")
+        prot = self.newProtocol(XmippProtCompareReprojections,
+                                symmetryGroup="d6", numberOfMpi=5, doEvaluateResiduals=True)
+        prot.inputSet.set(self.protProjMatch.outputParticles)
+        prot.inputVolume.set(self.protImportVol.outputVolume)
+        self.launchProtocol(prot)
+        self.assertIsNotNone(prot.reprojections, "There was a problem with Compare Reprojections from projections"
+                                                 " with angles evaluating residuals without user mask")
 
 class TestXmippCreateGallery(TestXmippBase):
     """This class check if the protocol create gallery in Xmipp works properly."""