From bdf197845c3a06fbf6cd16f5de899e04b07b783f Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Wed, 13 Sep 2023 11:12:51 +0200
Subject: [PATCH 01/11] small fixes

---
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index f3c93369..7e2871be 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -53,7 +53,7 @@ def _defineParams(self, form):
         form.addSection(label='Input')
         form.addParam('inputSubtomos', PointerParam, pointerClass='SetOfSubTomograms', label="Subtomograms ",
                       help='Select the SetOfSubTomograms with transform matrix which will be subtracted.')
-        form.addParam('average', PointerParam, pointerClass='SetOfSubTomograms', label="Average subtomogram ",
+        form.addParam('average', PointerParam, pointerClass='SubTomogram', label="Average subtomogram ",
                       help='Select an average subtomogram to be subtracted.')
         form.addParam('maskBool', BooleanParam, label='Mask subtomograms?', default=True,
                       help='The mask are not mandatory but highly recommendable.')
@@ -150,7 +150,7 @@ def applyAlignStep(self):
 
     def subtractionStep(self):
         """Subtract reference to each of the subtomogram in the input Set"""
-        average = self.average.get().getFirstItem().getFileName()
+        average = self.average.get().getFileName()
         if average.endswith('.mrc'):
             average += ':mrc'
         resol = self.resol.get()
@@ -206,7 +206,7 @@ def _summary(self):
             if self.maskBool:
                 summary.append("Mask: %s" % self.mask.get().getFileName())
             if self.resol.get() != 0:
-                summary.append("Subtraction at resolution %f A" % self.resol.get())
+                summary.append("Subtraction at resolution %0.2f A" % self.resol.get())
         return summary
 
     def _methods(self):
@@ -217,7 +217,7 @@ def _methods(self):
             methods.append("Subtraction of average %s performed to %s subtomograms" %
                            (self.average.get().getFileName(), self.inputSubtomos.get().getSize()))
             if self.resol.get() != 0:
-                methods.append(" at resolution %f A" % self.resol.get())
+                methods.append(" at resolution %0.2f A" % self.resol.get())
         return methods
 
     def _validate(self):
@@ -229,8 +229,6 @@ def _validate(self):
             if not subtomo.hasTransform():
                 validateMsgs.append(
                     'Please provide subtomograms which have transformation matrix.')
-        if self.average.get().getSize() != 1:
-            validateMsgs.append('Average subtomogram must contain just 1 item.')
         return validateMsgs
 
     # --------------------------- UTLIS functions --------------------------------------------

From 7e56ca46487952daf964c7641ff6cbb43540fd66 Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Fri, 15 Sep 2023 11:05:41 +0200
Subject: [PATCH 02/11] changes for MPI version

---
 .../protocols/protocol_subtraction_subtomo.py   | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 7e2871be..7a0c05a1 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -26,7 +26,7 @@
 # *
 # **************************************************************************
 
-from pyworkflow import BETA
+from pyworkflow import BETA, UPDATED, NEW, PROD
 from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 
@@ -36,8 +36,11 @@
 from pwem.protocols import EMProtocol
 
 from tomo.protocols import ProtTomoBase
+from xmipp3.convert import writeSetOfParticles, readSetOfParticles
 from xmipp3.convert import xmippToLocation, writeSetOfVolumes, alignmentToRow
 
+OUTPUT = "output_subtomograms.xmd"
+INPUT_PARTICLES = "input_subtomograms.xmd"
 
 class XmippProtSubtractionSubtomo(EMProtocol, ProtTomoBase):
     """ This protocol subtracts a subtomogram average to a SetOfSubtomograms, which are internally aligned and
@@ -46,7 +49,7 @@ class XmippProtSubtractionSubtomo(EMProtocol, ProtTomoBase):
     determined region."""
 
     _label = 'subtomo subtraction'
-    _devStatus = BETA
+    _devStatus = UPDATED
 
     # --------------------------- DEFINE param functions --------------------------------------------
     def _defineParams(self, form):
@@ -77,14 +80,18 @@ def _defineParams(self, form):
                       expertLevel=LEVEL_ADVANCED, help='Save input volume 1 (first subtomogram of the set) filtered '
                                                        'and input volume 2 (average) adjusted, which are the volumes '
                                                        'that are really subtracted.')
+        form.addParallelSection(threads=0, mpi=4)
 
     # --------------------------- INSERT steps functions --------------------------------------------
     def _insertAllSteps(self):
-        self._insertFunctionStep('applyAlignStep')
+        self._insertFunctionStep('convertStep')
         self._insertFunctionStep('subtractionStep')
         self._insertFunctionStep('createOutputStep')
 
     # --------------------------- STEPS functions --------------------------------------------
+    def convertStep(self):
+        writeSetOfParticles(self.inputSubtomos.get(), self._getExtraPath(self.INPUT_PARTICLES))
+
     def applyAlignStep(self):
         """Align subtomograms to be in the same orientation as the reference"""
         inputSet = self.inputSubtomos.get()
@@ -156,7 +163,7 @@ def subtractionStep(self):
         resol = self.resol.get()
         iter = self.iter.get()
         program = "xmipp_volume_subtraction"
-        args = '--i2 %s --iter %s --lambda %s --sub' % \
+        args = '--ref %s --iter %s --lambda %s --sub' % \
                (average, iter, self.rfactor.get())
         if resol:
             fc = self.inputSubtomos.get().getSamplingRate() / resol
@@ -191,7 +198,7 @@ def createOutputStep(self):
         inputSubtomos = self.inputSubtomos.get()
         outputSet = self._createSetOfSubTomograms()
         outputSet.copyInfo(inputSubtomos)
-        outputSet.copyItems(inputSubtomos, updateItemCallback=self._updateItemOutput)
+        readSetOfParticles(self._getExtraPath(OUTPUT), outputSet)
         self._defineOutputs(outputSubtomograms=outputSet)
         self._defineSourceRelation(inputSubtomos, outputSet)
 

From 0464958447c7c17d52f04351b911bbd7a1b1a82c Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Fri, 15 Sep 2023 12:20:23 +0200
Subject: [PATCH 03/11] remove padding and apply geo (now in program)

---
 .../protocols/protocol_subtraction_subtomo.py | 110 ++----------------
 1 file changed, 7 insertions(+), 103 deletions(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 7a0c05a1..095531b1 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -29,18 +29,12 @@
 from pyworkflow import BETA, UPDATED, NEW, PROD
 from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
-
-from pwem import ALIGN_3D
-from pwem.emlib import MDL_IMAGE, MDL_ITEM_ID
-import pwem.emlib.metadata as md
 from pwem.protocols import EMProtocol
-
 from tomo.protocols import ProtTomoBase
 from xmipp3.convert import writeSetOfParticles, readSetOfParticles
-from xmipp3.convert import xmippToLocation, writeSetOfVolumes, alignmentToRow
 
 OUTPUT = "output_subtomograms.xmd"
-INPUT_PARTICLES = "input_subtomograms.xmd"
+INPUT = "input_subtomograms.xmd"
 
 class XmippProtSubtractionSubtomo(EMProtocol, ProtTomoBase):
     """ This protocol subtracts a subtomogram average to a SetOfSubtomograms, which are internally aligned and
@@ -90,70 +84,7 @@ def _insertAllSteps(self):
 
     # --------------------------- STEPS functions --------------------------------------------
     def convertStep(self):
-        writeSetOfParticles(self.inputSubtomos.get(), self._getExtraPath(self.INPUT_PARTICLES))
-
-    def applyAlignStep(self):
-        """Align subtomograms to be in the same orientation as the reference"""
-        inputSet = self.inputSubtomos.get()
-        outputFn = self._getExtraPath('input_subtomos.xmd')
-        if inputSet.getFirstItem().getFileName().endswith('.mrc') or \
-                inputSet.getFirstItem().getFileName().endswith('.map'):
-            S = self._createSetOfSubTomograms()
-            S.setSamplingRate(inputSet.getSamplingRate())
-            for subtomo in inputSet:
-                s = subtomo.clone()
-                s.setFileName(subtomo.getFileName() + ':mrc')
-                S.append(s)
-            writeSetOfVolumes(S, outputFn, alignType=ALIGN_3D)
-        else:
-            writeSetOfVolumes(inputSet, outputFn, alignType=ALIGN_3D)
-
-        # Window subtomograms twice their size
-        windowedStk = self._getExtraPath('windowed_subtomograms.stk')
-        self.runJob('xmipp_transform_window', '-i %s -o %s --size %d --save_metadata_stack' %
-                    (outputFn, windowedStk, 2 * inputSet.getFirstItem().getDim()[0]), numberOfMpi=1)
-
-        # Add input transform matrix to md generated by xmipp_transform_window
-        mdWindow = md.MetaData(self._getExtraPath('windowed_subtomograms.xmd'))
-        mdWindowTransform = md.MetaData()
-        mdInput = md.MetaData(self._getExtraPath('input_subtomos.xmd'))
-        idList = list(inputSet.getIdSet())
-        for rowWin,rowInp in zip(md.iterRows(mdWindow), md.iterRows(mdInput)):
-            rowOut = md.Row()
-            rowOut.copyFromRow(rowInp)
-            idx = rowWin.getValue(MDL_IMAGE)
-            idx = idx.split('@')[0]
-            idx = idx.strip('0')
-            alignmentToRow(inputSet[(idList[int(idx) - 1])].getTransform(), rowOut, ALIGN_3D)
-            rowOut.addToMd(mdWindowTransform)
-        mdWindowTransform.write(self._getExtraPath("window_with_original_geometry.xmd"))
-
-        # Align subtomograms
-        self.runJob('xmipp_transform_geometry', '-i %s -o %s --apply_transform --dont_wrap' %
-                    (self._getExtraPath("window_with_original_geometry.xmd"),
-                     self._getExtraPath('aligned_subtomograms.stk')))
-
-        # Window subtomograms to their original size
-        alignMd = self._getExtraPath('aligned_subtomograms.xmd')
-        outputStk = self._getExtraPath('output_subtomograms.stk')
-        self.runJob('xmipp_transform_window', '-i %s -o %s --size %d ' %
-                    (alignMd, outputStk, inputSet.getFirstItem().getDim()[0]),
-                    numberOfMpi=1)
-                      
-        #outputMd = md.MetaData(self._getExtraPath('output_subtomograms.stk'))
-        outputMd = md.MetaData()
-        mdAli = md.MetaData(self._getExtraPath('aligned_subtomograms.xmd'))
-        for rowAli in md.iterRows(mdAli):
-            rowOutp = md.Row()
-            rowOutp.copyFromRow(rowAli)
-            rowOutp.addToMd(outputMd)
-        outputMd.write(self._getExtraPath("output_subtomograms.xmd"))
-        self.alignedSet = self._createSetOfSubTomograms()
-        self.alignedSet.copyInfo(inputSet)
-        inputMd = self._getExtraPath('output_subtomograms.xmd')
-        self.alignedSet.copyItems(inputSet,
-                             updateItemCallback=self._updateItemAlign,
-                             itemDataIterator=md.iterRows(inputMd, sortByLabel=MDL_ITEM_ID))             
+        writeSetOfParticles(self.inputSubtomos.get(), self._getExtraPath(INPUT))
 
     def subtractionStep(self):
         """Subtract reference to each of the subtomogram in the input Set"""
@@ -163,8 +94,8 @@ def subtractionStep(self):
         resol = self.resol.get()
         iter = self.iter.get()
         program = "xmipp_volume_subtraction"
-        args = '--ref %s --iter %s --lambda %s --sub' % \
-               (average, iter, self.rfactor.get())
+        args = '-i %s -o % s --ref %s --iter %s --lambda %s --sub --subtomos' % \
+               (self._getExtraPath(INPUT), self._getExtraPath(OUTPUT), average, iter, self.rfactor.get())
         if resol:
             fc = self.inputSubtomos.get().getSamplingRate() / resol
             args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())
@@ -173,27 +104,10 @@ def subtractionStep(self):
             maskSub = self.maskSub.get()
             if maskSub:
                 args += ' --maskSub %s' % maskSub.getFileName()
+        if self.saveFiles:
+            args += ' --saveV1 %s --saveV2 %s' % (self._getExtraPath('v1.mrc'), self._getExtraPath('v2.mrc'))
+        self.runJob(program, args)
 
-        mdOrig = md.MetaData(self._getExtraPath('window_with_original_geometry.xmd'))
-        mdOut = md.MetaData(self._getExtraPath('output_subtomograms.xmd'))
-        for subtomo, row, rowOut in zip(self.alignedSet.iterItems(), md.iterRows(mdOrig), md.iterRows(mdOut)):
-                if subtomo.getObjId() == row.getValue(MDL_ITEM_ID):
-                        ix = subtomo.getIndex()
-                        fnOutSubtomo = self._getExtraPath("output_subtomo%06d.mrc" % row.getValue(MDL_ITEM_ID))
-                        argsSubtomo = ' --i1 %s -o % s' % (rowOut.getValue(MDL_IMAGE), fnOutSubtomo)
-                        if self.saveFiles and ix == 1:
-                                argsSubtomo += ' --saveV1 %s --saveV2 %s' % (self._getExtraPath('vol1_filtered.mrc'), 
-                                self._getExtraPath('average_adjusted.mrc'))
-                        print('\n-----Subtomogram %d-----' % ix)
-                        self.runJob(program, args + argsSubtomo)
-                        # Apply inverse transform for the output to have the original orientation
-                        self.runJob('xmipp_transform_geometry', 
-                        '-i %s -o %s --rotate_volume euler %d %d %d --shift %d %d %d --inverse --dont_wrap' %
-                        (fnOutSubtomo, fnOutSubtomo, row.getValue('angleRot'), row.getValue('angleTilt'), 
-                        row.getValue('anglePsi'), row.getValue('shiftX'), row.getValue('shiftY'), row.getValue('shiftZ')))        
-                else:
-                        continue            
-        
     def createOutputStep(self):
         inputSubtomos = self.inputSubtomos.get()
         outputSet = self._createSetOfSubTomograms()
@@ -237,13 +151,3 @@ def _validate(self):
                 validateMsgs.append(
                     'Please provide subtomograms which have transformation matrix.')
         return validateMsgs
-
-    # --------------------------- UTLIS functions --------------------------------------------
-    def _updateItemAlign(self, item, row):
-        newFn = row.getValue(md.MDL_IMAGE)
-        newLoc = xmippToLocation(newFn)
-        item.setLocation(newLoc)
-        item.setTransform(None)
-
-    def _updateItemOutput(self, item, row):
-        item.setFileName(self._getExtraPath("output_subtomo%06d.mrc" % item.getObjId()))

From 015650fd5979c31d81bac5df4e5f89cd63f13cd5 Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Thu, 5 Oct 2023 13:17:00 +0200
Subject: [PATCH 04/11] implement writeSetOfSubtomos

---
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 095531b1..3889736a 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -30,8 +30,10 @@
 from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 from pwem.protocols import EMProtocol
+from pwem import ALIGN_3D
 from tomo.protocols import ProtTomoBase
-from xmipp3.convert import writeSetOfParticles, readSetOfParticles
+from xmipptomo.convert import writeSetOfSubtomograms
+from xmipp3.convert import readSetOfParticles
 
 OUTPUT = "output_subtomograms.xmd"
 INPUT = "input_subtomograms.xmd"
@@ -84,7 +86,7 @@ def _insertAllSteps(self):
 
     # --------------------------- STEPS functions --------------------------------------------
     def convertStep(self):
-        writeSetOfParticles(self.inputSubtomos.get(), self._getExtraPath(INPUT))
+        writeSetOfSubtomograms(self.inputSubtomos.get(), self._getExtraPath(INPUT), aligntype=ALIGN_3D)
 
     def subtractionStep(self):
         """Subtract reference to each of the subtomogram in the input Set"""

From e8ccd6a758bffced9e70b9a3a3529f80f50e4ac8 Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Thu, 5 Oct 2023 13:17:25 +0200
Subject: [PATCH 05/11] implement writeSetOfSubtomograms

---
 xmipptomo/convert/__init__.py |  27 ++++
 xmipptomo/convert/convert.py  | 271 ++++++++++++++++++++++++++++++++++
 2 files changed, 298 insertions(+)
 create mode 100644 xmipptomo/convert/__init__.py
 create mode 100644 xmipptomo/convert/convert.py

diff --git a/xmipptomo/convert/__init__.py b/xmipptomo/convert/__init__.py
new file mode 100644
index 00000000..c8371661
--- /dev/null
+++ b/xmipptomo/convert/__init__.py
@@ -0,0 +1,27 @@
+# **************************************************************************
+# *
+# * Authors:     Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *
+# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
+# *
+# * This program is free software; you can redistribute it and/or modify
+# * it under the terms of the GNU General Public License as published by
+# * the Free Software Foundation; either version 2 of the License, or
+# * (at your option) any later version.
+# *
+# * This program is distributed in the hope that it will be useful,
+# * but WITHOUT ANY WARRANTY; without even the implied warranty of
+# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# * GNU General Public License for more details.
+# *
+# * You should have received a copy of the GNU General Public License
+# * along with this program; if not, write to the Free Software
+# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+# * 02111-1307  USA
+# *
+# *  All comments concerning this program package may be sent to the
+# *  e-mail address 'scipion@cnb.csic.es'
+# *
+# **************************************************************************
+
+from .convert import *
diff --git a/xmipptomo/convert/convert.py b/xmipptomo/convert/convert.py
new file mode 100644
index 00000000..face1e7b
--- /dev/null
+++ b/xmipptomo/convert/convert.py
@@ -0,0 +1,271 @@
+# **************************************************************************
+# *
+# * Authors:     Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
+# *
+# * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
+# *
+# * This program is free software; you can redistribute it and/or modify
+# * it under the terms of the GNU General Public License as published by
+# * the Free Software Foundation; either version 2 of the License, or
+# * (at your option) any later version.
+# *
+# * This program is distributed in the hope that it will be useful,
+# * but WITHOUT ANY WARRANTY; without even the implied warranty of
+# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# * GNU General Public License for more details.
+# *
+# * You should have received a copy of the GNU General Public License
+# * along with this program; if not, write to the Free Software
+# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+# * 02111-1307  USA
+# *
+# *  All comments concerning this program package may be sent to the
+# *  e-mail address 'scipion@cnb.csic.es'
+# *
+# **************************************************************************
+from pwem import emlib
+import pwem.emlib.metadata as md
+from pwem.emlib.image import ImageHandler
+from pwem.constants import (NO_INDEX, ALIGN_NONE, ALIGN_PROJ, ALIGN_2D, ALIGN_3D)
+import numpy as np
+from collections import OrderedDict
+from xmipp3.base import getLabelPythonType, iterMdRows
+
+def prefixAttribute(attribute):
+    return '_xmipp_%s' % attribute
+
+if not getattr(emlib, "GHOST_ACTIVATED", False):
+    """ Some of MDL may not exist when Ghost is activated
+    """
+    # This dictionary will be used to map
+    # between CTFModel properties and Xmipp labels
+    ACQUISITION_DICT = OrderedDict([
+           ("_amplitudeContrast", emlib.MDL_CTF_Q0),
+           ("_sphericalAberration", emlib.MDL_CTF_CS),
+           ("_voltage", emlib.MDL_CTF_VOLTAGE)
+           ])
+
+    COOR_DICT = OrderedDict([
+        ("_x", emlib.MDL_XCOOR),
+        ("_y", emlib.MDL_YCOOR),
+        ("_z", emlib.MDL_ZCOOR)
+    ])
+
+    ALIGNMENT_DICT = OrderedDict([
+           (prefixAttribute("shiftX"), emlib.MDL_SHIFT_X),
+           (prefixAttribute("shiftY"), emlib.MDL_SHIFT_Y),
+           (prefixAttribute("shiftZ"), emlib.MDL_SHIFT_Z),
+           (prefixAttribute("flip"), emlib.MDL_FLIP),
+           (prefixAttribute("anglePsi"), emlib.MDL_ANGLE_PSI),
+           (prefixAttribute("angleRot"), emlib.MDL_ANGLE_ROT),
+           (prefixAttribute("angleTilt"), emlib.MDL_ANGLE_TILT),
+           ])
+
+def writeSetOfSubtomograms(imgSet, filename, blockName='Particles', **kwargs):
+    writeSetOfImages(imgSet, filename, particleToRow, blockName, **kwargs)
+
+def writeSetOfImages(imgSet, filename, imgToFunc,
+                     blockName='Images', **kwargs):
+    """ This function will write a SetOfImages as a Xmipp metadata.
+    Params:
+        imgSet: the set of images to be written (particles,
+        micrographs or volumes)
+        filename: the filename where to write the metadata.
+        rowFunc: this function can be used to setup the row before
+            adding to metadata.
+    """
+    mdSet = emlib.MetaData()
+
+    setOfImagesToMd(imgSet, mdSet, imgToFunc, **kwargs)
+    mdSet.write('%s@%s' % (blockName, filename))
+
+def particleToRow(part, partRow, **kwargs):
+    """ Set labels values from Particle to md row. """
+    imageToRow(part, partRow, emlib.MDL_IMAGE, **kwargs)
+    coord = part.getCoordinate3D()
+    if coord is not None:
+        coordinateToRow(coord, partRow, copyId=False)
+
+def setOfImagesToMd(imgSet, mdIn, imgToFunc, **kwargs):
+    """ This function will fill Xmipp metadata from a SetOfMicrographs
+    Params:
+        imgSet: the set of images to be converted to metadata
+        mdIn: metadata to be filled
+        rowFunc: this function can be used to setup the row before
+            adding to metadata.
+    """
+
+    if 'alignType' not in kwargs:
+        kwargs['alignType'] = imgSet.getAlignment()
+
+    if 'where' in kwargs:
+        where = kwargs['where']
+        for img in imgSet.iterItems(where=where):
+            objId = mdIn.addObject()
+            imgRow = md.Row()
+            imgToFunc(img, imgRow, **kwargs)
+            imgRow.writeToMd(mdIn, objId)
+    else:
+        for img in imgSet:
+            objId = mdIn.addObject()
+            imgRow = md.Row()
+            imgToFunc(img, imgRow, **kwargs)
+            imgRow.writeToMd(mdIn, objId)
+
+def imageToRow(img, imgRow, imgLabel, **kwargs):
+    # Provide a hook to be used if something is needed to be
+    # done for special cases before converting image to row
+    preprocessImageRow = kwargs.get('preprocessImageRow', None)
+    if preprocessImageRow:
+        preprocessImageRow(img, imgRow)
+
+    setRowId(imgRow, img)  # Set the id in the metadata as MDL_ITEM_ID
+    index, filename = img.getLocation()
+    fn = locationToXmipp(index, filename)
+    imgRow.setValue(imgLabel, fn)
+
+    # alignment is mandatory at this point, it shoud be check
+    # and detected defaults if not passed at readSetOf.. level
+    alignType = kwargs.get('alignType')
+
+    if alignType != ALIGN_NONE:
+        alignmentToRow(img.getTransform(), imgRow, alignType)
+
+    if kwargs.get('writeAcquisition', True) and img.hasAcquisition():
+        acquisitionToRow(img.getAcquisition(), imgRow)
+
+    # Write all extra labels to the row
+    objectToRow(img, imgRow, {})
+
+    # Provide a hook to be used if something is needed to be
+    # done for special cases before converting image to row
+    postprocessImageRow = kwargs.get('postprocessImageRow', None)
+    if postprocessImageRow:
+        postprocessImageRow(img, imgRow)
+
+def coordinateToRow(coord, coordRow, copyId=True):
+    """ Set labels values from Coordinate coord to md row. """
+    if copyId:
+        setRowId(coordRow, coord)
+    objectToRow(coord, coordRow, COOR_DICT)
+    if coord.getTomoId():
+        coordRow.setValue(emlib.MDL_MICROGRAPH, str(coord.getTomoId()))
+
+def setRowId(mdRow, obj, label=emlib.MDL_ITEM_ID):
+    mdRow.setValue(label, int(obj.getObjId()))
+
+
+def objectToRow(obj, row, attrDict, extraLabels=[]):
+    """ This function will convert an EMObject into a Row.
+    Params:
+        obj: the EMObject instance (input)
+        row: the Row instance (output)  -see emlib.metadata.utils.Row()-
+        attrDict: dictionary with the map between obj attributes(keys) and
+            row MDLabels in Xmipp (values).
+        extraLabels: a list with extra labels that could be included
+            as _xmipp_labelName
+    """
+    if obj.isEnabled():
+        enabled = 1
+    else:
+        enabled = -1
+    row.setValue(emlib.MDL_ENABLED, enabled)
+
+    for attr, label in attrDict.items():
+        if hasattr(obj, attr):
+            valueType = getLabelPythonType(label)
+            value = getattr(obj, attr).get()
+            try:
+                row.setValue(label, valueType(value))
+            except Exception as e:
+                print(e)
+                print("Problems found converting metadata: ")
+                print("Label id = %s" % label)
+                print("Attribute = %s" % attr)
+                print("Value = %s" % value)
+                print("Value type = %s" % valueType)
+                raise e
+            row.setValue(label, valueType(getattr(obj, attr).get()))
+
+    attrLabels = attrDict.values()
+
+    for label in extraLabels:
+        attrName = prefixAttribute(emlib.label2Str(label))
+        if label not in attrLabels and hasattr(obj, attrName):
+            value = obj.getAttributeValue(attrName)
+            row.setValue(label, value)
+
+def locationToXmipp(index, filename):
+    """ Convert an index and filename location
+    to a string with @ as expected in Xmipp.
+    """
+    return ImageHandler.locationToXmipp((index, filename))
+
+def alignmentToRow(alignment, alignmentRow, alignType):
+    """
+    is2D == True-> matrix is 2D (2D images alignment)
+            otherwise matrix is 3D (3D volume alignment or projection)
+    invTransform == True  -> for xmipp implies projection
+                          -> for xmipp implies alignment
+    """
+    if alignment is None:
+        return
+
+    is2D = alignType == ALIGN_2D
+    inverseTransform = alignType == ALIGN_PROJ
+    # only flip is meaninfull if 2D case
+    # in that case the 2x2 determinant is negative
+    flip = False
+    matrix = alignment.getMatrix()
+    if alignType == ALIGN_2D:
+        # get 2x2 matrix and check if negative
+        flip = bool(np.linalg.det(matrix[0:2, 0:2]) < 0)
+        if flip:
+            matrix[0, :2] *= -1.  # invert only the first two columns keep x
+            matrix[2, 2] = 1.  # set 3D rot
+        else:
+            pass
+
+    elif alignType == ALIGN_3D:
+        flip = bool(np.linalg.det(matrix[0:3, 0:3]) < 0)
+        if flip:
+            matrix[0, :4] *= -1.  # now, invert first line including x
+            matrix[3, 3] = 1.  # set 3D rot
+        else:
+            pass
+
+    else:
+        flip = bool(np.linalg.det(matrix[0:3, 0:3]) < 0)
+        if flip:
+            raise Exception("the det of the transformation matrix is "
+                            "negative. This is not a valid transformation "
+                            "matrix for Scipion.")
+    shifts, angles = geometryFromMatrix(matrix, inverseTransform)
+    alignmentRow.setValue(emlib.MDL_SHIFT_X, shifts[0])
+    alignmentRow.setValue(emlib.MDL_SHIFT_Y, shifts[1])
+
+    if is2D:
+        angle = angles[0] + angles[2]
+        alignmentRow.setValue(emlib.MDL_ANGLE_PSI,  angle)
+    else:
+        # if alignType == ALIGN_3D:
+        alignmentRow.setValue(emlib.MDL_SHIFT_Z, shifts[2])
+        alignmentRow.setValue(emlib.MDL_ANGLE_ROT,  angles[0])
+        alignmentRow.setValue(emlib.MDL_ANGLE_TILT, angles[1])
+        alignmentRow.setValue(emlib.MDL_ANGLE_PSI,  angles[2])
+    alignmentRow.setValue(emlib.MDL_FLIP, flip)
+
+def acquisitionToRow(acquisition, ctfRow):
+    """ Set labels values from acquisition to md row. """
+    objectToRow(acquisition, ctfRow, ACQUISITION_DICT)
+
+def geometryFromMatrix(matrix, inverseTransform):
+    from pwem.convert.transformations import (translation_from_matrix,
+                                              euler_from_matrix)
+    if inverseTransform:
+        matrix = np.linalg.inv(matrix)
+        shifts = -translation_from_matrix(matrix)
+    else:
+        shifts = translation_from_matrix(matrix)
+    angles = -np.rad2deg(euler_from_matrix(matrix, axes='szyz'))
+    return shifts, angles
\ No newline at end of file

From 1d2138b138fc5e0500749507f66642de3302363d Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Mon, 9 Oct 2023 11:58:47 +0200
Subject: [PATCH 06/11] convert alignment fix

---
 xmipptomo/convert/convert.py                        | 5 +++--
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 4 ++--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/xmipptomo/convert/convert.py b/xmipptomo/convert/convert.py
index face1e7b..8309515f 100644
--- a/xmipptomo/convert/convert.py
+++ b/xmipptomo/convert/convert.py
@@ -26,7 +26,8 @@
 from pwem import emlib
 import pwem.emlib.metadata as md
 from pwem.emlib.image import ImageHandler
-from pwem.constants import (NO_INDEX, ALIGN_NONE, ALIGN_PROJ, ALIGN_2D, ALIGN_3D)
+from pwem.constants import (NO_INDEX, ALIGN_NONE, ALIGN_PROJ, ALIGN_2D,
+                            ALIGN_3D)
 import numpy as np
 from collections import OrderedDict
 from xmipp3.base import getLabelPythonType, iterMdRows
@@ -74,8 +75,8 @@ def writeSetOfImages(imgSet, filename, imgToFunc,
         rowFunc: this function can be used to setup the row before
             adding to metadata.
     """
-    mdSet = emlib.MetaData()
 
+    mdSet = emlib.MetaData()
     setOfImagesToMd(imgSet, mdSet, imgToFunc, **kwargs)
     mdSet.write('%s@%s' % (blockName, filename))
 
diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 3889736a..063fe5e6 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -30,7 +30,7 @@
 from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam, FloatParam
 from pyworkflow.protocol.constants import LEVEL_ADVANCED
 from pwem.protocols import EMProtocol
-from pwem import ALIGN_3D
+from pwem.constants import ALIGN_3D
 from tomo.protocols import ProtTomoBase
 from xmipptomo.convert import writeSetOfSubtomograms
 from xmipp3.convert import readSetOfParticles
@@ -86,7 +86,7 @@ def _insertAllSteps(self):
 
     # --------------------------- STEPS functions --------------------------------------------
     def convertStep(self):
-        writeSetOfSubtomograms(self.inputSubtomos.get(), self._getExtraPath(INPUT), aligntype=ALIGN_3D)
+        writeSetOfSubtomograms(self.inputSubtomos.get(), self._getExtraPath(INPUT), alignType=ALIGN_3D)
 
     def subtractionStep(self):
         """Subtract reference to each of the subtomogram in the input Set"""

From 04709668272e45b90814944b684fb5ea45f8a34c Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Wed, 18 Oct 2023 10:25:39 +0200
Subject: [PATCH 07/11] add --oroot

---
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 063fe5e6..8a7044a0 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -96,8 +96,9 @@ def subtractionStep(self):
         resol = self.resol.get()
         iter = self.iter.get()
         program = "xmipp_volume_subtraction"
-        args = '-i %s -o % s --ref %s --iter %s --lambda %s --sub --subtomos' % \
-               (self._getExtraPath(INPUT), self._getExtraPath(OUTPUT), average, iter, self.rfactor.get())
+        args = '-i %s -o % s --ref %s --oroot %s --iter %s --lambda %s --sub --subtomos' % \
+               (self._getExtraPath(INPUT), self._getExtraPath(OUTPUT), average,
+                self._getExtraPath("subtracted_subtomo"), iter, self.rfactor.get())
         if resol:
             fc = self.inputSubtomos.get().getSamplingRate() / resol
             args += ' --cutFreq %f --sigma %d' % (fc, self.sigma.get())

From 0bf37d316defde8fb3bd853846f4e288c181a1f2 Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Wed, 18 Oct 2023 12:18:01 +0200
Subject: [PATCH 08/11] create read part of convert + use it in subtraction

---
 xmipptomo/convert/convert.py                  | 398 +++++++++++++++++-
 .../protocols/protocol_subtraction_subtomo.py |   5 +-
 2 files changed, 399 insertions(+), 4 deletions(-)

diff --git a/xmipptomo/convert/convert.py b/xmipptomo/convert/convert.py
index 8309515f..8bf5d04c 100644
--- a/xmipptomo/convert/convert.py
+++ b/xmipptomo/convert/convert.py
@@ -28,9 +28,12 @@
 from pwem.emlib.image import ImageHandler
 from pwem.constants import (NO_INDEX, ALIGN_NONE, ALIGN_PROJ, ALIGN_2D,
                             ALIGN_3D)
+from pwem.objects import Acquisition, Transform
+from pyworkflow.object import ObjectWrap, String
 import numpy as np
 from collections import OrderedDict
 from xmipp3.base import getLabelPythonType, iterMdRows
+from tomo.objects import SubTomogram, Coordinate3D, CTFTomo
 
 def prefixAttribute(attribute):
     return '_xmipp_%s' % attribute
@@ -62,6 +65,73 @@ def prefixAttribute(attribute):
            (prefixAttribute("angleTilt"), emlib.MDL_ANGLE_TILT),
            ])
 
+    CTF_DICT = OrderedDict([
+           ("_defocusU", emlib.MDL_CTF_DEFOCUSU),
+           ("_defocusV", emlib.MDL_CTF_DEFOCUSV),
+           ("_defocusAngle", emlib.MDL_CTF_DEFOCUS_ANGLE),
+           ("_resolution", emlib.MDL_CTF_CRIT_MAXFREQ),
+           ("_fitQuality", emlib.MDL_CTF_CRIT_FITTINGSCORE)
+           ])
+
+    # TODO: remove next dictionary when all
+    # cTFmodel has resolution and fitQuality
+    CTF_DICT_NORESOLUTION = OrderedDict([
+            ("_defocusU", emlib.MDL_CTF_DEFOCUSU),
+            ("_defocusV", emlib.MDL_CTF_DEFOCUSV),
+            ("_defocusAngle", emlib.MDL_CTF_DEFOCUS_ANGLE)
+            ])
+
+    CTF_PSD_DICT = OrderedDict([
+           ("_psdFile", emlib.MDL_PSD),
+           (prefixAttribute("enhanced_psd"), emlib.MDL_PSD_ENHANCED),
+           (prefixAttribute("ctfmodel_quadrant"), emlib.MDL_IMAGE1),
+           (prefixAttribute("ctfmodel_halfplane"), emlib.MDL_IMAGE1)
+           ])
+
+    CTF_EXTRA_LABELS = [
+        emlib.MDL_CTF_CA,
+        emlib.MDL_CTF_ENERGY_LOSS,
+        emlib.MDL_CTF_LENS_STABILITY,
+        emlib.MDL_CTF_CONVERGENCE_CONE,
+        emlib.MDL_CTF_LONGITUDINAL_DISPLACEMENT,
+        emlib.MDL_CTF_TRANSVERSAL_DISPLACEMENT,
+        emlib.MDL_CTF_K,
+        emlib.MDL_CTF_BG_GAUSSIAN_K,
+        emlib.MDL_CTF_BG_GAUSSIAN_SIGMAU,
+        emlib.MDL_CTF_BG_GAUSSIAN_SIGMAV,
+        emlib.MDL_CTF_BG_GAUSSIAN_CU,
+        emlib.MDL_CTF_BG_GAUSSIAN_CV,
+        emlib.MDL_CTF_BG_SQRT_K,
+        emlib.MDL_CTF_BG_SQRT_U,
+        emlib.MDL_CTF_BG_SQRT_V,
+        emlib.MDL_CTF_BG_SQRT_ANGLE,
+        emlib.MDL_CTF_BG_BASELINE,
+        emlib.MDL_CTF_BG_GAUSSIAN2_K,
+        emlib.MDL_CTF_BG_GAUSSIAN2_SIGMAU,
+        emlib.MDL_CTF_BG_GAUSSIAN2_SIGMAV,
+        emlib.MDL_CTF_BG_GAUSSIAN2_CU,
+        emlib.MDL_CTF_BG_GAUSSIAN2_CV,
+        emlib.MDL_CTF_BG_GAUSSIAN2_ANGLE,
+        emlib.MDL_CTF_CRIT_FITTINGCORR13,
+        emlib.MDL_CTF_CRIT_ICENESS,
+        emlib.MDL_CTF_VPP_RADIUS,
+        emlib.MDL_CTF_DOWNSAMPLE_PERFORMED,
+        emlib.MDL_CTF_CRIT_PSDVARIANCE,
+        emlib.MDL_CTF_CRIT_PSDPCA1VARIANCE,
+        emlib.MDL_CTF_CRIT_PSDPCARUNSTEST,
+        emlib.MDL_CTF_CRIT_FIRSTZEROAVG,
+        emlib.MDL_CTF_CRIT_DAMPING,
+        emlib.MDL_CTF_CRIT_FIRSTZERORATIO,
+        emlib.MDL_CTF_CRIT_PSDCORRELATION90,
+        emlib.MDL_CTF_CRIT_PSDRADIALINTEGRAL,
+        emlib.MDL_CTF_CRIT_NORMALITY,
+        # In xmipp the ctf also contains acquisition information
+        emlib.MDL_CTF_Q0,
+        emlib.MDL_CTF_CS,
+        emlib.MDL_CTF_VOLTAGE,
+        emlib.MDL_CTF_SAMPLING_RATE
+        ]
+
 def writeSetOfSubtomograms(imgSet, filename, blockName='Particles', **kwargs):
     writeSetOfImages(imgSet, filename, particleToRow, blockName, **kwargs)
 
@@ -244,6 +314,7 @@ def alignmentToRow(alignment, alignmentRow, alignType):
     shifts, angles = geometryFromMatrix(matrix, inverseTransform)
     alignmentRow.setValue(emlib.MDL_SHIFT_X, shifts[0])
     alignmentRow.setValue(emlib.MDL_SHIFT_Y, shifts[1])
+    alignmentRow.setValue(emlib.MDL_SHIFT_Z, shifts[2])
 
     if is2D:
         angle = angles[0] + angles[2]
@@ -269,4 +340,329 @@ def geometryFromMatrix(matrix, inverseTransform):
     else:
         shifts = translation_from_matrix(matrix)
     angles = -np.rad2deg(euler_from_matrix(matrix, axes='szyz'))
-    return shifts, angles
\ No newline at end of file
+    return shifts, angles
+
+def readSetOfSubtomograms(filename, partSet, **kwargs):
+    readSetOfImages(filename, partSet, rowToParticle, **kwargs)
+
+def readSetOfImages(filename, imgSet, rowToFunc, **kwargs):
+    """read from Xmipp image metadata.
+        filename: The metadata filename where the image are.
+        imgSet: the SetOfParticles that will be populated.
+        rowToFunc: this function will be used to convert the row to Object
+    """
+    imgMd = emlib.MetaData(filename)
+
+    # By default remove disabled items from metadata
+    # be careful if you need to preserve the original number of items
+    if kwargs.get('removeDisabled', True):
+        imgMd.removeDisabled()
+
+    # If the type of alignment is not sent through the kwargs
+    # try to deduced from the metadata labels
+    if 'alignType' not in kwargs:
+        imgRow = rowFromMd(imgMd, imgMd.firstObject())
+        if _containsAny(imgRow, ALIGNMENT_DICT):
+            if imgRow.containsLabel(emlib.MDL_ANGLE_TILT):
+                kwargs['alignType'] = ALIGN_PROJ
+            else:
+                kwargs['alignType'] = ALIGN_2D
+        else:
+            kwargs['alignType'] = ALIGN_NONE
+
+    if imgMd.size() > 0:
+        for objId in imgMd:
+            imgRow = rowFromMd(imgMd, objId)
+            img = rowToFunc(imgRow, **kwargs)
+            imgSet.append(img)
+
+        imgSet.setHasCTF(img.hasCTF())
+        imgSet.setAlignment(kwargs['alignType'])
+
+def rowToParticle(partRow, **kwargs):
+    return _rowToParticle(partRow, SubTomogram, **kwargs)
+
+def _rowToParticle(partRow, particleClass, **kwargs):
+    """ Create a Particle from a row of a metadata. """
+    # Since postprocessImage is intended to be after the object is
+    # setup, we need to intercept it here and call it at the end
+    postprocessImageRow = kwargs.get('postprocessImageRow', None)
+    if postprocessImageRow:
+        del kwargs['postprocessImageRow']
+
+    img = rowToImage(partRow, emlib.MDL_IMAGE, particleClass, **kwargs)
+    img.setCoordinate3D(rowToCoordinate(partRow))
+    # copy micId if available
+    # if not copy micrograph name if available
+    try:
+        if partRow.hasLabel(emlib.MDL_MICROGRAPH_ID):
+            img.setMicId(partRow.getValue(emlib.MDL_MICROGRAPH_ID))
+#        elif partRow.hasLabel(emlib.MDL_MICROGRAPH):
+#            micName = partRow.getValue(emlib.MDL_MICROGRAPH)
+#            img._micrograph = micName
+#            print("setting micname as %s" % micName)
+#            img.printAll()
+#            print("getAttributes1 %s" % img._micrograph)
+#            print("getAttributes2 %s" % getattr(img, "_micrograph", 'kk')
+#        else:
+#            print("WARNING: No micname")
+    except Exception as e:
+        print("Warning:", e.message)
+
+    if postprocessImageRow:
+        postprocessImageRow(img, partRow)
+
+    return img
+
+def rowFromMd(mdIn, objId):
+    row = md.Row()
+    row.readFromMd(mdIn, objId)
+    return row
+
+def _containsAny(row, labels):
+    """ Check if the labels (values) in labelsDict
+    are present in the row.
+    """
+    values = labels.values() if isinstance(labels, dict) else labels
+    return any(row.containsLabel(l) for l in values)
+
+def rowToImage(imgRow, imgLabel, imgClass, **kwargs):
+    """ Create an Image from a row of a metadata. """
+    img = imgClass()
+
+    # Provide a hook to be used if something is needed to be
+    # done for special cases before converting image to row
+    preprocessImageRow = kwargs.get('preprocessImageRow', None)
+    if preprocessImageRow:
+        preprocessImageRow(img, imgRow)
+
+    # Decompose Xmipp filename
+    index, filename = xmippToLocation(imgRow.getValue(imgLabel))
+    img.setLocation(index, filename)
+
+    if imgRow.containsLabel(emlib.MDL_REF):
+        img.setClassId(imgRow.getValue(emlib.MDL_REF))
+    elif imgRow.containsLabel(emlib.MDL_REF3D):
+        img.setClassId(imgRow.getValue(emlib.MDL_REF3D))
+
+    if kwargs.get('readCtf', True):
+        img.setCTF(rowToCtfModel(imgRow))
+
+    # alignment is mandatory at this point, it shoud be check
+    # and detected defaults if not passed at readSetOf.. level
+    alignType = kwargs.get('alignType')
+
+    if alignType != ALIGN_NONE:
+        img.setTransform(rowToAlignment(imgRow, alignType))
+
+    if kwargs.get('readAcquisition', True):
+        img.setAcquisition(rowToAcquisition(imgRow))
+
+    if kwargs.get('magnification', None):
+        img.getAcquisition().setMagnification(kwargs.get("magnification"))
+
+    setObjId(img, imgRow)
+    # Read some extra labels
+    rowToObject(imgRow, img, {}, [])
+
+    # Provide a hook to be used if something is needed to be
+    # done for special cases before converting image to row
+    postprocessImageRow = kwargs.get('postprocessImageRow', None)
+    if postprocessImageRow:
+        postprocessImageRow(img, imgRow)
+
+    return img
+
+def rowToCoordinate(coordRow):
+    """ Create a Coordinate from a row of a metadata. """
+    # Check that all required labels are present in the row
+    if _containsAll(coordRow, COOR_DICT):
+        coord = Coordinate3D()
+        rowToObject(coordRow, coord, COOR_DICT)
+
+        # Setup the micId if is integer value
+        try:
+            coord.setTomoId(int(coordRow.getValue(emlib.MDL_MICROGRAPH_ID)))
+        except Exception:
+            pass
+    else:
+        coord = None
+
+    return coord
+
+def xmippToLocation(xmippFilename):
+    """ Return a location (index, filename) given
+    a Xmipp filename with the index@filename structure. """
+    if '@' in xmippFilename:
+        return emlib.FileName(xmippFilename).decompose()
+    else:
+        return NO_INDEX, str(xmippFilename)
+
+def rowToCtfModel(ctfRow):
+    """ Create a CTFModel from a row of a metadata. """
+    # Check if the row has CTF values, this could be called from a xmipp
+    # particles metadata
+    if _containsAll(ctfRow, CTF_DICT_NORESOLUTION):
+
+        # for compatibility reason ignore resolution and fitQuality
+        # Instantiate Scipion CTF Model
+        ctfModel = CTFTomo()
+
+        # Case for metadata coming with Xmipp resolution label
+        # Populate Scipion CTF from metadata row (using mapping dictionary
+        # plus extra labels
+        if ctfRow.hasLabel(md.MDL_CTF_PHASE_SHIFT):
+            ctfModel.setPhaseShift(ctfRow.getValue(md.MDL_CTF_PHASE_SHIFT, 0))
+        if ctfRow.containsLabel(emlib.label2Str(emlib.MDL_CTF_CRIT_MAXFREQ)):
+            rowToObject(ctfRow, ctfModel, CTF_DICT,
+                        extraLabels=CTF_EXTRA_LABELS)
+        else:
+            rowToObject(ctfRow, ctfModel, CTF_DICT_NORESOLUTION)
+
+        # Standarize defocus values
+        ctfModel.standardize()
+        # Set psd file names
+        setPsdFiles(ctfModel, ctfRow)
+        # ctfModel.setPhaseShift(0.0)  # for consistency with ctfModel
+
+    else:
+        ctfModel = None
+
+    return ctfModel
+
+def rowToAlignment(alignmentRow, alignType):
+    """
+    is2D == True-> matrix is 2D (2D images alignment)
+            otherwise matrix is 3D (3D volume alignment or projection)
+    invTransform == True  -> for xmipp implies projection
+        """
+    is2D = alignType == ALIGN_2D
+    inverseTransform = alignType == ALIGN_PROJ
+
+    if _containsAny(alignmentRow, ALIGNMENT_DICT):
+        alignment = Transform()
+        angles = np.zeros(3)
+        shifts = np.zeros(3)
+        flip = alignmentRow.getValue(emlib.MDL_FLIP)
+
+        shifts[0] = alignmentRow.getValue(emlib.MDL_SHIFT_X, 0.)
+        shifts[1] = alignmentRow.getValue(emlib.MDL_SHIFT_Y, 0.)
+        if not is2D:
+            angles[0] = alignmentRow.getValue(emlib.MDL_ANGLE_ROT, 0.)
+            angles[1] = alignmentRow.getValue(emlib.MDL_ANGLE_TILT, 0.)
+            shifts[2] = alignmentRow.getValue(emlib.MDL_SHIFT_Z, 0.)
+            angles[2] = alignmentRow.getValue(emlib.MDL_ANGLE_PSI, 0.)
+            if flip:
+                angles[1] = angles[1]+180  # tilt + 180
+                angles[2] = - angles[2]    # - psi, COSS: this is mirroring X
+                shifts[0] = -shifts[0]     # -x
+        else:
+            psi = alignmentRow.getValue(emlib.MDL_ANGLE_PSI, 0.)
+            rot = alignmentRow.getValue(emlib.MDL_ANGLE_ROT, 0.)
+            if rot != 0. and psi != 0:
+                print("HORROR rot and psi are different from zero in 2D case")
+            angles[0] = \
+                alignmentRow.getValue(emlib.MDL_ANGLE_PSI, 0.)\
+                + alignmentRow.getValue(emlib.MDL_ANGLE_ROT, 0.)
+
+        matrix = matrixFromGeometry(shifts, angles, inverseTransform)
+
+        if flip:
+            if alignType == ALIGN_2D:
+                matrix[0, :2] *= -1.  # invert only the first two columns
+                # keep x
+                matrix[2, 2] = -1.  # set 3D rot
+            elif alignType == ALIGN_3D:
+                matrix[0, :3] *= -1.  # now, invert first line excluding x
+                matrix[3, 3] *= -1.
+            elif alignType == ALIGN_PROJ:
+                pass
+
+        alignment.setMatrix(matrix)
+
+        # FIXME: now are also storing the alignment parameters since
+        # the conversions to the Transform matrix have not been extensively
+        # tested.
+        # After this, we should only keep the matrix
+        # for paramName, label in ALIGNMENT_DICT.iter():
+        #    if alignmentRow.hasLabel(label):
+        #        setattr(alignment, paramName,
+        #                alignmentRow.getValueAsObject(label))
+    else:
+        alignment = None
+
+    return alignment
+
+def rowToAcquisition(acquisitionRow):
+    """ Create an acquisition from a row of a metadata. """
+    if _containsAll(acquisitionRow, ACQUISITION_DICT):
+        acquisition = Acquisition()
+        rowToObject(acquisitionRow, acquisition, ACQUISITION_DICT)
+    else:
+        acquisition = None
+
+    return acquisition
+
+def rowToObject(row, obj, attrDict, extraLabels=[]):
+    """ This function will convert from a Row to an EMObject.
+    Params:
+        row: the Row instance (input)  -see emlib.metadata.utils.Row()-
+        obj: the EMObject instance (output)
+        attrDict: dictionary with the map between obj attributes(keys) and
+            row MDLabels in Xmipp (values).
+        extraLabels: a list with extra labels that could be included
+            as _xmipp_labelName
+    """
+    obj.setEnabled(row.getValue(emlib.MDL_ENABLED, 1) > 0)
+
+    for attr, label in attrDict.items():
+        value = row.getValue(label)
+        if not hasattr(obj, attr):
+            setattr(obj, attr, ObjectWrap(value))
+        else:
+            getattr(obj, attr).set(value)
+
+    attrLabels = attrDict.values()
+
+    for label in extraLabels:
+        if label not in attrLabels and row.hasLabel(label):
+            labelStr = emlib.label2Str(label)
+            setattr(obj, prefixAttribute(labelStr), row.getValueAsObject(label))
+
+def setObjId(obj, mdRow, label=emlib.MDL_ITEM_ID):
+    if mdRow.containsLabel(label):
+        obj.setObjId(mdRow.getValue(label))
+    else:
+        obj.setObjId(None)
+
+def _containsAll(row, labels):
+    """ Check if the labels (values) in labelsDict
+    are present in the row.
+    """
+    values = labels.values() if isinstance(labels, dict) else labels
+    return all(row.containsLabel(l) for l in values)
+
+def setPsdFiles(ctfModel, ctfRow):
+    """ Set the PSD files of CTF estimation related
+    to this ctfModel. The values will be read from
+    the ctfRow if present.
+    """
+    for attr, label in CTF_PSD_DICT.items():
+        if ctfRow.containsLabel(label):
+            setattr(ctfModel, attr, String(ctfRow.getValue(label)))
+
+def matrixFromGeometry(shifts, angles, inverseTransform):
+    """ Create the transformation matrix from a given
+    2D shifts in X and Y...and the 3 euler angles.
+    """
+    from pwem.convert import euler_matrix
+    radAngles = -np.deg2rad(angles)
+
+    M = euler_matrix(radAngles[0], radAngles[1], radAngles[2], 'szyz')
+    if inverseTransform:
+        M[:3, 3] = -shifts[:3]
+        M = np.linalg.inv(M)
+    else:
+        M[:3, 3] = shifts[:3]
+
+    return M
\ No newline at end of file
diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index 8a7044a0..bb70d84f 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -32,8 +32,7 @@
 from pwem.protocols import EMProtocol
 from pwem.constants import ALIGN_3D
 from tomo.protocols import ProtTomoBase
-from xmipptomo.convert import writeSetOfSubtomograms
-from xmipp3.convert import readSetOfParticles
+from xmipptomo.convert import writeSetOfSubtomograms, readSetOfSubtomograms
 
 OUTPUT = "output_subtomograms.xmd"
 INPUT = "input_subtomograms.xmd"
@@ -115,7 +114,7 @@ def createOutputStep(self):
         inputSubtomos = self.inputSubtomos.get()
         outputSet = self._createSetOfSubTomograms()
         outputSet.copyInfo(inputSubtomos)
-        readSetOfParticles(self._getExtraPath(OUTPUT), outputSet)
+        readSetOfSubtomograms(self._getExtraPath(OUTPUT), outputSet, alignType=ALIGN_3D)
         self._defineOutputs(outputSubtomograms=outputSet)
         self._defineSourceRelation(inputSubtomos, outputSet)
 

From abbaab6e3d4413b717a09ec0a717aa519a2384ae Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Wed, 18 Oct 2023 13:12:08 +0200
Subject: [PATCH 09/11] update program name

---
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index bb70d84f..e042fe42 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -94,7 +94,7 @@ def subtractionStep(self):
             average += ':mrc'
         resol = self.resol.get()
         iter = self.iter.get()
-        program = "xmipp_volume_subtraction"
+        program = "xmipp_subtomo_subtraction"
         args = '-i %s -o % s --ref %s --oroot %s --iter %s --lambda %s --sub --subtomos' % \
                (self._getExtraPath(INPUT), self._getExtraPath(OUTPUT), average,
                 self._getExtraPath("subtracted_subtomo"), iter, self.rfactor.get())

From cc6b27dd5a36973b042de5deaa3e2339008509dc Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Wed, 18 Oct 2023 13:25:27 +0200
Subject: [PATCH 10/11] solve sonar cloud bug

---
 xmipptomo/convert/convert.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xmipptomo/convert/convert.py b/xmipptomo/convert/convert.py
index 8bf5d04c..ed47c98e 100644
--- a/xmipptomo/convert/convert.py
+++ b/xmipptomo/convert/convert.py
@@ -559,7 +559,7 @@ def rowToAlignment(alignmentRow, alignType):
         else:
             psi = alignmentRow.getValue(emlib.MDL_ANGLE_PSI, 0.)
             rot = alignmentRow.getValue(emlib.MDL_ANGLE_ROT, 0.)
-            if rot != 0. and psi != 0:
+            if rot != 0 and psi != 0:
                 print("HORROR rot and psi are different from zero in 2D case")
             angles[0] = \
                 alignmentRow.getValue(emlib.MDL_ANGLE_PSI, 0.)\

From 070c969187ebb6181d4de6fa209d133e3e09f47b Mon Sep 17 00:00:00 2001
From: EstrellaFG <estrellafergi@gmail.com>
Date: Thu, 19 Oct 2023 15:19:36 +0200
Subject: [PATCH 11/11] remove subtomo option as it always

---
 xmipptomo/protocols/protocol_subtraction_subtomo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xmipptomo/protocols/protocol_subtraction_subtomo.py b/xmipptomo/protocols/protocol_subtraction_subtomo.py
index e042fe42..ac09d363 100644
--- a/xmipptomo/protocols/protocol_subtraction_subtomo.py
+++ b/xmipptomo/protocols/protocol_subtraction_subtomo.py
@@ -95,7 +95,7 @@ def subtractionStep(self):
         resol = self.resol.get()
         iter = self.iter.get()
         program = "xmipp_subtomo_subtraction"
-        args = '-i %s -o % s --ref %s --oroot %s --iter %s --lambda %s --sub --subtomos' % \
+        args = '-i %s -o % s --ref %s --oroot %s --iter %s --lambda %s --sub' % \
                (self._getExtraPath(INPUT), self._getExtraPath(OUTPUT), average,
                 self._getExtraPath("subtracted_subtomo"), iter, self.rfactor.get())
         if resol: