diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index bc6fdd8d..a6676261 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -1,48 +1,14 @@ -# Workflow to send master to pypi and tag the branch: -# You need to edit FOLDER_WITH_VERSION with the folder that has the __version__ value. - +# Workflow to send master to pypi and tag the branch name: master to pypi with comments and tag - -# Controls when the action will run. Triggers the workflow on push or pull request -# events but only for the master branch +# Triggers the workflow on push to the master branch on: push: branches: [ master ] -env: - FOLDER_WITH_VERSION: imod -# A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: - deploy: - - runs-on: ubuntu-latest - - steps: - - uses: actions/checkout@v2 - - name: Set up MPI - uses: mpi4py/setup-mpi@v1 - - name: Set up Python - uses: actions/setup-python@v2 - with: - python-version: '3.8' - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install setuptools wheel twine - pip install scipion-pyworkflow - pip install scipion-em - - name: Build and publish - env: - TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }} - TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }} - run: | - python setup.py sdist - twine upload dist/* -c "${{ secrets.PYPI_COMMENT }}" - - name: Get version and tag - env: - FOLDER_WITH_VERSION: imod - run: | - export PACKAGE_VERSION=$(python -c "import $FOLDER_WITH_VERSION; print('VERSION', 'v'+$FOLDER_WITH_VERSION.__version__)" | grep VERSION | sed "s/VERSION //g") - git tag $PACKAGE_VERSION - git push origin $PACKAGE_VERSION + call-publish-workflow: + uses: scipion-em/.github/.github/workflows/publish_and_tag.yml@master + with: + folder: imod + secrets: inherit diff --git a/CHANGES.txt b/CHANGES.txt index 4aec7ef1..c2cb802a 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,3 +1,19 @@ +3.5.0: + Users: + - Protocols' help added. + - Help of the parameters was enhanced. + - New feature: patch tracking alignment added. + - Trimming options in fiducial model added + - The order of the parameters in the forms was slightly modified according to the user actions. + - Protocol CTF correction: + * Now it can deal with excluded views in the CTF and/or in the TS. + * The same when entering with interpolated or not aligned TS. + * The xf, defocus, and tlt files are generated inside each TS extra sub-directory. + Developers: + - Update the acquisition order in the CTFTomo objects (field added to that class in scipion-em-tomo v3.7.0). + - Refactorization: avoid the excesive usage of getFirstItem. + - File generation methods (most of them in utils) adapted to the CTF-TS intersection functionality. + - Fixing sorted dictionary in ctf correction. How to access the defocus list in utils.py 3.3.0: - bugfix for import ctf, set missing defocus flag - move plugin-specific import CTF protocol to the core diff --git a/MANIFEST.in b/MANIFEST.in index c3d3b670..607c4755 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -4,7 +4,7 @@ include *.rst # Include the license file include LICENSE -# Include the changelog +# Include the changelog and requirements include *.txt # Include conf file diff --git a/README.rst b/README.rst index 8110fd97..901ced79 100644 --- a/README.rst +++ b/README.rst @@ -24,12 +24,6 @@ This plugin provides wrappers for several programs of `IMOD SetOfTiltSeries: """ Method to generate output classes of set of tilt-series""" - if self.TiltSeries: - self.TiltSeries.enableAppend() + outputSetOfTiltSeries = getattr(self, OUTPUT_TILTSERIES_NAME, None) + if outputSetOfTiltSeries: + outputSetOfTiltSeries.enableAppend() else: outputSetOfTiltSeries = self._createSetOfTiltSeries() @@ -270,7 +479,7 @@ def getOutputSetOfTiltSeries(self, inputSet, binning=1) -> SetOfTiltSeries: self._defineOutputs(**{OUTPUT_TILTSERIES_NAME: outputSetOfTiltSeries}) self._defineSourceRelation(inputSet, outputSetOfTiltSeries) - return self.TiltSeries + return outputSetOfTiltSeries def getOutputInterpolatedSetOfTiltSeries(self, inputSet): """ Method to generate output interpolated classes of set of tilt-series""" @@ -409,8 +618,6 @@ def getOutputSetOfTomograms(self, inputSet, binning=1): if self.Tomograms: getattr(self, OUTPUT_TOMOGRAMS_NAME).enableAppend() - - else: outputSetOfTomograms = self._createSetOfTomograms() @@ -450,7 +657,8 @@ def getOutputSetOfCTFTomoSeries(self, outputSetName): return outputSetOfCTFTomoSeries - def parseTSDefocusFile(self, inputTs, defocusFilePath, newCTFTomoSeries): + @staticmethod + def parseTSDefocusFile(inputTs, defocusFilePath, newCTFTomoSeries): """ Parse tilt-series ctf estimation file. :param inputTs: input tilt-series :param defocusFilePath: input *.defocus file to be parsed @@ -490,57 +698,57 @@ def parseTSDefocusFile(self, inputTs, defocusFilePath, newCTFTomoSeries): f"Defocus file flag {defocusFileFlag} is not supported. Only supported formats " "correspond to flags 0, 1, 4, 5, and 37.") - excludedViews = inputTs.getExcludedViewsIndex() - ids = inputTs.getIdSet() - for index in ids: + for ti in inputTs: + tiObjId = ti.getObjId() newCTFTomo = CTFTomo() - newCTFTomo.setIndex(index) + newCTFTomo.setAcquisitionOrder(ti.getAcquisitionOrder()) + newCTFTomo.setIndex(ti.getIndex()) - if index not in defocusUDict.keys() and index not in excludedViews: + if tiObjId not in defocusUDict.keys() and not ti.isEnabled(): raise IndexError("ERROR IN TILT-SERIES %s: NO CTF ESTIMATED FOR VIEW %d, TILT ANGLE %f" % ( - inputTs.getTsId(), index, inputTs[index].getTiltAngle())) + inputTs.getTsId(), tiObjId, inputTs[tiObjId].getTiltAngle())) " Plain estimation (any defocus flag)" newCTFTomo._defocusUList = CsvList(pType=float) - newCTFTomo.setDefocusUList(defocusUDict.get(index, [0.])) + newCTFTomo.setDefocusUList(defocusUDict.get(tiObjId, [0.])) if defocusFileFlag == 1: " Astigmatism estimation " newCTFTomo._defocusVList = CsvList(pType=float) - newCTFTomo.setDefocusVList(defocusVDict.get(index, [0.])) + newCTFTomo.setDefocusVList(defocusVDict.get(tiObjId, [0.])) newCTFTomo._defocusAngleList = CsvList(pType=float) - newCTFTomo.setDefocusAngleList(defocusAngleDict.get(index, [0.])) + newCTFTomo.setDefocusAngleList(defocusAngleDict.get(tiObjId, [0.])) elif defocusFileFlag == 4: " Phase-shift information " newCTFTomo._phaseShiftList = CsvList(pType=float) - newCTFTomo.setPhaseShiftList(phaseShiftDict.get(index, [0.])) + newCTFTomo.setPhaseShiftList(phaseShiftDict.get(tiObjId, [0.])) elif defocusFileFlag == 5: " Astigmatism and phase shift estimation " newCTFTomo._defocusVList = CsvList(pType=float) - newCTFTomo.setDefocusVList(defocusVDict.get(index, [0.])) + newCTFTomo.setDefocusVList(defocusVDict.get(tiObjId, [0.])) newCTFTomo._defocusAngleList = CsvList(pType=float) - newCTFTomo.setDefocusAngleList(defocusAngleDict.get(index, [0.])) + newCTFTomo.setDefocusAngleList(defocusAngleDict.get(tiObjId, [0.])) newCTFTomo._phaseShiftList = CsvList(pType=float) - newCTFTomo.setPhaseShiftList(phaseShiftDict.get(index, [0.])) + newCTFTomo.setPhaseShiftList(phaseShiftDict.get(tiObjId, [0.])) elif defocusFileFlag == 37: " Astigmatism, phase shift and cut-on frequency estimation " newCTFTomo._defocusVList = CsvList(pType=float) - newCTFTomo.setDefocusVList(defocusVDict.get(index, [0.])) + newCTFTomo.setDefocusVList(defocusVDict.get(tiObjId, [0.])) newCTFTomo._defocusAngleList = CsvList(pType=float) - newCTFTomo.setDefocusAngleList(defocusAngleDict.get(index, [0.])) + newCTFTomo.setDefocusAngleList(defocusAngleDict.get(tiObjId, [0.])) newCTFTomo._phaseShiftList = CsvList(pType=float) - newCTFTomo.setPhaseShiftList(phaseShiftDict.get(index, [0.])) + newCTFTomo.setPhaseShiftList(phaseShiftDict.get(tiObjId, [0.])) newCTFTomo._cutOnFreqList = CsvList(pType=float) - newCTFTomo.setCutOnFreqList(cutOnFreqDict.get(index, [0.])) + newCTFTomo.setCutOnFreqList(cutOnFreqDict.get(tiObjId, [0.])) newCTFTomo.completeInfoFromList() newCTFTomoSeries.append(newCTFTomo) @@ -550,21 +758,20 @@ def parseTSDefocusFile(self, inputTs, defocusFilePath, newCTFTomoSeries): newCTFTomoSeries.calculateDefocusUDeviation(defocusUTolerance=20) newCTFTomoSeries.calculateDefocusVDeviation(defocusVTolerance=20) - def createOutputFailedSet(self, tsObjId): + def createOutputFailedSet(self, ts, presentAcqOrders=None): # Check if the tilt-series ID is in the failed tilt-series # list to add it to the set - if tsObjId in self._failedTs: - ts = self._getTiltSeries(tsObjId) + tsId = ts.getTsId() + if tsId in self._failedTs: tsSet = self._getSetOfTiltSeries() - tsId = ts.getTsId() - output = self.getOutputFailedSetOfTiltSeries(tsSet) - - newTs = TiltSeries(tsId=tsId) + newTs = ts.clone() newTs.copyInfo(ts) output.append(newTs) - for index, tiltImage in enumerate(ts): + for tiltImage in ts: + if presentAcqOrders and tiltImage.getAcquisitionOrder() not in presentAcqOrders: + continue newTi = TiltImage() newTi.copyInfo(tiltImage, copyId=True, copyTM=True) newTi.setAcquisition(tiltImage.getAcquisition()) diff --git a/imod/protocols/protocol_ctfCorrection.py b/imod/protocols/protocol_ctfCorrection.py index 2d730c1e..9c719158 100644 --- a/imod/protocols/protocol_ctfCorrection.py +++ b/imod/protocols/protocol_ctfCorrection.py @@ -1,6 +1,7 @@ # ***************************************************************************** # * # * Authors: Federico P. de Isidro Gomez (fp.deisidro@cnb.csic.es) [1] +# * Scipion Team (scipion@cnb.csic.es) [1] # * # * [1] Centro Nacional de Biotecnologia, CSIC, Spain # * @@ -23,18 +24,14 @@ # * e-mail address 'scipion@cnb.csic.es' # * # ***************************************************************************** - -import os - from pyworkflow import BETA -from pyworkflow.object import Set +from pyworkflow.object import Set, String import pyworkflow.protocol.params as params -import pyworkflow.utils.path as path -import tomo.objects as tomoObj from pwem.emlib.image import ImageHandler - +from tomo.objects import TiltSeries, TiltImage +from tomo.utils import getCommonTsAndCtfElements from .. import Plugin, utils -from .protocol_base import ProtImodBase, EXT_MRCS_TS_ODD_NAME, EXT_MRCS_TS_EVEN_NAME +from .protocol_base import ProtImodBase, DEFOCUS_EXT, TLT_EXT, XF_EXT, ODD, EVEN class ProtImodCtfCorrection(ProtImodBase): @@ -42,11 +39,43 @@ class ProtImodCtfCorrection(ProtImodBase): CTF correction of a set of input tilt-series using the IMOD procedure. More info: https://bio3d.colorado.edu/imod/doc/man/ctfphaseflip.html + + This program will correct the CTF of an input tilt series by phase + flipping, with an option to attenuate frequencies near the zeros of the + CTF. + + Ctfphaseflip corrects each view strip by strip. A strip is defined as + an image region whose defocus difference is less than a user specified + value, the defocus tolerance. Normally, the strips are vertically ori- + ented and defocus is assumed to be the same along a vertical line. + Thus, the tilt series must be aligned so that the tilt axis is vertical + before applying this correction. The original thinking was that an + image region with defocus difference less than the tolerance could be + considered to have constant defocus and could be corrected as one + strip. However, the validity of the correction at the center of the + strip probably does not depend on whether it contains material beyond + this focus range, since only vertical lines near or at the center are + used in the corrected image. The program may limit the width further + to reduce computation time, or expand it to retain enough resolution + between successive zeros in the X direction of frequency space. + + Through most of the image, each strip is corrected based on the defocus + at the center of the strip. However, the strips at the left and right + edges of the image may be corrected repeatedly, at different defocus + values, in order to extend the correction close enough to the edges of + the image. + + """ _label = 'CTF correction' _devStatus = BETA + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.tsDict = None + self.ctfDict = None + # -------------------------- DEFINE param functions ----------------------- def _defineParams(self, form): form.addSection('Input') @@ -135,74 +164,54 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self._failedTs = [] - - for ts in self.inputSetOfTiltSeries.get().iterItems(): - self._insertFunctionStep(self.convertInputStep, ts.getObjId()) - self._insertFunctionStep(self.generateDefocusFile, ts.getObjId()) - self._insertFunctionStep(self.ctfCorrection, ts.getObjId()) - self._insertFunctionStep(self.createOutputStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputFailedSet, ts.getObjId()) + nonMatchingTsIds = [] + self._initialize() + for tsId in self.tsDict.keys(): # Stores the steps serializing the tsId instead of the whole ts object + matchingCtfTomoSeries = self.ctfDict.get(tsId, None) + if matchingCtfTomoSeries: + presentAcqOrders = getCommonTsAndCtfElements(self.tsDict[tsId], self.ctfDict[tsId]) + self._insertFunctionStep(self.convertInputsStep, tsId, presentAcqOrders) + self._insertFunctionStep(self.ctfCorrection, tsId) + self._insertFunctionStep(self.createOutputStep, tsId, presentAcqOrders) + self._insertFunctionStep(self.createOutputFailedStep, tsId, presentAcqOrders) + else: + nonMatchingTsIds.append(tsId) self._insertFunctionStep(self.closeOutputSetsStep) + if nonMatchingTsIds: + self.matchingMsg.set(f'WARNING! No CTFTomoSeries found for the tilt-series {nonMatchingTsIds}') + self._store(self.matchingMsg) # --------------------------- STEPS functions ----------------------------- - def convertInputStep(self, tsObjId, **kwargs): - oddEvenFlag = self.applyToOddEven(self.inputSetOfTiltSeries.get()) - - # Considering swapXY is required to make tilt axis vertical - super().convertInputStep(tsObjId, doSwap=True, oddEven=oddEvenFlag) - - def tsToProcess(self, tsObjId) -> bool: - tsId = self.inputSetOfTiltSeries.get()[tsObjId].getTsId() - ctfTomoSeries = self.getCtfTomoSeriesFromTsId(tsId) - return ctfTomoSeries - - def generateDefocusFile(self, tsObjId): - if self.tsToProcess(tsObjId) is None: - return - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - self.debug(f"Generating defocus file for {tsObjId} (ObjId), {tsId} (TsId)") - - # Compose the defocus file path - defocusFilePath = self.getDefocusFileName(ts) - - """Generate defocus file""" - ctfTomoSeries = self.getCtfTomoSeriesFromTsId(tsId) - utils.generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, inputTiltSeries=ts) - - def getDefocusFileName(self, ts): - """ Returns the path of the defocus filename based on - the tilt series and creates the folder/s""" - - tmpPrefix = self._getTmpPath(ts.getTsId()) - path.makePath(tmpPrefix) - defocusFn = ts.getFirstItem().parseFileName(extension=".defocus") - defocusFilePath = os.path.join(tmpPrefix, defocusFn) - return defocusFilePath + def _initialize(self): + self.matchingMsg = String() # Update the msg for each protocol execution to avoid duplicities in the summary + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + self.ctfDict = {ctf.getTsId(): ctf.clone(ignoreAttrs=[]) for ctf in self.inputSetOfCtfTomoSeries.get()} + + def convertInputsStep(self, tsId, presentAcqOrders): + # Generate the alignment-related files: xf, tlt, and a possible mrc + super().convertInputStep(tsId, # Considering swapXY is required to make tilt axis vertical + doSwap=True, + oddEven=self.applyToOddEven(self.inputSetOfTiltSeries.get()), + presentAcqOrders=presentAcqOrders) + # Generate the defocus file + self.generateDefocusFile(tsId, presentAcqOrders=presentAcqOrders) @ProtImodBase.tryExceptDecorator - def ctfCorrection(self, tsObjId): - if self.tsToProcess(tsObjId) is None: - return - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - firstItem = ts.getFirstItem() + def ctfCorrection(self, tsId): + ts = self.tsDict[tsId] + acq = ts.getAcquisition() """Run ctfphaseflip IMOD program""" paramsCtfPhaseFlip = { - 'inputStack': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'angleFile': os.path.join(tmpPrefix, firstItem.parseFileName(extension=".tlt")), - 'outputFileName': os.path.join(extraPrefix, firstItem.parseFileName()), - 'defocusFile': self.getDefocusFileName(ts), - 'voltage': self.inputSetOfTiltSeries.get().getAcquisition().getVoltage(), - 'sphericalAberration': self.inputSetOfTiltSeries.get().getAcquisition().getSphericalAberration(), + 'inputStack': self.getTmpOutFile(tsId), + 'angleFile': self.getExtraOutFile(tsId, ext=TLT_EXT), + 'outputFileName': self.getExtraOutFile(tsId), + 'defocusFile': self.getExtraOutFile(tsId, ext=DEFOCUS_EXT), + 'voltage': acq.getVoltage(), + 'sphericalAberration': acq.getSphericalAberration(), 'defocusTol': self.defocusTol.get(), - 'pixelSize': self.inputSetOfTiltSeries.get().getSamplingRate() / 10, - 'amplitudeContrast': self.inputSetOfTiltSeries.get().getAcquisition().getAmplitudeContrast(), + 'pixelSize': ts.getSamplingRate() / 10, # nm/px + 'amplitudeContrast': acq.getAmplitudeContrast(), 'interpolationWidth': self.interpolationWidth.get() } @@ -222,77 +231,63 @@ def ctfCorrection(self, tsObjId): "-ActionIfGPUFails 2,2 " if ts.getFirstItem().hasTransform(): - paramsCtfPhaseFlip['xformFile'] = os.path.join(tmpPrefix, firstItem.parseFileName(extension=".xf")) + paramsCtfPhaseFlip['xformFile'] = self.getExtraOutFile(tsId, ext=XF_EXT) argsCtfPhaseFlip += "-TransformFile %(xformFile)s " Plugin.runImod(self, 'ctfphaseflip', argsCtfPhaseFlip % paramsCtfPhaseFlip) if self.applyToOddEven(ts): - #ODD - oddFn = self.getTmpTSFile(tsId, tmpPrefix=tmpPrefix, suffix=EXT_MRCS_TS_ODD_NAME) - paramsCtfPhaseFlip['inputStack'] = oddFn - paramsCtfPhaseFlip['outputFileName'] = os.path.join(extraPrefix, tsId+EXT_MRCS_TS_ODD_NAME) + # ODD + paramsCtfPhaseFlip['inputStack'] = self.getTmpOutFile(tsId, suffix=ODD) + paramsCtfPhaseFlip['outputFileName'] = self.getExtraOutFile(tsId, suffix=ODD) Plugin.runImod(self, 'ctfphaseflip', argsCtfPhaseFlip % paramsCtfPhaseFlip) # EVEN - evenFn = self.getTmpTSFile(tsId, tmpPrefix=tmpPrefix, suffix=EXT_MRCS_TS_EVEN_NAME) - paramsCtfPhaseFlip['inputStack'] = evenFn - paramsCtfPhaseFlip['outputFileName'] = os.path.join(extraPrefix, tsId+EXT_MRCS_TS_EVEN_NAME) + paramsCtfPhaseFlip['inputStack'] = self.getTmpOutFile(tsId, suffix=EVEN) + paramsCtfPhaseFlip['outputFileName'] = self.getExtraOutFile(tsId, suffix=EVEN) Plugin.runImod(self, 'ctfphaseflip', argsCtfPhaseFlip % paramsCtfPhaseFlip) - def createOutputStep(self, tsObjId): - if self.tsToProcess(tsObjId) is None: - return - if tsObjId not in self._failedTs: - inputTs = self.inputSetOfTiltSeries.get() - output = self.getOutputSetOfTiltSeries(inputTs) - hasAlign = inputTs.getFirstItem().getFirstItem().hasTransform() + def createOutputStep(self, tsId, presentAcqOrders): + if tsId not in self._failedTs: + inTsSet = self.inputSetOfTiltSeries.get() + outputSetOfTs = self.getOutputSetOfTiltSeries(inTsSet) - ts = inputTs[tsObjId] - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - - newTs = tomoObj.TiltSeries(tsId=tsId) + newTs = TiltSeries(tsId=tsId) + ts = self.tsDict[tsId] newTs.copyInfo(ts) newTs.setCtfCorrected(True) newTs.setInterpolated(True) - output.append(newTs) + acq = newTs.getAcquisition() + acq.setTiltAxisAngle(0.) # 0 because TS is aligned + newTs.setAcquisition(acq) + outputSetOfTs.append(newTs) ih = ImageHandler() - for index, tiltImage in enumerate(ts): - newTi = tomoObj.TiltImage() - newTi.copyInfo(tiltImage, copyId=True, copyTM=False) - acq = tiltImage.getAcquisition() - if hasAlign: - acq.setTiltAxisAngle(0.) - newTi.setAcquisition(acq) - newTi.setLocation(index + 1, - (os.path.join(extraPrefix, - tiltImage.parseFileName()))) - if self.applyToOddEven(ts): - locationOdd = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_ODD_NAME)) - locationEven = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_EVEN_NAME)) - newTi.setOddEven([ih.locationToXmipp(locationOdd), ih.locationToXmipp(locationEven)]) - else: - newTi.setOddEven([]) - - newTs.append(newTi) - - if hasAlign: - acq = newTs.getAcquisition() - acq.setTiltAxisAngle(0.) # 0 because TS is aligned - newTs.setAcquisition(acq) + for index, inTi in enumerate(ts): + if inTi.getAcquisitionOrder() in presentAcqOrders: + newTi = TiltImage() + newTi.copyInfo(inTi, copyId=True, copyTM=False) + acq = inTi.getAcquisition() + acq.setTiltAxisAngle(0.) # Is interpolated + newTi.setAcquisition(acq) + newTi.setLocation(index + 1, self.getExtraOutFile(tsId)) + if self.applyToOddEven(ts): + locationOdd = index + 1, self.getExtraOutFile(tsId, suffix=ODD) + locationEven = index + 1, self.getExtraOutFile(tsId, suffix=EVEN) + newTi.setOddEven([ih.locationToXmipp(locationOdd), ih.locationToXmipp(locationEven)]) + else: + newTi.setOddEven([]) + newTs.append(newTi) newTs.write(properties=False) - output.update(newTs) - output.write() - self._store() + outputSetOfTs.update(newTs) + outputSetOfTs.write() + self._store(outputSetOfTs) - def createOutputFailedSet(self, tsObjId): - if self.tsToProcess(tsObjId) is None: - return - super().createOutputFailedSet(tsObjId) + def createOutputFailedStep(self, tsId, presentAcqOrders): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts, presentAcqOrders=presentAcqOrders) def closeOutputSetsStep(self): for _, output in self.iterOutputAttributes(): @@ -301,11 +296,17 @@ def closeOutputSetsStep(self): self._store() # --------------------------- UTILS functions ----------------------------- - def getCtfTomoSeriesFromTsId(self, tsId): - for ctfTomoSeries in self.inputSetOfCtfTomoSeries.get(): - if tsId == ctfTomoSeries.getTsId(): - return ctfTomoSeries - return None + def generateDefocusFile(self, tsId, presentAcqOrders=None): + ts = self.tsDict[tsId] + ctfTomoSeries = self.ctfDict[tsId] + + self.debug(f"Generating defocus file for {tsId} (ObjId), {tsId} (TsId)") + # Compose the defocus file path + defocusFilePath = self.getExtraOutFile(tsId, ext=DEFOCUS_EXT) + """Generate defocus file""" + utils.generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, + inputTiltSeries=ts, + presentAcqOrders=presentAcqOrders) # --------------------------- INFO functions ------------------------------ def _warnings(self): @@ -328,6 +329,8 @@ def _summary(self): self.TiltSeries.getSize())) else: summary.append("Outputs are not ready yet.") + if self.matchingMsg.get(): + summary.append(self.matchingMsg.get()) return summary def _methods(self): diff --git a/imod/protocols/protocol_ctfEstimation_automatic.py b/imod/protocols/protocol_ctfEstimation_automatic.py index 6b57507f..14fbefcf 100644 --- a/imod/protocols/protocol_ctfEstimation_automatic.py +++ b/imod/protocols/protocol_ctfEstimation_automatic.py @@ -32,7 +32,7 @@ import tomo.objects as tomoObj from .. import Plugin, utils -from .protocol_base import ProtImodBase, OUTPUT_CTF_SERIE +from .protocol_base import ProtImodBase, OUTPUT_CTF_SERIE, TLT_EXT, DEFOCUS_EXT class ProtImodAutomaticCtfEstimation(ProtImodBase): @@ -49,13 +49,18 @@ class ProtImodAutomaticCtfEstimation(ProtImodBase): defocusVTolerance = 20 _interactiveMode = False + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.sRate = None + self.acq = None + # -------------------------- DEFINE param functions ----------------------- def _defineParams(self, form): form.addSection('Input') form.addParam('inputSet', params.PointerParam, pointerClass='SetOfTiltSeries, SetOfCTFTomoSeries', - label='Input set of tilt-series', + label='Tilt-series', help='This should be a *raw stack*, not an aligned stack, ' 'because the interpolation used to make ' 'an aligned stack attenuates high frequencies and ' @@ -292,45 +297,41 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self._failedTs = [] - - # This assignment is needed to use methods from base class - self.inputSetOfTiltSeries = self._getSetOfTiltSeries() + self._initialize() expDefoci = self.getExpectedDefocus() - - for item in self.inputSet.get(): - tsObjId = item.getObjId() - self._insertFunctionStep(self.convertInputStep, tsObjId) - self._insertFunctionStep(self.ctfEstimation, tsObjId, expDefoci) - self._insertFunctionStep(self.createOutputStep, tsObjId) - self._insertFunctionStep(self.createOutputFailedSet, tsObjId) + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + self._insertFunctionStep(self.ctfEstimation, tsId, expDefoci) + self._insertFunctionStep(self.createOutputStep, tsId) + self._insertFunctionStep(self.createOutputFailedStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def convertInputStep(self, tsObjId, **kwargs): + def _initialize(self): + self._failedTs = [] + tsSet = self._getSetOfTiltSeries() + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in tsSet} + self.sRate = tsSet.getSamplingRate() + self.acq = tsSet.getAcquisition() + + def convertInputStep(self, tsId, **kwargs): """ Implement the convertStep to cancel interpolation of the tilt series.""" - super().convertInputStep(tsObjId, imodInterpolation=None) + super().convertInputStep(tsId, imodInterpolation=None) @ProtImodBase.tryExceptDecorator - def ctfEstimation(self, tsObjId, expDefoci): + def ctfEstimation(self, tsId, expDefoci): """Run ctfplotter IMOD program""" - ts = self._getTiltSeries(tsObjId) - tsSet = self._getSetOfTiltSeries() - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - - firstTi = ts.getFirstItem() + ts = self.tsDict[tsId] paramsCtfPlotter = { - 'inputStack': os.path.join(tmpPrefix, firstTi.parseFileName()), - 'angleFile': os.path.join(tmpPrefix, firstTi.parseFileName(extension=".tlt")), - 'defocusFile': os.path.join(extraPrefix, firstTi.parseFileName(extension=".defocus")), + 'inputStack': self.getTmpOutFile(tsId), + 'angleFile': self.getExtraOutFile(tsId, ext=TLT_EXT), + 'defocusFile': self.getExtraOutFile(tsId, ext=DEFOCUS_EXT), 'axisAngle': ts.getAcquisition().getTiltAxisAngle(), - 'pixelSize': tsSet.getSamplingRate() / 10, - 'voltage': tsSet.getAcquisition().getVoltage(), - 'sphericalAberration': tsSet.getAcquisition().getSphericalAberration(), - 'amplitudeContrast': tsSet.getAcquisition().getAmplitudeContrast(), + 'pixelSize': self.sRate / 10, # nm + 'voltage': self.acq.getVoltage(), + 'sphericalAberration': self.acq.getSphericalAberration(), + 'amplitudeContrast': self.acq.getAmplitudeContrast(), 'defocusTol': self.defocusTol.get(), 'psResolution': 101, 'leftDefTol': self.leftDefTol.get(), @@ -434,13 +435,9 @@ def ctfEstimation(self, tsObjId, expDefoci): Plugin.runImod(self, 'ctfplotter', argsCtfPlotter % paramsCtfPlotter) - def createOutputStep(self, tsObjId, outputSetName=OUTPUT_CTF_SERIE): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - defocusFilePath = os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".defocus")) + def createOutputStep(self, tsId, outputSetName=OUTPUT_CTF_SERIE): + ts = self.tsDict[tsId] + defocusFilePath = self.getExtraOutFile(tsId, ext=DEFOCUS_EXT) if os.path.exists(defocusFilePath): output = self.getOutputSetOfCTFTomoSeries(outputSetName) defocusFileFlag = utils.getDefocusFileFlag(defocusFilePath) @@ -448,7 +445,7 @@ def createOutputStep(self, tsObjId, outputSetName=OUTPUT_CTF_SERIE): newCTFTomoSeries = tomoObj.CTFTomoSeries() newCTFTomoSeries.copyInfo(ts) newCTFTomoSeries.setTiltSeries(ts) - newCTFTomoSeries.setObjId(tsObjId) + # newCTFTomoSeries.setObjId(tsObjId) newCTFTomoSeries.setTsId(tsId) newCTFTomoSeries.setIMODDefocusFileFlag(defocusFileFlag) newCTFTomoSeries.setNumberOfEstimationsInRange(None) @@ -464,6 +461,10 @@ def createOutputStep(self, tsObjId, outputSetName=OUTPUT_CTF_SERIE): output.write() self._store() + def createOutputFailedStep(self, tsId): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts) + def closeOutputSetsStep(self): for _, output in self.iterOutputAttributes(): output.setStreamState(Set.STREAM_CLOSED) @@ -498,6 +499,7 @@ def getExpectedDefocus(self): if self.expectedDefocusOrigin.get() == 1: with open(self.expectedDefocusFile.get()) as f: lines = f.readlines() + lines = filter(lambda x: x.strip(), lines) result = {line.split()[0]: line.split()[1] for line in lines} return result else: diff --git a/imod/protocols/protocol_ctfEstimation_manual.py b/imod/protocols/protocol_ctfEstimation_manual.py index 10b9b74d..95fb7b07 100644 --- a/imod/protocols/protocol_ctfEstimation_manual.py +++ b/imod/protocols/protocol_ctfEstimation_manual.py @@ -40,6 +40,13 @@ class ProtImodManualCtfEstimation(ProtImodAutomaticCtfEstimation): More info: https://bio3d.colorado.edu/imod/doc/man/ctfplotter.html + This GUI program will plot the logarithm of a rotationally averaged + power spectrum of an input tilt series after subtracting the noise + floor. The method is based on periodogram averaging; namely, averaging + of spectra from small, overlapping areas of the images, referred to as + tiles. The user can interactively choose which projection views are + included in the averaging. It is also possible to run the program non- + interactively for automatic fitting. """ _label = 'CTF estimation (manual)' diff --git a/imod/protocols/protocol_doseFilter.py b/imod/protocols/protocol_doseFilter.py index fa2d1318..29e93a41 100644 --- a/imod/protocols/protocol_doseFilter.py +++ b/imod/protocols/protocol_doseFilter.py @@ -23,18 +23,14 @@ # * e-mail address 'scipion@cnb.csic.es' # * # ***************************************************************************** - -import os - from pyworkflow import BETA from pyworkflow.object import Set import pyworkflow.protocol.params as params -import pyworkflow.utils.path as path import tomo.objects as tomoObj from pwem.emlib.image import ImageHandler from .. import Plugin, utils -from .protocol_base import ProtImodBase, EXT_MRCS_TS_EVEN_NAME, EXT_MRCS_TS_ODD_NAME +from .protocol_base import ProtImodBase, ODD, EVEN SCIPION_IMPORT = 0 FIXED_DOSE = 1 @@ -45,6 +41,16 @@ class ProtImodDoseFilter(ProtImodBase): Tilt-series dose filtering based on the IMOD procedure. More info: https://bio3d.colorado.edu/imod/doc/man/mtffilter.html + + A specialized filter can be applied to perform dose weight-filtering of + cryoEM images, particularly ones from tilt series. The filter is as + described in Grant and Grigorieff, 2015 (DOI: 10.7554/eLife.06980) and + the implementation follows that in their "unblur" program. At any fre- + quency, the filter follows an exponential decay with dose, where the + exponential is of the dose divided by 2 times a "critical dose" for + that frequency. This critical dose was empirically found to be approx- + imated by a * k^b + c, where k is frequency; the values of a, b, c in + that paper are used by default. """ _label = 'Dose filter' @@ -58,13 +64,15 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTiltSeries', important=True, - label='Input set of tilt-series') + label='Tilt Series', + help='This input tilt-series will be low pass filtered according' + 'to their cumulated dose') form.addParam('initialDose', params.FloatParam, default=0.0, expertLevel=params.LEVEL_ADVANCED, - label='Initial dose (e/sq. Å)', + label='Initial dose (e/Å^2)', help='Dose applied before any of the images in the ' 'input file were taken; this value will be ' 'added to all the dose values.') @@ -81,12 +89,12 @@ def _defineParams(self, form): 'during import of the tilt-series\n' '- Fixed dose: manually input fixed dose ' 'for each image of the input file, ' - 'in electrons/square Ångstrom.') + 'in electrons/Å^2.') form.addParam('fixedImageDose', params.FloatParam, default=FIXED_DOSE, - label='Fixed dose (e/sq Å)', + label='Fixed dose (e/Å^2)', condition='inputDoseType == %i' % FIXED_DOSE, help='Fixed dose for each image of the input file, ' 'in electrons/square Ångstrom.') @@ -101,29 +109,26 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.doseFilterStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputStep, ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.doseFilterStep, tsId) + self._insertFunctionStep(self.createOutputStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def doseFilterStep(self, tsObjId): - """Apply the dose filter to every tilt series""" - - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() + def _initialize(self): + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - - path.makePath(tmpPrefix) - path.makePath(extraPrefix) + def doseFilterStep(self, tsId): + """Apply the dose filter to every tilt series""" + self.genTsPaths(tsId) + ts = self.tsDict[tsId] firstItem = ts.getFirstItem() paramsMtffilter = { 'input': firstItem.getFileName(), - 'output': os.path.join(extraPrefix, firstItem.parseFileName()), + 'output': self.getExtraOutFile(tsId), 'pixsize': ts.getSamplingRate(), 'voltage': ts.getAcquisition().getVoltage(), } @@ -137,9 +142,7 @@ def doseFilterStep(self, tsObjId): argsMtffilter += f"-InitialDose {self.initialDose.get():f} " if self.inputDoseType.get() == SCIPION_IMPORT: - outputDoseFilePath = os.path.join(tmpPrefix, - firstItem.parseFileName(extension=".dose")) - + outputDoseFilePath = self.getExtraOutFile(tsId, ext="dose") utils.generateDoseFileFromDoseTS(ts, outputDoseFilePath) paramsMtffilter.update({ @@ -162,27 +165,22 @@ def doseFilterStep(self, tsObjId): if self.applyToOddEven(ts): oddFn = firstItem.getOdd().split('@')[1] paramsMtffilter['input'] = oddFn - paramsMtffilter['output'] = os.path.join(extraPrefix, tsId+EXT_MRCS_TS_ODD_NAME) + paramsMtffilter['output'] = self.getExtraOutFile(tsId, suffix=ODD) Plugin.runImod(self, 'mtffilter', argsMtffilter % paramsMtffilter) evenFn = firstItem.getEven().split('@')[1] paramsMtffilter['input'] = evenFn - paramsMtffilter['output'] = os.path.join(extraPrefix, tsId+EXT_MRCS_TS_EVEN_NAME) + paramsMtffilter['output'] = self.getExtraOutFile(tsId, suffix=EVEN) Plugin.runImod(self, 'mtffilter', argsMtffilter % paramsMtffilter) - def createOutputStep(self, tsObjId): + def createOutputStep(self, tsId): """Generate output filtered tilt series""" - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - + ts = self.tsDict[tsId] output = self.getOutputSetOfTiltSeries(self.inputSetOfTiltSeries.get()) newTs = tomoObj.TiltSeries(tsId=tsId) newTs.copyInfo(ts) - output.append(newTs) - ih = ImageHandler() for index, tiltImage in enumerate(ts): @@ -190,14 +188,13 @@ def createOutputStep(self, tsObjId): newTi.copyInfo(tiltImage, copyId=True, copyTM=True) newTi.setAcquisition(tiltImage.getAcquisition()) if self.applyToOddEven(ts): - locationOdd = index + 1, (os.path.join(extraPrefix, tsId+EXT_MRCS_TS_ODD_NAME)) - locationEven = index + 1, (os.path.join(extraPrefix, tsId+EXT_MRCS_TS_EVEN_NAME)) + locationOdd = index + 1, self.getExtraOutFile(tsId, suffix=ODD) + locationEven = index + 1, self.getExtraOutFile(tsId, suffix=EVEN) newTi.setOddEven([ih.locationToXmipp(locationOdd), ih.locationToXmipp(locationEven)]) else: newTi.setOddEven([]) - locationTi = index + 1, (os.path.join(extraPrefix, - tiltImage.parseFileName())) + locationTi = index + 1, self.getExtraOutFile(tsId) newTi.setLocation(locationTi) newTs.append(newTi) newTs.update(newTi) diff --git a/imod/protocols/protocol_etomo.py b/imod/protocols/protocol_etomo.py index 0a05ef2c..cd9dd90d 100644 --- a/imod/protocols/protocol_etomo.py +++ b/imod/protocols/protocol_etomo.py @@ -46,6 +46,17 @@ class ProtImodEtomo(ProtImodBase): More info: https://bio3d.colorado.edu/imod/doc/etomoTutorial.html + + Etomo is software tool for assisting users in the tomographic reconstruction + process of both single and dual axis tilt series. Throughout this procedure, + eTomo executes numerous program commands and frequently initiates 3dmod + and Midas to enable users to make precise adjustments. Some of the main features + are:\n + - Xray eraser\n + - dose filtering\n + - Tilt series alignment\n + - Gold beads detection and eraser\n + - Tomogram reconstruction\n """ _label = 'Etomo interactive' @@ -56,7 +67,7 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self.PrealignedTiltSeries = None self.FullTomograms = None - self.PostProcessedTomograms = None + self.PostProcessTomograms = None # -------------------------- DEFINE param functions ----------------------- def _defineParams(self, form): @@ -73,7 +84,9 @@ def _defineParams(self, form): params.FloatParam, default=10, label='Fiducial markers diameter (nm)', - help='Diameter of gold beads in nanometers.') + help='Diameter of gold beads in nanometers. Note that fiducials are' + 'small gold beads that are used as marker elements in the images.' + 'They can be used as reference points to align the tilt series') form.addParam('applyAlignment', params.BooleanParam, @@ -632,7 +645,7 @@ def _summary(self): if self.FullTomograms: summary.append("- Raw reconstructed tomogram") - if self.PostProcessedTomograms: + if self.PostProcessTomograms: summary.append("- Post-processed tomogram") if len(summary) == 1: diff --git a/imod/protocols/protocol_fiducialAlignment.py b/imod/protocols/protocol_fiducialAlignment.py index 0110a2d5..473eb15f 100644 --- a/imod/protocols/protocol_fiducialAlignment.py +++ b/imod/protocols/protocol_fiducialAlignment.py @@ -36,7 +36,7 @@ TiltSeries, TiltSeriesCoordinate) from .. import Plugin, utils -from .protocol_base import ProtImodBase +from .protocol_base import ProtImodBase, TLT_EXT, XF_EXT, FID_EXT, TXT_EXT, XYZ_EXT, MOD_EXT, SFID_EXT class ProtImodFiducialAlignment(ProtImodBase): @@ -45,10 +45,70 @@ class ProtImodFiducialAlignment(ProtImodBase): on the IMOD procedure. More info: https://bio3d.colorado.edu/imod/doc/man/tiltalign.html - https://bio3d.colorado.edu/imod/doc/man/model2point.html - https://bio3d.colorado.edu/imod/doc/man/imodtrans.html - https://bio3d.colorado.edu/imod/doc/man/newstack.html - https://bio3d.colorado.edu/imod/doc/man/ccderaser.html + + This program will solve for the displacements, rotations, tilts, and + magnification differences relating a set of tilted views of an object. + It uses a set of fiducial points that have been identified in a series + of views. These input data are read from a model in which each fiducial + point is a separate contour. + + This program has several notable features: + + 1) Any given fiducial point need not be present in every view. Thus, + one can track each fiducial point only through the set of views in + which it can be reliably identified, and one can even skip views in the + middle of that set. + + 2) The program can solve for distortion (stretching) in the plane of + the section. + + 3) It is possible to constrain several views to have the same unknown + value of rotation, tilt angle, magnification, compression, or distor- + tion. This can reduce the number of unknowns and can give more accu- + rate overall solutions. + + 4) If the fiducial points are supposed to lie in one or two planes, + then after the minimization procedure is complete, the program can ana- + lyze the solved point positions and determine the slope of this plane. + It uses this slope to estimate how to adjust tilt angles so as to make + the planes be horizontal in a reconstruction. + + 5) The program can use a robust fitting method to give different + weights to different modeled points based on their individual fitting + errors. Points with the most extreme errors are eliminated from the + fit, and ones with high but less extreme errors are down-weighted. + This fitting provides a substitute for fixing many modeled point posi- + tions based on their errors. + + _On the alignment model_:\n + The program implements the following model for the imaging of the spec- + imen in each individual view: + 1) The specimen itself changes by + a) an isotropic size change (magnification variable); + b) additional thinning in the Z dimension (compression variable); and + c) linear stretch along one axis in the specimen plane, implemented by + variables representing stretch along the X axis and skew between + the X and Y axes; + 2) The specimen is tilted slightly around the X axis (X tilt variable) + 3) The specimen is tilted around the X axis by the negative of the beam + tilt, if any (one variable for all views) + 4) The specimen is tilted around the Y axis (tilt variable) + 5) The specimen is tilted back around the X axis by the beam tilt, if any + 6) The projected image rotates in the plane of the camera (rotation + variable) + 7) The projected image may stretch along an axis midway between the + original X and Y axes (one variable for all views) + 8) The image shifts on the camera + + The complete model is summarized in: + Mastronarde, D. N. 2008. Correction for non-perpendicularity of beam + and tilt axis in tomographic reconstructions with the IMOD package. J. + Microsc. 230: 212-217. + The version of the model prior to the addition of beam tilt is described + in more detail in: + Mastronarde, D. N. 2007. Fiducial marker and hybrid alignment methods + for single- and double-axis tomography. In: Electron Tomography, Ed. + J. Frank, 2nd edition, pp 163-185. Springer, New York. """ _label = 'Fiducial alignment' @@ -63,7 +123,7 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfLandmarkModels', important=True, - label='Input set of fiducial models.') + label='Fiducial model') # TODO: Allow for a different set of tilt-series input source than the one from the landmark model. This is not # TODO: possible due to a change of data type when applying the transformation with scipion applyTransform @@ -91,17 +151,14 @@ def _defineParams(self, form): params.EnumParam, choices=['Yes', 'No'], default=1, - label='Find beads on two surfaces?', + label='Assume beads on two surfaces?', display=params.EnumParam.DISPLAY_HLIST, - help="Track fiducials differentiating in which side of " - "the sample are located.\nIMPORTANT: It is highly " - "recommended to match the option selected in the " - "generation of the fiducial models. In case they " - "do not match, it is not intended to fail but could " - "be missing the whole potential of the algorithm. " - "In case the algorithm used fot he calculation" - "of the fiducial models does not consider this " - "option it is algo recomended to set this " + help="Track fiducials differentiating in which side of the sample are located.\n" + "IMPORTANT: It is highly recommended to match the option selected in the " + "generation of the fiducial models. In case they do not match, it is not " + "intended to fail but could be missing the whole potential of the algorithm. " + "In case the algorithm used for the calculation of the fiducial models does " + "not consider this option it is algo recommended to set this " "option to 'No'.") form.addParam('computeAlignment', @@ -111,8 +168,12 @@ def _defineParams(self, form): label='Generate interpolated tilt-series?', important=True, display=params.EnumParam.DISPLAY_HLIST, - help='Generate and save the interpolated tilt-series ' - 'applying the obtained transformation matrices.') + help='Generate and save the interpolated tilt-series applying the obtained transformation ' + 'matrices.\n' + 'By default, the output of this protocol will be a tilseries that will have associated' + 'the alignment information as a transformation matrix. When this option is set as Yes, ' + 'then a second output, called interpolated tilt series, is generated. The interpolated tilt ' + 'series should be used for visualization purpose but not for image processing') groupInterpolation = form.addGroup('Interpolated tilt-series', condition='computeAlignment==0') @@ -120,12 +181,14 @@ def _defineParams(self, form): groupInterpolation.addParam('binning', params.IntParam, default=1, - label='Binning', - help='Binning to be applied to the ' - 'interpolated tilt-series in IMOD ' - 'convention. Images will be binned ' - 'by the given factor. Must be an ' - 'integer bigger than 1') + label='Binning for the interpolated', + help='Binning to be applied to the interpolated tilt-series in IMOD ' + 'convention. \n' + 'Binning is an scaling factor given by an integer greater than 1. ' + 'IMOD uses ordinary binning to reduce images in size by the given factor. ' + 'The value of a binned pixel is the average of pixel values in each block ' + 'of pixels being binned. Binning is applied before all other image ' + 'transformations.') form.addSection('Global variables') @@ -136,14 +199,19 @@ def _defineParams(self, form): default=3, label='Rotation solution type', display=params.EnumParam.DISPLAY_HLIST, - help='Type of rotation solution.') + help='Type of rotation solution: See rotOption in tiltalign IMOD command \n' + '* No rotation: The in-plane rotation will not be estimated\n' + '* One rotation: To solve for a single rotation variable \n' + '* Group rotations: Group views to solve for fewer rotations variables. Automapping of ' + 'rotation variables linearly changing values\n' + '* Solve for all rotations: for each view having an independent rotation\n') form.addParam('groupRotationSize', params.IntParam, default=5, condition='rotationSolutionType==2', label='Group size', - help='Size of the rotation group') + help='Default group size when automapping rotation variables') form.addParam('magnificationSolutionType', params.EnumParam, @@ -153,14 +221,18 @@ def _defineParams(self, form): default=1, label='Magnification solution type', display=params.EnumParam.DISPLAY_HLIST, - help='Type of magnification solution.') + help='Type of magnification solution: See MagOption in tiltaling IMOD command\n' + '* Fixed magnification: Do not solve magnification. This fixes all magnifications at 1.0.\n' + '* Group magnifications: Group views to solve for fewer magnifications variables. ' + 'Automapping of variables (linearly changing values) \n' + '* Solve for all magnifications: to vary all magnifications of each view independently\n') form.addParam('groupMagnificationSize', params.IntParam, default=4, condition='magnificationSolutionType==1', label='Group size', - help='Size of the magnification group') + help='Group size when automapping magnification variables') form.addParam('tiltAngleSolutionType', params.EnumParam, @@ -169,14 +241,18 @@ def _defineParams(self, form): default=1, label='Tilt angle solution type', display=params.EnumParam.DISPLAY_HLIST, - help='Type of tilt angle solution.') + help='Type of tilt angle solution: See TiltOption in tiltalign IMOD command\n' + ' * Fixed tilt angles: To fix all tilt angles at their initial (input) values \n' + ' * Group tilt angles: To automap groups of tilt angles (linearly changing values) \n' + ' * Solve for all except minimum tilt:to solve for all tilt angles except for the view ' + 'at minimum tilt \n') form.addParam('groupTiltAngleSize', params.IntParam, default=5, condition='tiltAngleSolutionType==1', label='Group size', - help='Size of the tilt angle group') + help='Average default group size when automapping tilt variables') form.addParam('distortionSolutionType', params.EnumParam, @@ -184,7 +260,12 @@ def _defineParams(self, form): default=0, label='Distortion solution type', display=params.EnumParam.DISPLAY_HLIST, - help='Type of distortion solution.') + help='Type of skew solution:' + '* 0 to fix all skew angles at 0.0 \n' + '* 1 to vary all skew angles independently\n ' + '* 2 to specify a mapping of skew variables, or\n ' + '* 3 or 4 for automapping of variables (3 for linearly changing values or 4 for values all ' + 'the same within a group)..') form.addParam('xStretchGroupSize', params.IntParam, @@ -236,80 +317,59 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self.inputSetOfTiltSeries = self.inputSetOfLandmarkModels.get().getSetOfTiltSeries(pointer=True) - - tsIds = self.inputSetOfLandmarkModels.get().aggregate(["COUNT"], "_tsId", ["_tsId"]) - tsIds = set([d['_tsId'] for d in tsIds]) - - tsIdsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in - self.inputSetOfTiltSeries.get() if - ts.getTsId() in tsIds} - - self._failedTs = [] - + self._initialize() for lm in self.inputSetOfLandmarkModels.get(): lmTsId = lm.getTsId() self.fiducialDiameterPixel = lm.getSize() - - tsObjId = tsIdsDict[lmTsId].getObjId() - self._insertFunctionStep(self.convertInputStep, tsObjId) - self._insertFunctionStep(self.computeFiducialAlignmentStep, tsObjId) - self._insertFunctionStep(self.translateFiducialPointModelStep, tsObjId) - self._insertFunctionStep(self.computeOutputStackStep, tsObjId) + self._insertFunctionStep(self.convertInputStep, lmTsId) + self._insertFunctionStep(self.computeFiducialAlignmentStep, lmTsId) + self._insertFunctionStep(self.translateFiducialPointModelStep, lmTsId) + self._insertFunctionStep(self.computeOutputStackStep, lmTsId) if self.computeAlignment.get() == 0 or self.eraseGoldBeads.get() == 0: self._insertFunctionStep(self.computeOutputInterpolatedStackStep, - tsObjId) + lmTsId) if self.eraseGoldBeads.get() == 0: - self._insertFunctionStep(self.eraseGoldBeadsStep, tsObjId) + self._insertFunctionStep(self.eraseGoldBeadsStep, lmTsId) - self._insertFunctionStep(self.computeOutputModelsStep, tsObjId) - self._insertFunctionStep(self.createOutputFailedSet, tsObjId) + self._insertFunctionStep(self.computeOutputModelsStep, lmTsId) + self._insertFunctionStep(self.createOutputFailedStep, lmTsId) self._insertFunctionStep(self.createOutputStep) # --------------------------- STEPS functions ----------------------------- - @ProtImodBase.tryExceptDecorator - def computeFiducialAlignmentStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() + def _initialize(self): + self.inputSetOfTiltSeries = self.inputSetOfLandmarkModels.get().getSetOfTiltSeries(pointer=True) - lm = self.inputSetOfLandmarkModels.get().getLandmarkModelFromTsId(tsId=tsId) + tsIds = self.inputSetOfLandmarkModels.get().aggregate(["COUNT"], "_tsId", ["_tsId"]) + tsIds = set([d['_tsId'] for d in tsIds]) - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in + self.inputSetOfTiltSeries.get() if + ts.getTsId() in tsIds} - firstItem = ts.getFirstItem() + self._failedTs = [] + + @ProtImodBase.tryExceptDecorator + def computeFiducialAlignmentStep(self, tsId): + ts = self.tsDict[tsId] + lm = self.inputSetOfLandmarkModels.get().getLandmarkModelFromTsId(tsId=tsId) paramsTiltAlign = { 'modelFile': lm.getModelName(), - 'imageFile': os.path.join(tmpPrefix, firstItem.parseFileName()), + 'imageFile': self.getTmpOutFile(tsId), 'imagesAreBinned': 1, 'unbinnedPixelSize': ts.getSamplingRate() / 10, - 'outputModelFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fidxyz", - extension=".mod")), - 'outputResidualFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_resid", - extension=".txt")), - 'outputFidXYZFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".xyz")), - 'outputTiltFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_interpolated", - extension=".tlt")), - 'outputXAxisTiltFile': os.path.join(extraPrefix, - firstItem.parseFileName(extension=".xtilt")), - 'outputTransformFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".xf")), - 'outputFilledInModel': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_noGaps", - extension=".fid")), + 'outputModelFile': self.getExtraOutFile(tsId, suffix="fidxyz", ext=MOD_EXT), + 'outputResidualFile': self.getExtraOutFile(tsId, suffix="resid", ext=TXT_EXT), + 'outputFidXYZFile': self.getExtraOutFile(tsId, suffix="fid", ext=XYZ_EXT), + 'outputTiltFile': self.getExtraOutFile(tsId, suffix="interpolated", ext=TLT_EXT), + 'outputXAxisTiltFile': self.getExtraOutFile(tsId, ext="xtilt"), + 'outputTransformFile': self.getExtraOutFile(tsId, suffix="fid", ext=XF_EXT), + 'outputFilledInModel': self.getExtraOutFile(tsId, suffix="noGaps", ext=FID_EXT), 'rotationAngle': ts.getAcquisition().getTiltAxisAngle(), - 'tiltFile': os.path.join(tmpPrefix, - firstItem.parseFileName(extension=".tlt")), + 'tiltFile': self.getExtraOutFile(tsId, ext=TLT_EXT), 'angleOffset': 0.0, 'rotOption': self.getRotationType(), 'rotDefaultGrouping': self.groupRotationSize.get(), @@ -334,9 +394,7 @@ def computeFiducialAlignmentStep(self, tsObjId): 'axisZShift': 0.0, 'shiftZFromOriginal': 1, 'localAlignments': 0, - 'outputLocalFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_local", - extension=".xf")), + 'outputLocalFile': self.getExtraOutFile(tsId, suffix="local", ext=XF_EXT), 'targetPatchSizeXandY': '700,700', 'minSizeOrOverlapXandY': '0.5,0.5', 'minFidsTotalAndEachSurface': '8,3', @@ -353,7 +411,7 @@ def computeFiducialAlignmentStep(self, tsObjId): 'localXStretchDefaultGrouping': 7, 'localSkewOption': 0, 'localSkewDefaultGrouping': 11, - 'outputTiltAlignFileText': os.path.join(extraPrefix, "align.log"), + 'outputTiltAlignFileText': self._getExtraPath("align.log"), } argsTiltAlign = "-ModelFile %(modelFile)s " \ @@ -420,28 +478,16 @@ def computeFiducialAlignmentStep(self, tsObjId): argsTiltAlign += "2>&1 | tee %(outputTiltAlignFileText)s " Plugin.runImod(self, 'tiltalign', argsTiltAlign % paramsTiltAlign) - Plugin.runImod(self, 'alignlog', '-s > taSolution.log', cwd=extraPrefix) + Plugin.runImod(self, 'alignlog', '-s > taSolution.log', cwd=self._getExtraPath()) @ProtImodBase.tryExceptDecorator - def translateFiducialPointModelStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - - firstItem = ts.getFirstItem() - + def translateFiducialPointModelStep(self, tsId): # Check that previous steps have been completed satisfactorily - if os.path.exists(os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_noGaps", - extension=".fid"))): + noGapsFid = self.getExtraOutFile(tsId, suffix="noGaps", ext=FID_EXT) + if os.path.exists(noGapsFid): paramsNoGapModel2Point = { - 'inputFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_noGaps", - extension=".fid")), - 'outputFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_noGaps_fid", - extension=".txt")) + 'inputFile': noGapsFid, + 'outputFile': self.getExtraOutFile(tsId, suffix="noGaps_fid", ext=TXT_EXT) } argsNoGapModel2Point = "-InputFile %(inputFile)s " \ "-OutputFile %(outputFile)s" @@ -449,32 +495,15 @@ def translateFiducialPointModelStep(self, tsObjId): Plugin.runImod(self, 'model2point', argsNoGapModel2Point % paramsNoGapModel2Point) @ProtImodBase.tryExceptDecorator - def computeOutputStackStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - - firstItem = ts.getFirstItem() + def computeOutputStackStep(self, tsId): + ts = self.tsDict[tsId] # Check that previous steps have been completed satisfactorily - tmpFileName = os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".xf")) - if os.path.exists(tmpFileName) and os.stat(tmpFileName).st_size != 0: - tltFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_interpolated", extension=".tlt") - ) + transformationMatricesFilePath = self.getExtraOutFile(tsId, suffix="fid", ext=XF_EXT) + if os.path.exists(transformationMatricesFilePath) and os.stat(transformationMatricesFilePath).st_size != 0: + tltFilePath = self.getExtraOutFile(tsId, suffix="interpolated", ext=TLT_EXT) tltList = utils.formatAngleList(tltFilePath) - - transformationMatricesFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_fid", extension=".xf") - ) - newTransformationMatricesList = utils.formatTransformationMatrix(transformationMatricesFilePath) - output = self.getOutputSetOfTiltSeries(self.inputSetOfTiltSeries.get()) newTs = TiltSeries(tsId=tsId) newTs.copyInfo(ts) @@ -514,30 +543,23 @@ def computeOutputStackStep(self, tsObjId): else: raise FileNotFoundError( "Error (computeOutputStackStep): \n Imod output file " - "%s does not exist or it is empty" % tmpFileName) + "%s does not exist or it is empty" % transformationMatricesFilePath) @ProtImodBase.tryExceptDecorator - def computeOutputInterpolatedStackStep(self, tsObjId): - tsIn = self.inputSetOfTiltSeries.get()[tsObjId] + def computeOutputInterpolatedStackStep(self, tsId): + tsIn = self.tsDict[tsId] tsId = tsIn.getTsId() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - firstItem = tsIn.getFirstItem() # Check that previous steps have been completed satisfactorily - tmpFileName = os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".xf")) + tmpFileName = self.getExtraOutFile(tsId, suffix="fid", ext=XF_EXT) if os.path.exists(tmpFileName) and os.stat(tmpFileName).st_size != 0: output = self.getOutputInterpolatedSetOfTiltSeries(self.inputSetOfTiltSeries.get()) paramsAlignment = { - 'input': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'output': os.path.join(extraPrefix, firstItem.parseFileName()), - 'xform': os.path.join(extraPrefix, firstItem.parseFileName(suffix="_fid", - extension=".xf")), + 'input': self.getTmpOutFile(tsId), + 'output': self.getExtraOutFile(tsId), + 'xform': tmpFileName, 'bin': self.binning.get(), 'imagebinned': 1.0} @@ -570,11 +592,7 @@ def computeOutputInterpolatedStackStep(self, tsObjId): newTs.setInterpolated(True) output.append(newTs) - tltFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_interpolated", extension=".tlt") - ) - + tltFilePath = self.getExtraOutFile(tsId, suffix="interpolated", ext=TLT_EXT) tltList = utils.formatAngleList(tltFilePath) if self.binning > 1: @@ -584,7 +602,7 @@ def computeOutputInterpolatedStackStep(self, tsObjId): newTi = TiltImage() newTi.copyInfo(tiltImage, copyId=True, copyTM=False) newTi.setAcquisition(tiltImage.getAcquisition()) - newTi.setLocation(index + 1, os.path.join(extraPrefix, tiltImage.parseFileName())) + newTi.setLocation(index + 1, self.getExtraOutFile(tsId)) newTi.setTiltAngle(float(tltList[index])) if self.binning > 1: newTi.setSamplingRate(tiltImage.getSamplingRate() * self.binning.get()) @@ -604,22 +622,15 @@ def computeOutputInterpolatedStackStep(self, tsObjId): "Imod output file %s does not exist or it is empty" % tmpFileName) @ProtImodBase.tryExceptDecorator - def eraseGoldBeadsStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] + def eraseGoldBeadsStep(self, tsId): + ts = self.tsDict[tsId] tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - - firstItem = ts.getFirstItem() - # Erase gold beads on aligned stack paramsCcderaser = { - 'inputFile': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'outputFile': os.path.join(extraPrefix, firstItem.parseFileName()), - 'modelFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_noGaps", - extension=".fid")), + 'inputFile': self.getTmpOutFile(tsId), + 'outputFile': self.getExtraOutFile(tsId), + 'modelFile': self.getExtraOutFile(tsId, suffix="noGaps", ext=FID_EXT), 'betterRadius': self.betterRadius.get() / 2, 'polynomialOrder': 0, 'circleObjects': "/" @@ -639,44 +650,19 @@ def eraseGoldBeadsStep(self, tsObjId): Plugin.runImod(self, 'ccderaser', argsCcderaser % paramsCcderaser) @ProtImodBase.tryExceptDecorator - def computeOutputModelsStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] + def computeOutputModelsStep(self, tsId): + ts = self.tsDict[tsId] tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - - firstItem = ts.getFirstItem() # Create the output set of landmark models with no gaps - if os.path.exists( - os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(suffix="_noGaps_fid", - extension=".txt"))): - + fiducialNoGapFilePath = self.getExtraOutFile(tsId, suffix="noGaps_fid", ext=TXT_EXT) + if os.path.exists(fiducialNoGapFilePath): output = self.getOutputFiducialModelNoGaps() - output.setSetOfTiltSeries(self.inputSetOfTiltSeries.get()) - - fiducialNoGapFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_noGaps_fid", extension=".txt") - ) - fiducialNoGapList = utils.formatFiducialList(fiducialNoGapFilePath) - - fiducialModelNoGapPath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_noGaps", extension=".fid") - ) - - landmarkModelNoGapsFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_noGaps", extension=".sfid") - ) - - landmarkModelNoGapsResidPath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_resid", extension=".txt") - ) + fiducialModelNoGapPath = self.getExtraOutFile(tsId, suffix="noGaps", ext=FID_EXT) + landmarkModelNoGapsFilePath = self.getExtraOutFile(tsId, suffix="noGaps", ext=SFID_EXT) + landmarkModelNoGapsResidPath = self.getExtraOutFile(tsId, suffix="resid", ext=TXT_EXT) fiducialNoGapsResidList = utils.formatFiducialResidList(landmarkModelNoGapsResidPath) @@ -689,10 +675,12 @@ def computeOutputModelsStep(self, tsObjId): prevTiltIm = 0 chainId = 0 indexFake = 0 + firstExec = True for fiducial in fiducialNoGapList: - if int(float(fiducial[2])) <= prevTiltIm: + if (int(float(fiducial[2])) <= prevTiltIm) or firstExec: chainId += 1 + firstExec = False prevTiltIm = int(float(fiducial[2])) if indexFake < len(fiducialNoGapsResidList) and fiducial[2] == fiducialNoGapsResidList[indexFake][2]: @@ -709,17 +697,15 @@ def computeOutputModelsStep(self, tsObjId): yCoor=fiducial[1], tiltIm=fiducial[2] + 1, chainId=chainId, - xResid='0', - yResid='0') + xResid=float('nan'), + yResid=float('nan')) output.append(landmarkModelNoGaps) output.update(landmarkModelNoGaps) output.write() # Create the output set of 3D coordinates - coordFilePath = os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".xyz")) + coordFilePath = self.getExtraOutFile(tsId, suffix="fid", ext=XYZ_EXT) if os.path.exists(coordFilePath): @@ -753,6 +739,10 @@ def createOutputStep(self): self._store() + def createOutputFailedStep(self, tsId): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts) + # --------------------------- UTILS functions ----------------------------- def getRotationType(self): if self.rotationSolutionType.get() == 0: diff --git a/imod/protocols/protocol_fiducialModel.py b/imod/protocols/protocol_fiducialModel.py index 714f6ed4..749c3b69 100644 --- a/imod/protocols/protocol_fiducialModel.py +++ b/imod/protocols/protocol_fiducialModel.py @@ -33,13 +33,14 @@ import tomo.objects as tomoObj from .. import Plugin, utils -from .protocol_base import ProtImodBase +from .protocol_base import ProtImodBase, TLT_EXT, XF_EXT, FID_EXT, TXT_EXT, SEED_EXT, SFID_EXT class ProtImodFiducialModel(ProtImodBase): """ Construction of a fiducial model and alignment of tilt-series based on the IMOD procedure. + More info: https://bio3d.colorado.edu/imod/doc/man/autofidseed.html https://bio3d.colorado.edu/imod/doc/man/beadtrack.html @@ -49,65 +50,125 @@ class ProtImodFiducialModel(ProtImodBase): _label = 'Generate fiducial model' _devStatus = BETA + FIDUCIAL_MODEL = 0 + PATCH_TRACKING = 1 + # -------------------------- DEFINE param functions ----------------------- def _defineParams(self, form): form.addSection('Input') + form.addParam('typeOfModel', + params.EnumParam, + choices=["Make seed and Track", "Patch Tracking"], + default=0, + important=True, + label='Model generation') + form.addParam('inputSetOfTiltSeries', params.PointerParam, pointerClass='SetOfTiltSeries', important=True, - label='Input set of tilt-series') - - form.addParam('fiducialDiameter', - params.FloatParam, - label='Fiducial diameter (nm)', - default='10', - important=True, - help="Fiducials diameter to be tracked for alignment.") - - form.addParam('twoSurfaces', - params.EnumParam, - choices=['Yes', 'No'], - default=1, - label='Find beads on two surfaces?', - display=params.EnumParam.DISPLAY_HLIST, - help="Track fiducials differentiating in which side " - "of the sample are located.") - - form.addParam('numberFiducial', - params.IntParam, - label='Number of fiducials', - default='25', - expertLevel=params.LEVEL_ADVANCED, - help="Number of fiducials to be tracked for alignment.") - - form.addParam('doTrackWithModel', params.BooleanParam, - default=True, - label="Track with fiducial model as seed", - help="Turn the tracked model into new seed and " - "repeat tracking.") - - form.addParam('shiftsNearZeroFraction', - params.FloatParam, - label='Shifts near zero fraction', - default='0.2', - expertLevel=params.LEVEL_ADVANCED, - help="Fraction of the tracking box size above which to " - "supply shifts near zero tilt to Beadtrack. The " - "dominant net shifts in the bead positions between " - "views are found as described above, and if one of " - "the shifts is larger than this fraction of the " - "-BoxSizeXandY entry to Beadtrack, then the shifts " - "are provided when running Beadtrack on the initial " - "seed models. Also, a command file will be written " - "with modified parameters, named as the root name " - "of the input command file followed by '_adjusted' " - "and its extension. Enter 0 or a large value to " - "disable this analysis.") + label='Tilt Series') + + self._patchTrackingForm(form, 'typeOfModel==%d' % self.PATCH_TRACKING) + self._fiducialSeedForm(form, 'typeOfModel==%d' % self.FIDUCIAL_MODEL) + + def _patchTrackingForm(self, form, condition, levelType=params.LEVEL_NORMAL): + patchtrack = form.addGroup('Patch Tracking', expertLevel=levelType, condition=condition) + + patchtrack.addParam('sizeOfPatches', params.NumericListParam, + label='Size of the patches (X,Y)', + default='100 100', + expertLevel=levelType, + help="Size of the patches to track by correlation. In imod documentation " + "(tiltxcorr: SizeOfPatchesXandY)") + + patchtrack.addParam('patchLayout', + params.EnumParam, + choices=['Fractional overlap of patches', + 'Number of patches'], + default=0, + label='Patch layout', + display=params.EnumParam.DISPLAY_HLIST, + help='To be added') + + patchtrack.addParam('overlapPatches', + params.NumericListParam, + default='0.33 0.33', + condition='patchLayout==0', + label='Fractional overlap of the patches (X,Y)', + help="Fractional overlap in X and Y to track by correlation. In imod documentation" + "(tiltxcorr: NumberOfPatchesXandY)") + + patchtrack.addParam('numberOfPatches', + params.NumericListParam, + condition='patchLayout==1', + label='Number of patches (X,Y)', + help="Number of patches in X and Y of the patches. In imod documentation" + "(tiltxcorr: OverlapOfPatchesXandY)") + + patchtrack.addParam('iterationsSubpixel', + params.IntParam, + default=1, + label='Iterations to increase subpixel accuracy', + help="Number of iteration of each correlation to reduce interpolation of the peak position" + "In imod documentation: (tiltxcorr: IterateCorrelations)") + + self.trimimgForm(patchtrack, pxTrimCondition='True', correlationCondition='True', + levelType=params.LEVEL_ADVANCED) + self.filteringParametersForm(form, condition=condition, levelType=params.LEVEL_ADVANCED) + + def _fiducialSeedForm(self, form, condition, levelType=params.LEVEL_NORMAL): + seedModel = form.addGroup('"Make seed and Track', expertLevel=levelType, condition=condition) + seedModel.addParam('fiducialDiameter', + params.FloatParam, + label='Fiducial diameter (nm)', + default='10', + important=True, + help="Fiducials diameter to be tracked for alignment.") + + seedModel.addParam('twoSurfaces', + params.EnumParam, + choices=['Yes', 'No'], + default=1, + label='Find beads on two surfaces?', + display=params.EnumParam.DISPLAY_HLIST, + help="Track fiducials differentiating in which side " + "of the sample are located.") + + seedModel.addParam('numberFiducial', + params.IntParam, + label='Number of fiducials', + default='25', + expertLevel=params.LEVEL_ADVANCED, + help="Number of fiducials to be tracked for alignment.") + + seedModel.addParam('doTrackWithModel', params.BooleanParam, + default=True, + label="Track with fiducial model as seed", + help="Turn the tracked model into new seed and " + "repeat tracking.") + + seedModel.addParam('shiftsNearZeroFraction', + params.FloatParam, + label='Shifts near zero fraction', + default='0.2', + expertLevel=params.LEVEL_ADVANCED, + help="Fraction of the tracking box size above which to " + "supply shifts near zero tilt to Beadtrack. The " + "dominant net shifts in the bead positions between " + "views are found as described above, and if one of " + "the shifts is larger than this fraction of the " + "-BoxSizeXandY entry to Beadtrack, then the shifts " + "are provided when running Beadtrack on the initial " + "seed models. Also, a command file will be written " + "with modified parameters, named as the root name " + "of the input command file followed by '_adjusted' " + "and its extension. Enter 0 or a large value to " + "disable this analysis.") groupGlobalVariables = form.addGroup('Filter variables', - expertLevel=params.LEVEL_ADVANCED) + expertLevel=params.LEVEL_ADVANCED, condition=condition) groupGlobalVariables.addParam('refineSobelFilter', params.EnumParam, @@ -136,43 +197,39 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self._failedTs = [] - - for ts in self.inputSetOfTiltSeries.get(): - tsObjId = ts.getObjId() - self._insertFunctionStep(self.convertInputStep, tsObjId) - self._insertFunctionStep(self.generateTrackComStep, tsObjId) - self._insertFunctionStep(self.generateFiducialSeedStep, tsObjId) - self._insertFunctionStep(self.generateFiducialModelStep, tsObjId) - self._insertFunctionStep(self.translateFiducialPointModelStep, tsObjId) - self._insertFunctionStep(self.computeOutputModelsStep, tsObjId) - self._insertFunctionStep(self.createOutputFailedSet, tsObjId) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + if self.typeOfModel == self.FIDUCIAL_MODEL: + self._insertFunctionStep(self.generateTrackComStep, tsId) + self._insertFunctionStep(self.generateFiducialSeedStep, tsId) + self._insertFunctionStep(self.generateFiducialModelStep, tsId) + else: + self._insertFunctionStep(self.xcorrStep, tsId) + self._insertFunctionStep(self.chopcontsStep, tsId) + self._insertFunctionStep(self.translateFiducialPointModelStep, tsId) + self._insertFunctionStep(self.computeOutputModelsStep, tsId) + self._insertFunctionStep(self.createOutputFailedSetStep, tsId) self._insertFunctionStep(self.createOutputStep) # --------------------------- STEPS functions ----------------------------- - def generateTrackComStep(self, tsObjId): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - - firstItem = ts.getFirstItem() + def _initialize(self): + self._failedTs = [] + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + def generateTrackComStep(self, tsId): + ts = self.tsDict[tsId] fiducialDiameterPixel = self.fiducialDiameter.get() / (self.inputSetOfTiltSeries.get().getSamplingRate() / 10) boxSizeXandY = max(3.3 * fiducialDiameterPixel + 2, 2 * fiducialDiameterPixel + 20, 32) boxSizeXandY = min(512, 2 * int(boxSizeXandY / 2)) scaling = fiducialDiameterPixel / 12.5 if fiducialDiameterPixel > 12.5 else 1 paramsDict = { - 'imageFile': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'inputSeedModel': os.path.join(extraPrefix, - firstItem.parseFileName(extension=".seed")), - 'outputModel': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps", - extension=".fid")), - 'tiltFile': os.path.join(tmpPrefix, firstItem.parseFileName(extension=".tlt")), + 'imageFile': self.getTmpOutFile(tsId), + 'inputSeedModel': self.getExtraOutFile(tsId, ext=SEED_EXT), + 'outputModel': self.getExtraOutFile(tsId, suffix="gaps", ext=FID_EXT), + 'tiltFile': self.getExtraOutFile(tsId, ext=TLT_EXT), 'rotationAngle': ts.getAcquisition().getTiltAxisAngle(), 'fiducialDiameter': fiducialDiameterPixel, 'samplingRate': self.inputSetOfTiltSeries.get().getSamplingRate() / 10, @@ -188,16 +245,9 @@ def generateTrackComStep(self, tsObjId): self.translateTrackCom(ts, paramsDict) @ProtImodBase.tryExceptDecorator - def generateFiducialSeedStep(self, tsObjId): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - + def generateFiducialSeedStep(self, tsId): paramsAutofidseed = { - 'trackCommandFile': os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(suffix="_track", - extension=".com")), + 'trackCommandFile': self.getExtraOutFile(tsId, suffix="track", ext="com"), 'minSpacing': 0.85, 'peakStorageFraction': 1.0, 'targetNumberOfBeads': self.numberFiducial.get(), @@ -222,13 +272,8 @@ def generateFiducialSeedStep(self, tsObjId): path.moveFile("autofidseed.info", self._getExtraPath(tsId)) @ProtImodBase.tryExceptDecorator - def generateFiducialModelStep(self, tsObjId): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) - + def generateFiducialModelStep(self, tsId): + ts = self.tsDict[tsId] firstItem = ts.getFirstItem() fiducialDiameterPixel = self.fiducialDiameter.get() / (self.inputSetOfTiltSeries.get().getSamplingRate() / 10) @@ -237,15 +282,11 @@ def generateFiducialModelStep(self, tsObjId): scaling = fiducialDiameterPixel / 12.5 if fiducialDiameterPixel > 12.5 else 1 paramsBeadtrack = { - 'inputSeedModel': os.path.join(extraPrefix, - firstItem.parseFileName(extension=".seed")), - 'outputModel': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps", - extension=".fid")), - 'imageFile': os.path.join(tmpPrefix, firstItem.parseFileName()), + 'inputSeedModel': self.getExtraOutFile(tsId, ext=SEED_EXT), + 'outputModel': self.getExtraOutFile(tsId, suffix="gaps", ext=FID_EXT), + 'imageFile': self.getTmpOutFile(tsId), 'imagesAreBinned': 1, - 'tiltFile': os.path.join(tmpPrefix, - firstItem.parseFileName(extension=".tlt")), + 'tiltFile': self.getExtraOutFile(tsId, ext=TLT_EXT), 'tiltDefaultGrouping': 7, 'magDefaultGrouping': 5, 'rotDefaultGrouping': 1, @@ -314,8 +355,7 @@ def generateFiducialModelStep(self, tsObjId): argsBeadtrack += f"-SkipViews {','.join(excludedViews)} " if firstItem.hasTransform(): - XfFileName = os.path.join(tmpPrefix, - firstItem.parseFileName(extension=".xf")) + XfFileName = self.getExtraOutFile(tsId, ext=XF_EXT) argsBeadtrack += f"-prexf {XfFileName} " Plugin.runImod(self, 'beadtrack', argsBeadtrack % paramsBeadtrack) @@ -323,74 +363,39 @@ def generateFiducialModelStep(self, tsObjId): if self.doTrackWithModel: # repeat tracking with the current model as seed path.copyFile(paramsBeadtrack['inputSeedModel'], - os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_orig", - extension=".seed"))) + self.getExtraOutFile(tsId, suffix="orig", ext=SEED_EXT)) path.moveFile(paramsBeadtrack['outputModel'], paramsBeadtrack['inputSeedModel']) Plugin.runImod(self, 'beadtrack', argsBeadtrack % paramsBeadtrack) - def translateFiducialPointModelStep(self, tsObjId): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - - firstItem = ts.getFirstItem() - + def translateFiducialPointModelStep(self, tsId): # Check that previous steps have been completed satisfactorily - if os.path.exists(os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps", - extension=".fid"))): + gapsFidFile = self.getExtraOutFile(tsId, suffix='gaps', ext=FID_EXT) + if os.path.exists(gapsFidFile): paramsGapModel2Point = { - 'inputFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps", - extension=".fid")), - 'outputFile': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps_fid", - extension=".txt")) + 'inputFile': gapsFidFile, + 'outputFile': self.getExtraOutFile(tsId, suffix="gaps_fid", ext=TXT_EXT) } argsGapModel2Point = "-InputFile %(inputFile)s " \ "-OutputFile %(outputFile)s" Plugin.runImod(self, 'model2point', argsGapModel2Point % paramsGapModel2Point) - def computeOutputModelsStep(self, tsObjId): - ts = self._getTiltSeries(tsObjId) - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - - firstItem = ts.getFirstItem() + def computeOutputModelsStep(self, tsId): + ts = self.tsDict[tsId] # Create the output set of landmark models with gaps # Check that previous steps have been completed satisfactorily - if os.path.exists( - os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_gaps", - extension=".fid"))): - + fiducialModelGapPath = self.getExtraOutFile(tsId, suffix='gaps', ext=FID_EXT) + if os.path.exists(fiducialModelGapPath): output = self.getOutputFiducialModelGaps() - - landmarkModelGapsFilePath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_gaps", extension=".sfid") - ) - - fiducialModelGapPath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_gaps", extension=".fid") - ) - - fiducialModelGapTxtPath = os.path.join( - extraPrefix, - firstItem.parseFileName(suffix="_gaps_fid", extension=".txt") - ) + landmarkModelGapsFilePath = self.getExtraOutFile(tsId, suffix='gaps', ext=SFID_EXT) + fiducialModelGapTxtPath = self.getExtraOutFile(tsId, suffix="gaps_fid", ext=TXT_EXT) fiducialGapList = utils.formatFiducialList(fiducialModelGapTxtPath) fiducialDiameterPixel = self.fiducialDiameter.get() / ( - self.inputSetOfTiltSeries.get().getSamplingRate() / 10) + self.inputSetOfTiltSeries.get().getSamplingRate() / 10) landmarkModelGaps = tomoObj.LandmarkModel(tsId=tsId, tiltSeriesPointer=ts, @@ -421,6 +426,102 @@ def computeOutputModelsStep(self, tsObjId): output.update(landmarkModelGaps) output.write() + def xcorrStep(self, tsId): + """ + Imod uses the next command line for the xcorr alignment + $tiltxcorr -StandardInput + InputFile cryo_preali.mrc + OutputFile cryo_pt.fid + RotationAngle -12.6 + TiltFile cryo.rawtlt + FilterRadius2 0.125 + FilterSigma1 0.03 + FilterSigma2 0.03 + BordersInXandY 102,102 + IterateCorrelations 1 + SizeOfPatchesXandY 680,680 + OverlapOfPatchesXandY 0.33,0.33 + PrealignmentTransformFile cryo.prexg + ImagesAreBinned 1 + """ + ts = self.tsDict[tsId] + angleFilePath = self.getExtraOutFile(tsId, ext=TLT_EXT) + xfFile = self.getExtraOutFile(tsId, ext=XF_EXT) + ts.writeXfFile(xfFile) + + borders = self.pxTrim.getListFromValues() + sizePatches = self.sizeOfPatches.getListFromValues() + + BordersInXandY = '%d,%d' % (borders[0], borders[1]) + SizeOfPatchesXandY = '%d,%d' % (sizePatches[0], sizePatches[1]) + + paramsTiltXCorr = { + 'inputFile': self.getTmpOutFile(tsId), + 'outputFile': self.getExtraOutFile(tsId, suffix="pt", ext=FID_EXT), + 'RotationAngle': ts.getAcquisition().getTiltAxisAngle(), + 'TiltFile': angleFilePath, + 'FilterRadius2': self.filterRadius2.get(), + 'FilterSigma1': self.filterSigma1.get(), + 'FilterSigma2': self.filterSigma2.get(), + 'BordersInXandY': BordersInXandY, + 'IterateCorrelations': self.iterationsSubpixel.get(), + 'SizeOfPatchesXandY': SizeOfPatchesXandY, + 'PrealignmentTransformFile': xfFile, + + 'ImagesAreBinned': 1, + } + argsTiltXCorr = " " \ + "-InputFile %(inputFile)s " \ + "-OutputFile %(outputFile)s " \ + "-RotationAngle %(RotationAngle)s " \ + "-TiltFile %(TiltFile)s " \ + "-FilterRadius2 %(FilterRadius2)s " \ + "-FilterSigma1 %(FilterSigma1)s " \ + "-FilterSigma2 %(FilterSigma2)s " \ + "-BordersInXandY %(BordersInXandY)s " \ + "-IterateCorrelations %(IterateCorrelations)s " \ + "-SizeOfPatchesXandY %(SizeOfPatchesXandY)s " \ + "-PrealignmentTransformFile %(PrealignmentTransformFile)s " \ + "-ImagesAreBinned %(ImagesAreBinned)s " + + if self.patchLayout.get() == 0: + patchesXY = self.overlapPatches.getListFromValues(caster=float) + OverlapOfPatchesXandY = '%f,%f' % (patchesXY[0], patchesXY[1]) + argsTiltXCorr += ' -OverlapOfPatchesXandY %s ' % OverlapOfPatchesXandY + else: + numberPatchesXY = self.numberOfPatches.getListFromValues() + argsTiltXCorr += ' -NumberOfPatchesXandY %d,%d ' % (numberPatchesXY[0], numberPatchesXY[1]) + + Plugin.runImod(self, 'tiltxcorr', argsTiltXCorr % paramsTiltXCorr) + + def chopcontsStep(self, tsId): + """ + $imodchopconts -StandardInput + InputModel cryo_pt.fid + OutputModel cryo.fid + MinimumOverlap 4 + AssignSurfaces 1 + """ + MinimumOverlap = 4 + AssignSurfaces = 1 + LengthOfPieces = -1 + + paramschopconts = { + 'inputFile': self.getExtraOutFile(tsId, suffix="pt", ext=FID_EXT), + 'outputFile': self.getExtraOutFile(tsId, suffix="gaps", ext=FID_EXT), + 'MinimumOverlap': MinimumOverlap, + 'AssignSurfaces': AssignSurfaces, + 'LengthOfPieces': LengthOfPieces + } + argschopconts = " " \ + "-InputModel %(inputFile)s " \ + "-OutputModel %(outputFile)s " \ + "-MinimumOverlap %(MinimumOverlap)s " \ + "-AssignSurfaces %(AssignSurfaces)s " \ + "-LengthOfPieces %(LengthOfPieces)s " + + Plugin.runImod(self, 'imodchopconts', argschopconts % paramschopconts) + def createOutputStep(self): if self.FiducialModelGaps: self.FiducialModelGaps.setStreamState(Set.STREAM_CLOSED) @@ -429,11 +530,14 @@ def createOutputStep(self): self._store() + def createOutputFailedSetStep(self, tsId): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts) + # --------------------------- UTILS functions ----------------------------- def translateTrackCom(self, ts, paramsDict): tsId = ts.getTsId() extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) firstItem = ts.getFirstItem() trackFilePath = os.path.join(extraPrefix, @@ -515,8 +619,7 @@ def translateTrackCom(self, ts, paramsDict): template += f"SkipViews {','.join(excludedViews)}" if firstItem.hasTransform(): - XfFileName = os.path.join(tmpPrefix, - firstItem.parseFileName(extension=".xf")) + XfFileName = self.getExtraOutFile(tsId, ext=XF_EXT) template += f"PrealignTransformFile {XfFileName}" with open(trackFilePath, 'w') as f: diff --git a/imod/protocols/protocol_goldBeadPicker3d.py b/imod/protocols/protocol_goldBeadPicker3d.py index 8abc6119..64ff9556 100644 --- a/imod/protocols/protocol_goldBeadPicker3d.py +++ b/imod/protocols/protocol_goldBeadPicker3d.py @@ -24,17 +24,14 @@ # * # ***************************************************************************** -import os - from pyworkflow import BETA import pyworkflow.protocol.params as params -from pyworkflow.utils import path, removeBaseExt from pyworkflow.object import Set import tomo.objects as tomoObj import tomo.constants as constants from .. import Plugin, utils -from .protocol_base import ProtImodBase +from .protocol_base import ProtImodBase, XYZ_EXT, MOD_EXT class ProtImodGoldBeadPicker3d(ProtImodBase): @@ -107,39 +104,39 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self.defineExecutionPararell() - allOutputId = [] - - for ts in self.inputSetOfTomograms.get(): + self._initialize() + for tsId in self.tomoDict.keys(): pickId = self._insertFunctionStep(self.pickGoldBeadsStep, - ts.getObjId(), + tsId, prerequisites=[]) convertId = self._insertFunctionStep(self.convertModelToCoordinatesStep, - ts.getObjId(), + tsId, prerequisites=[pickId]) outputID = self._insertFunctionStep(self.createOutputStep, - ts.getObjId(), + tsId, prerequisites=[convertId]) allOutputId.append(outputID) - self._insertFunctionStep('closeOutputSetStep', + self._insertFunctionStep(self.closeOutputSetStep, prerequisites=allOutputId) # --------------------------- STEPS functions ----------------------------- - def pickGoldBeadsStep(self, tsObjId): - tomo = self.inputSetOfTomograms.get()[tsObjId] - fileName = removeBaseExt(tomo.getFileName()) - extraPrefix = self._getExtraPath(fileName) - path.makePath(extraPrefix) + def _initialize(self): + self.defineExecutionPararell() + self.tomoDict = {tomo.getTsId(): tomo.clone() for tomo in self.inputSetOfTomograms.get()} + + def pickGoldBeadsStep(self, tsId): + self.genTsPaths(tsId) + tomo = self.tomoDict[tsId] """ Run findbeads3d IMOD program """ paramsFindbeads3d = { 'inputFile': tomo.getFileName(), - 'outputFile': os.path.join(extraPrefix, "%s.mod" % fileName), + 'outputFile': self.getExtraOutFile(tsId, ext=MOD_EXT), 'beadSize': self.beadDiameter.get(), 'minRelativeStrength': self.minRelativeStrength.get(), 'minSpacing': self.minSpacing.get(), @@ -157,17 +154,11 @@ def pickGoldBeadsStep(self, tsObjId): Plugin.runImod(self, 'findbeads3d', argsFindbeads3d % paramsFindbeads3d) - def convertModelToCoordinatesStep(self, tsObjId): - tomo = self.inputSetOfTomograms.get()[tsObjId] - location = tomo.getFileName() - fileName, _ = os.path.splitext(location) - - extraPrefix = self._getExtraPath(os.path.basename(fileName)) - + def convertModelToCoordinatesStep(self, tsId): """ Run model2point IMOD program """ paramsModel2Point = { - 'inputFile': os.path.join(extraPrefix, "%s.mod" % os.path.basename(fileName)), - 'outputFile': os.path.join(extraPrefix, "%s.xyz" % os.path.basename(fileName)), + 'inputFile': self.getExtraOutFile(tsId, ext=MOD_EXT), + 'outputFile': self.getExtraOutFile(tsId, ext=XYZ_EXT), } argsModel2Point = "-InputFile %(inputFile)s " \ @@ -175,20 +166,14 @@ def convertModelToCoordinatesStep(self, tsObjId): Plugin.runImod(self, 'model2point', argsModel2Point % paramsModel2Point) - def createOutputStep(self, tsObjId): - tomo = self.inputSetOfTomograms.get()[tsObjId] - location = tomo.getFileName() - fileName, _ = os.path.splitext(location) - - extraPrefix = self._getExtraPath(os.path.basename(fileName)) + def createOutputStep(self, tsId): + tomo = self.tomoDict[tsId] """ Create the output set of coordinates 3D from gold beads detected """ output = self.getOutputSetOfCoordinates3Ds(self.inputSetOfTomograms.get(), self.inputSetOfTomograms.get()) - coordFilePath = os.path.join(extraPrefix, - "%s.xyz" % os.path.basename(fileName)) - + coordFilePath = self.getExtraOutFile(tsId, ext=XYZ_EXT) coordList = utils.formatGoldBead3DCoordinatesList(coordFilePath) with self._lock: @@ -199,7 +184,7 @@ def createOutputStep(self, tsObjId): newCoord3D.setY(element[1], constants.BOTTOM_LEFT_CORNER) newCoord3D.setZ(element[2], constants.BOTTOM_LEFT_CORNER) - newCoord3D.setVolId(tsObjId) + # newCoord3D.setVolId(tsObjId) output.append(newCoord3D) output.update(newCoord3D) diff --git a/imod/protocols/protocol_importSetOfTM.py b/imod/protocols/protocol_importSetOfTM.py index d35f16d6..48fc2256 100644 --- a/imod/protocols/protocol_importSetOfTM.py +++ b/imod/protocols/protocol_importSetOfTM.py @@ -83,25 +83,23 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self.matchBinningFactor = self.binningTM.get() / self.binningTS.get() - - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.generateTransformFileStep, - ts.getObjId()) - self._insertFunctionStep(self.assignTransformationMatricesStep, - ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.generateTransformFileStep, tsId) + self._insertFunctionStep(self.assignTransformationMatricesStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def generateTransformFileStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() + def _initialize(self): + self.matchBinningFactor = self.binningTM.get() / self.binningTS.get() + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + def generateTransformFileStep(self, tsId): + self.genTsPaths(tsId) + ts = self.tsDict[tsId] tsFileName = ts.getFirstItem().parseFileName(extension='') - extraPrefix = self._getExtraPath(tsId) - path.makePath(extraPrefix) outputTransformFile = os.path.join(extraPrefix, ts.getFirstItem().parseFileName(extension=".xf")) @@ -147,12 +145,9 @@ def generateTransformFileStep(self, tsObjId): else: path.createLink(tmFilePath, outputTransformFile) - def assignTransformationMatricesStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - + def assignTransformationMatricesStep(self, tsId): + ts = self.tsDict[tsId] extraPrefix = self._getExtraPath(tsId) - outputTransformFile = os.path.join(extraPrefix, ts.getFirstItem().parseFileName(extension=".xf")) diff --git a/imod/protocols/protocol_tomoNormalization.py b/imod/protocols/protocol_tomoPreprocess.py similarity index 69% rename from imod/protocols/protocol_tomoNormalization.py rename to imod/protocols/protocol_tomoPreprocess.py index f1b0985f..737a4451 100644 --- a/imod/protocols/protocol_tomoNormalization.py +++ b/imod/protocols/protocol_tomoPreprocess.py @@ -23,25 +23,39 @@ # * e-mail address 'scipion@cnb.csic.es' # * # ***************************************************************************** - -import os - from pyworkflow import BETA from pyworkflow.object import Set import pyworkflow.protocol.params as params import pyworkflow.utils.path as pwpath from tomo.objects import Tomogram, SetOfTomograms - from .. import Plugin -from .protocol_base import ProtImodBase, EXT_MRC_ODD_NAME, EXT_MRC_EVEN_NAME, OUTPUT_TOMOGRAMS_NAME +from .protocol_base import ProtImodBase, OUTPUT_TOMOGRAMS_NAME, MRC_EXT, ODD, EVEN class ProtImodTomoNormalization(ProtImodBase): """ Normalize input tomogram and change its storing formatting. + More info: https://bio3D.colorado.edu/imod/doc/newstack.html - https://bio3D.colorado.edu/imod/doc/binvol.html + https://bio3d.colorado.edu/imod/doc/man/binvol.html + + IMOD tilt series preprocess makes use of the Newstack and + binvol commands. In particular, three functionalities are possible:\n + + _1 Binning_: Binvol will bin down a volume in all three dimensions, + with the binning done isotropically. Binning means summing (actually + averaging) all of the values in a block of voxels (e.g., 2x2x2 + or 1x1x3) in the input volume to create one voxel in the output volume. + The output file will have appropriately larger pixel spacings + in its header.\n + _2 Normalization_: This protocol allows to scale the gray values + of the tomograms, also called normalization, to a common range or + mean of density. The most used normalization consists in zero + mean and standard deviation one.\n + + _3 storage format_: IMOD is able to modify the number of bit of + the stored data in order to reduce the disc occupancy. """ _label = 'Tomo preprocess' @@ -55,50 +69,117 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTomograms', important=True, - label='Input set of tomograms') + label='Input set of tomograms', + help='Introduce the set tomograms to be normalized') form.addParam('binning', params.IntParam, default=1, label='Binning', important=True, - help='Binning to be applied to the normalized tomograms ' - 'in IMOD convention. Volumes will be binned by the ' - 'given factor. Must be an integer bigger than 1') + help='Binning is an scaling factor for the output tomograms. ' + 'Must be an integer greater than 1. IMOD uses ordinary' + 'binning to reduce tomograms in size by the given factor. ' + 'The value of a binned pixel is the average of pixel ' + 'values in each block of pixels being binned. Binning ' + 'is applied before all') form.addParam('floatDensities', params.EnumParam, - choices=['default', '1', '2', '3', '4'], + choices=['No adjust', + 'range between min and max', + 'scaled to common mean and standard deviation', + 'shifted to a common mean without scaling', + 'shifted to mean and rescaled to a min and max'], default=0, label='Adjust densities mode', - display=params.EnumParam.DISPLAY_HLIST, + display=params.EnumParam.DISPLAY_COMBO, help='Adjust densities of sections individually:\n' '-Default: no adjustment performed\n' - '-Mode 1: sections fill the data range\n' - '-Mode 2: sections scaled to common mean and standard deviation.\n' - '-Mode 3: sections shifted to a common mean without scaling\n' - '-Mode 4: sections shifted to a common mean and then ' - 'rescale the resulting minimum and maximum densities ' - 'to the Min and Max values specified') + + '-Range between min and max: This option will scale the gray values' + 'to be in a range given by a minimum and a maximum values.' + 'This is the mode 1 in newstack flag -floatDensities.\n' + + '-Scaled to common mean and standard deviation: This is the most ' + 'common normalization procedure. The new tomogramss will have' + 'a mean and a standard deviation introduced by the user. Generaly,' + 'a zero mean and a standard deviation one is a good choice.' + 'This is the mode 2 in newstack flag -floatDensities.\n' + + '-Shifted to a common mean without scaling: This option only' + 'add an offset to the gray values of the tomograms. The offset will' + 'be calculated such as the new tomograms will present a mean gray value' + 'introduced by the user. This is the mode 3 in newstack flag ' + 'floatDensities.\n' + + '-shifted to mean and rescaled to a min and max: In this case, an ' + 'offset is added to the tomograms in order to achieve a mean gray value' + ' then they are rescale the resulting minimum and maximum densities ' + 'to the Min and Max values specified. This is the mode 4 in newstack' + ' flag -floatDensities.\n') + + groupMeanSd = form.addGroup('Mean and SD', + condition='floatDensities==2', + help='Scale all images to the given mean ' + 'and standard deviation. This option ' + 'implies -float 2 and is incompatible ' + 'with all other scaling options. If no ' + 'values are set, mean=0 and SD=1 by default') + + groupMeanSd.addParam('meanSdToggle', + params.EnumParam, + choices=['Yes', 'No'], + default=1, + label='Set mean and SD?', + display=params.EnumParam.DISPLAY_HLIST, + help='Set mean and SD values') + + groupMeanSd.addParam('scaleMean', + params.FloatParam, + default=0, + label='Mean', + help='Mean value for the rescaling') + + groupMeanSd.addParam('scaleSd', + params.FloatParam, + default=1, + label='SD', + help='Standard deviation value for the rescaling') + + groupScale = form.addGroup('Scaling values', + condition='floatDensities==4') + + groupScale.addParam('scaleMax', + params.FloatParam, + default=255, + label='Max.', + help='Maximum value for the rescaling') + + groupScale.addParam('scaleMin', + params.FloatParam, + default=0, + label='Min.', + help='Minimum value for the rescaling') + form.addParam('modeToOutput', params.EnumParam, - choices=['default', '4-bit', 'byte', 'signed 16-bit', + expertLevel=params.LEVEL_ADVANCED, + choices=['default', '4-bit', '8-bit', 'signed 16-bit', 'unsigned 16-bit', '32-bit float'], default=0, label='Storage data type', - display=params.EnumParam.DISPLAY_HLIST, - help='Apply one density scaling to all sections to map ' - 'current min and max to the given Min and Max. The ' - 'storage mode of the output file. The default is ' - 'the mode of the first input file, except for a ' - '4-bit input file, where the default is to output ' - 'as bytes') + display=params.EnumParam.DISPLAY_COMBO, + help='The storage mode of the output file. The ' + 'default is the mode of the first input file, ' + 'except for a 4-bit input file, where the default ' + 'is to output as bytes') form.addParam('scaleRangeToggle', params.EnumParam, choices=['Yes', 'No'], - condition="floatDensities==0 or floatDensities==1 or floatDensities==3", - default=1, + condition="floatDensities==1 or floatDensities==3", + default=0, label='Set scaling range values?', display=params.EnumParam.DISPLAY_HLIST, help='This option will rescale the densities of all ' @@ -108,14 +189,14 @@ def _defineParams(self, form): form.addParam('scaleRangeMax', params.FloatParam, - condition="(floatDensities==0 or floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", + condition="(floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", default=255, label='Max.', help='Maximum value for the rescaling') form.addParam('scaleRangeMin', params.FloatParam, - condition="(floatDensities==0 or floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", + condition="(floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", default=0, label='Min.', help='Minimum value for the rescaling') @@ -123,13 +204,14 @@ def _defineParams(self, form): form.addParam('antialias', params.EnumParam, choices=['None', 'Blackman', 'Triangle', 'Mitchell', - 'Lanczos 2', 'Lanczos 3'], + 'Lanczos 2 lobes', 'Lanczos 3 lobes'], default=5, label='Antialias method:', - display=params.EnumParam.DISPLAY_HLIST, + expertLevel=params.LEVEL_ADVANCED, + display=params.EnumParam.DISPLAY_COMBO, help='Type of antialiasing filter to use when reducing images.\n' 'The available types of filters are:\n\n' - 'None\n' + 'None - Antialias will not be applied\n' 'Blackman - fast but not as good at antialiasing as slower filters\n' 'Triangle - fast but smooths more than Blackman\n' 'Mitchell - good at antialiasing, smooths a bit\n' @@ -141,49 +223,6 @@ def _defineParams(self, form): 'based on images of natural scenes where there are ' 'sharp edges.') - groupMeanSd = form.addGroup('Mean and SD', - condition='floatDensities==2', - help='Scale all images to the given mean ' - 'and standard deviation. This option ' - 'implies -float 2 and is incompatible ' - 'with all other scaling options. If no ' - 'values are set, mean=0 and SD=1 by default') - - groupMeanSd.addParam('meanSdToggle', - params.EnumParam, - choices=['Yes', 'No'], - default=1, - label='Set mean and SD?', - display=params.EnumParam.DISPLAY_HLIST, - help='Set mean and SD values') - - groupMeanSd.addParam('scaleMean', - params.FloatParam, - default=0, - label='Mean', - help='Mean value for the rescaling') - - groupMeanSd.addParam('scaleSd', - params.FloatParam, - default=1, - label='SD', - help='Standard deviation value for the rescaling') - - groupScale = form.addGroup('Scaling values', - condition='floatDensities==4') - - groupScale.addParam('scaleMax', - params.FloatParam, - default=255, - label='Max.', - help='Maximum value for the rescaling') - - groupScale.addParam('scaleMin', - params.FloatParam, - default=0, - label='Min.', - help='Minimum value for the rescaling') - form.addParam('processOddEven', params.BooleanParam, expertLevel=params.LEVEL_ADVANCED, @@ -194,23 +233,20 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - for tomo in self.inputSetOfTomograms.get(): - self._insertFunctionStep(self.generateOutputStackStep, - tomo.getObjId()) + self._initialize() + for tsId in self.tomoDict.keys(): + self._insertFunctionStep(self.generateOutputStackStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def generateOutputStackStep(self, tsObjId): - tomo = self.inputSetOfTomograms.get()[tsObjId] - location = tomo.getFileName() - fileName = os.path.splitext(location)[0] + def _initialize(self): + self.tomoDict = {tomo.getTsId(): tomo.clone() for tomo in self.inputSetOfTomograms.get()} - extraPrefix = self._getExtraPath(os.path.basename(fileName)) - outputFileNameBase = os.path.basename(pwpath.replaceExt(location, "mrc")) - outputFile = os.path.join(extraPrefix, outputFileNameBase) - tmpPrefix = self._getTmpPath(os.path.basename(fileName)) - pwpath.makePath(extraPrefix) - pwpath.makePath(tmpPrefix) + def generateOutputStackStep(self, tsId): + self.genTsPaths(tsId) + tomo = self.tomoDict[tsId] + location = tomo.getFileName() + outputFile = self.getExtraOutFile(tsId, ext=MRC_EXT) runNewstack = False @@ -253,11 +289,11 @@ def generateOutputStackStep(self, tsObjId): if self.applyToOddEven(tomo): oddFn, evenFn = tomo.getHalfMaps().split(',') paramsNewstack['input'] = oddFn - oddEvenOutput[0] = os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_ODD_NAME) + oddEvenOutput[0] = self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT) paramsNewstack['output'] = oddEvenOutput[0] Plugin.runImod(self, 'newstack', argsNewstack % paramsNewstack) paramsNewstack['input'] = evenFn - oddEvenOutput[1] = os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_EVEN_NAME) + oddEvenOutput[1] = self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT) paramsNewstack['output'] = oddEvenOutput[1] Plugin.runImod(self, 'newstack', argsNewstack % paramsNewstack) @@ -265,17 +301,16 @@ def generateOutputStackStep(self, tsObjId): if binning != 1: if runNewstack: - baseLoc = os.path.basename(outputFile) - tmpPath = os.path.join(tmpPrefix, baseLoc) - pwpath.moveFile(os.path.join(extraPrefix, baseLoc), tmpPath) + tmpPath = self.getTmpOutFile(tsId) + pwpath.moveFile(outputFile, tmpPath) inputTomoPath = tmpPath if self.applyToOddEven(tomo): pwpath.moveFile(oddEvenOutput[0], tmpPath) pwpath.moveFile(oddEvenOutput[1], tmpPath) inputTomoPath = tmpPath - inputOdd, inputEven = (os.path.join(tmpPrefix, tomo.getTsId() + EXT_MRC_ODD_NAME), - os.path.join(tmpPrefix, tomo.getTsId() + EXT_MRC_EVEN_NAME)) + inputOdd, inputEven = (self.getTmpOutFile(tsId, suffix=ODD), + self.getTmpOutFile(tsId, suffix=EVEN)) else: inputTomoPath = location if self.applyToOddEven(tomo): @@ -297,10 +332,10 @@ def generateOutputStackStep(self, tsObjId): if self.applyToOddEven(tomo): paramsBinvol['input'] = inputOdd - paramsBinvol['output'] = os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_ODD_NAME) + paramsBinvol['output'] = self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT) Plugin.runImod(self, 'binvol', argsBinvol % paramsBinvol) paramsBinvol['input'] = inputEven - paramsBinvol['output'] = os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_EVEN_NAME) + paramsBinvol['output'] = self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT) Plugin.runImod(self, 'binvol', argsBinvol % paramsBinvol) output = self.getOutputSetOfTomograms(self.inputSetOfTomograms.get(), @@ -329,8 +364,8 @@ def generateOutputStackStep(self, tsObjId): newTomogram.copyAttributes(tomo, '_origin') if self.applyToOddEven(tomo): - halfMapsList = [os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_ODD_NAME), - os.path.join(extraPrefix, tomo.getTsId() + EXT_MRC_EVEN_NAME)] + halfMapsList = [self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT), + self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT)] newTomogram.setHalfMaps(halfMapsList) output.append(newTomogram) diff --git a/imod/protocols/protocol_tomoProjection.py b/imod/protocols/protocol_tomoProjection.py index e1df8f76..eaf1d659 100644 --- a/imod/protocols/protocol_tomoProjection.py +++ b/imod/protocols/protocol_tomoProjection.py @@ -24,18 +24,15 @@ # * # ***************************************************************************** -import os - from pwem.objects import Transform from pyworkflow import BETA from pyworkflow.object import Set import pyworkflow.protocol.params as params -import pyworkflow.utils.path as path from pwem.emlib.image import ImageHandler import tomo.objects as tomoObj from .. import Plugin -from .protocol_base import ProtImodBase +from .protocol_base import ProtImodBase, MRCS_EXT class ProtImodTomoProjection(ProtImodBase): @@ -43,6 +40,14 @@ class ProtImodTomoProjection(ProtImodBase): Re-project a tomogram given a geometric description (axis and angles). More info: https://bio3d.colorado.edu/imod/doc/man/xyzproj.html + + This program will compute projections of a tomogram at a series of + tilts around either the X, the Y or the Z axis.\n + + A projection along a ray line is simply the average of the pixels in + the block along that line. However, rather than taking the values of + the pixels that lie near the ray, interpolation is used to sample den- + sity at points evenly spaced at one pixel intervals along the ray. """ _label = 'Tomo projection' @@ -59,25 +64,26 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTomograms', important=True, - label='Input set of tomograms') + label='Input set of tomograms to be projected') - form.addParam('minAngle', + line = form.addLine('Tilt angles (deg)', + help='Starting, ending, and increment tilt angle. Enter the same value for ' + 'starting and ending angle to get only one image') + + line.addParam('minAngle', params.FloatParam, default=-60.0, - label='Minimum angle of rotation', - help='Minimum angle of the projection range') + label='Minimum rotation') - form.addParam('maxAngle', + line.addParam('maxAngle', params.FloatParam, default=60.0, - label='Maximum angle of rotation', - help='Maximum angle of the projection range') + label='Maximum rotation') - form.addParam('stepAngle', + line.addParam('stepAngle', params.FloatParam, default=2.0, - label='Step angle', - help='Step angle of the projection range') + label='Step angle') form.addParam('rotationAxis', params.EnumParam, @@ -91,26 +97,23 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - for tomo in self.inputSetOfTomograms.get(): - self._insertFunctionStep(self.projectTomogram, tomo.getObjId()) - self._insertFunctionStep(self.generateOutputStackStep, tomo.getObjId()) + self._initialize() + for tsId in self.tomoDict.keys(): + self._insertFunctionStep(self.projectTomogram, tsId) + self._insertFunctionStep(self.generateOutputStackStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def projectTomogram(self, tomoObjId): - tomo = self.inputSetOfTomograms.get()[tomoObjId] - - tomoId = os.path.splitext(os.path.basename(tomo.getFileName()))[0] + def _initialize(self): + self.tomoDict = {tomo.getTsId(): tomo.clone() for tomo in self.inputSetOfTomograms.get()} - extraPrefix = self._getExtraPath(tomoId) - tmpPrefix = self._getTmpPath(tomoId) - path.makePath(tmpPrefix) - path.makePath(extraPrefix) + def projectTomogram(self, tsId): + self.genTsPaths(tsId) + tomo = self.tomoDict[tsId] paramsXYZproj = { 'input': tomo.getFileName(), - 'output': os.path.join(extraPrefix, - os.path.basename(tomo.getFileName())), + 'output': self.getExtraOutFile(tsId, ext=MRCS_EXT), 'axis': self.getRotationAxis(), 'angles': str(self.minAngle.get()) + ',' + str(self.maxAngle.get()) + ',' + @@ -124,35 +127,25 @@ def projectTomogram(self, tomoObjId): Plugin.runImod(self, 'xyzproj', argsXYZproj % paramsXYZproj) - def generateOutputStackStep(self, tomoObjId): - tomo = self.inputSetOfTomograms.get()[tomoObjId] - - tomoId = os.path.splitext(os.path.basename(tomo.getFileName()))[0] - - extraPrefix = self._getExtraPath(tomoId) - + def generateOutputStackStep(self, tsId): + tomo = self.tomoDict[tsId] output = self.getOutputSetOfTiltSeries(self.inputSetOfTomograms.get()) - - newTs = tomoObj.TiltSeries(tsId=tomoId) - + newTs = tomoObj.TiltSeries(tsId=tsId) newTs.setTsId(tomo.getTsId()) newTs.setAcquisition(tomo.getAcquisition()) - newTs.setTsId(tomoId) # Add origin to output tilt-series output.append(newTs) - tiltAngleList = self.getTiltAngleList() - + sRate = self.inputSetOfTomograms.get().getSamplingRate() for index in range(self.getProjectionRange()): newTi = tomoObj.TiltImage() newTi.setTiltAngle(tiltAngleList[index]) - newTi.setTsId(tomoId) + newTi.setTsId(tsId) newTi.setAcquisitionOrder(index + 1) newTi.setLocation(index + 1, - os.path.join(extraPrefix, - os.path.basename(tomo.getFileName()))) - newTi.setSamplingRate(self.inputSetOfTomograms.get().getSamplingRate()) + self.getExtraOutFile(tsId, ext=MRCS_EXT)) + newTi.setSamplingRate(sRate) newTs.append(newTi) ih = ImageHandler() diff --git a/imod/protocols/protocol_tomoReconstruction.py b/imod/protocols/protocol_tomoReconstruction.py index af84c07d..7bbd8e18 100644 --- a/imod/protocols/protocol_tomoReconstruction.py +++ b/imod/protocols/protocol_tomoReconstruction.py @@ -32,8 +32,7 @@ from tomo.objects import Tomogram from .. import Plugin -from .protocol_base import (ProtImodBase, EXT_MRC_ODD_NAME, EXT_MRC_EVEN_NAME, - EXT_MRCS_TS_EVEN_NAME, EXT_MRCS_TS_ODD_NAME) +from .protocol_base import ProtImodBase, TLT_EXT, ODD, EVEN, MRC_EXT class ProtImodTomoReconstruction(ProtImodBase): @@ -42,6 +41,23 @@ class ProtImodTomoReconstruction(ProtImodBase): More info: https://bio3d.colorado.edu/imod/doc/man/tilt.html + + It makes use of the IMOD tilt program. Tilt is a program for reconstructing + a three-dimensional object (tomogram) from a series of 2D projections. + The projections are assumed to arise from rotation about a fixed tilt axis, + subject to minor variations from this scheme.\n + + The program uses a number of different numerical strategies depending on + the complexity of the alignments needed to reconstruct the tomogram + and on whether the computation is being done by the central processing + unit (CPU) or the graphical processing unit (GPU). If there are no + local alignments being applied, then for processing on the CPU, the + program will do a preliminary stretching of each input line by the + cosine of the tilt angle. This stretching speeds up the direct back- + projection because each stretched input line is in register with the + lines of the output planes. When computing on the GPU, the program + does not use cosine stretching, thus avoiding the consequences of + interpolating the data twice.\n """ _label = 'Tomo reconstruction' @@ -55,64 +71,102 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTiltSeries', important=True, - label='Input set of tilt-series') + label='Tilt Series') form.addParam('tomoThickness', params.FloatParam, default=1000, label='Tomogram thickness (voxels)', important=True, - help='Size in voxels of the tomogram in the z ' + help='Size in voxels of the tomogram along the z ' 'axis (beam direction).') - form.addParam('tomoShiftX', - params.FloatParam, - default=0, - label='Tomogram shift in X', - help='This entry allows one to shift the reconstructed ' - 'slice in X before it is output. If the X shift ' - 'is positive, the slice will be shifted to the ' - 'right, and the output will contain the left part ' - 'of the whole potentially reconstructable area.') form.addParam('tomoWidth', params.IntParam, default=0, label='Tomogram width (voxels)', help='Number of pixels to cut out in X, centered on the middle in X. Leave 0 for default X.') - form.addParam('tomoShiftZ', - params.FloatParam, - default=0, - label='Tomogram shift in Z', - help='This entry allows one to shift the reconstructed ' - 'slice in Z before it is output. If the Z ' - 'shift is positive, the slice is shifted upward. ' - 'The Z entry is optional and defaults to 0 ' - 'when omitted.') - - form.addParam('angleOffset', - params.FloatParam, - default=0, - label='Angle offset', - help='Apply an angle offset in degrees to all tilt ' - 'angles. This offset positively rotates the ' - 'reconstructed sections anticlockwise.') - form.addParam('tiltAxisOffset', - params.FloatParam, - default=0, - label='Tilt axis offset', - help='Apply an offset to the tilt axis in a stack of ' - 'full-sized projection images, cutting the ' - 'X-axis at NX/2. + offset instead of NX/2. The ' - 'DELXX entry is optional and defaults to 0 ' - 'when omitted.') + lineShift = form.addLine('Tomogram shift (Å)', + expertLevel=params.LEVEL_ADVANCED, + help="This entry allows one to shift the reconstructed" + " slice in X or Z before it is output. If the " + " X shift is positive, the slice will be shifted to " + " the right, and the output will contain the left " + " part of the whole potentially reconstructable area. " + " If the Z shift is positive, the slice is shifted " + " upward. The Z entry is optional and defaults to 0 when " + " omitted.") + + lineShift.addParam('tomoShiftX', + params.FloatParam, + default=0, + label=' in X ') + + lineShift.addParam('tomoShiftZ', + params.FloatParam, + default=0, + label=' in Z ') + + lineoffSet = form.addLine('Offset (deg) of the ', + expertLevel=params.LEVEL_ADVANCED, + help="* Tilt angle offset: pply an angle offset in degrees to all tilt " + "angles. This offset positively rotates the reconstructed sections anticlockwise.\n" + "* Tilt axis offset: Apply an offset to the tilt axis in a stack of full-sized " + "projection images, cutting the X-axis at NX/2 + offset instead of NX/2.") + + lineoffSet.addParam('angleOffset', + params.FloatParam, + default=0, + label='Tilt angles ', + help='Apply an angle offset in degrees to all tilt ' + 'angles. This offset positively rotates the ' + 'reconstructed sections anticlockwise.') + + lineoffSet.addParam('tiltAxisOffset', + params.FloatParam, + default=0, + label='Tilt axis', + help='Apply an offset to the tilt axis in a stack of ' + 'full-sized projection images, cutting the ' + 'X-axis at NX/2 + offset instead of NX/2. The ' + 'DELXX entry is optional and defaults to 0 ' + 'when omitted.') form.addParam('superSampleFactor', params.IntParam, default=2, label='Super-sampling factor', expertLevel=params.LEVEL_ADVANCED, - help='Compute slices in pixels smaller by this factor to reduce artifacts.') + help='Compute slices in pixels smaller by this factor to reduce artifacts.' + 'Super-sampling refers to computing the back projection in a slice ' + 'larger by an integer factor in each dimension, which is done here in ' + 'one of two ways: by interpolating the projection data at smaller ' + 'intervals during backprojection, or by expanding the input lines by ' + 'the factor using sync interpolation (padding in Fourier space). ' + 'Super-sampling will reduce the rays along the projection angles that ' + 'appear to reflect from the edges of a Fourier transform of an X/Z slice.' + 'These rays result from back-projecting into discrete pixels and represent' + 'inappropriate information generated by the transitions between successive ' + 'pixels along a backprojection ray. Super-sampling by 2 will remove most ' + 'of these rays, especially oblique ones. The additional benefit (amount ' + 'of change in the image) of going from 2 to 3 is about 10% as large as ' + 'that of super-sampling by 2; it is about 3% going from 3 to 4 and 1.5% ' + 'going from 4 to 8 (the highest allowed value). The super-sampled slice ' + 'is reduced to the original resolution by cropping its Fourier transform. ' + 'Super-sampling alone with the first method does not reduce the rays in ' + 'the corners of the Fourier transform past 0.5 cycle/pixel. These rays ' + 'are also inappropriate since they originate from lower-frequency ' + 'information in the projection images, so they are removed from the ' + 'cropped Fourier transform before inverting it. This removal has an ' + 'added benefit about 1/3 as large as the benefit from supersampling by 2. ' + 'The effect of these removals is a subtle change in the noise, and a ' + 'benefit may show up only with subvolume averaging.\n' + 'The additional effect of expanding the input lines is to avoid attenu-' + 'ating frequencies past half Nyquist (0.25/pixel) any more than is ' + 'achieved with the falloff of the radial filter. This will be noticeable ' + 'in a tomogram and would be particularly helpful if setting the radial ' + 'filter cutoff higher than the default for subvolume averaging.\n') form.addParam('fakeInteractionsSIRT', params.IntParam, @@ -125,31 +179,34 @@ def _defineParams(self, form): 'Gaussian filter is applied at the high-frequency ' 'end of the filter. The functioning of this filter ' 'is described in: \n\t' - 'https://bio3d.colorado.edu/imod/doc/man/tilt.html') + 'https://bio3d.colorado.edu/imod/doc/man/tilt.html' + 'This entry corresponds to the imod parameter – FakeSIRTiterations') groupRadialFrequencies = form.addGroup('Radial filtering', - help='This entry controls low-pass ' - 'filtering with the radial weighting ' - 'function. The radial weighting ' - 'function is linear away from the ' - 'origin out to the distance in ' - 'reciprocal space specified by the ' - 'first value, followed by a Gaussian ' - 'fall-off determined by the ' + help='This entry controls low-pass filtering with the radial weighting ' + 'function. The radial weighting function is linear away from the ' + 'origin out to the distance in reciprocal space specified by the ' + 'first value, followed by a Gaussian fall-off determined by the ' 'second value.', expertLevel=params.LEVEL_ADVANCED) groupRadialFrequencies.addParam('radialFirstParameter', params.FloatParam, default=0.35, - label='First parameter', - help='Linear region value') + label='Cutoff linear region', + help='Linear region of the filter. This is the first' + 'parameter of the -RADIAL flag in IMOD tilt command.' + 'If the cutoff is greater than 1 the distances are interpreted ' + 'as pixels in Fourier space; otherwise they are treated as ' + 'frequencies in cycles/pixel, which range from 0 to 0.5. ' + 'Use a cutoff of 0.5 for no low-pass filtering.') groupRadialFrequencies.addParam('radialSecondParameter', params.FloatParam, default=0.035, - label='Second parameter', - help='Gaussian fall-off parameter') + label='Radial fall-off', + help='Gaussian fall-off parameter. The sigma (or standard deviation) ' + 'of the Gaussian is the second value times 0.707.') form.addHidden(params.USE_GPU, params.BooleanParam, @@ -177,35 +234,31 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self._failedTs = [] - - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.convertInputStep, ts.getObjId()) - self._insertFunctionStep(self.computeReconstructionStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputFailedSet, ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + self._insertFunctionStep(self.computeReconstructionStep, tsId) + self._insertFunctionStep(self.createOutputStep, tsId) + self._insertFunctionStep(self.createOutputFailedStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def convertInputStep(self, tsObjId, **kwargs): + def _initialize(self): + self._failedTs = [] + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + + def convertInputStep(self, tsId, **kwargs): # Considering swapXY is required to make tilt axis vertical oddEvenFlag = self.applyToOddEven(self.inputSetOfTiltSeries.get()) - super().convertInputStep(tsObjId, doSwap=True, oddEven=oddEvenFlag) - - @ProtImodBase.tryExceptDecorator - def computeReconstructionStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - firstItem = ts.getFirstItem() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + super().convertInputStep(tsId, doSwap=True, oddEven=oddEvenFlag) + # @ProtImodBase.tryExceptDecorator + def computeReconstructionStep(self, tsId): + ts = self.tsDict[tsId] paramsTilt = { - 'InputProjections': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'OutputFile': os.path.join(tmpPrefix, firstItem.parseFileName(extension=".rec")), - 'TiltFile': os.path.join(tmpPrefix, firstItem.parseFileName(extension=".tlt")), + 'InputProjections': self.getTmpOutFile(tsId), + 'OutputFile': self.getTmpOutFile(tsId, ext=MRC_EXT), + 'TiltFile': self.getExtraOutFile(tsId, ext=TLT_EXT), 'Thickness': self.tomoThickness.get(), 'FalloffIsTrueSigma': 1, 'Radial': str(self.radialFirstParameter.get()) + "," + str(self.radialSecondParameter.get()), @@ -255,21 +308,19 @@ def getArgs(): oddEvenTmp = [[], []] if self.applyToOddEven(ts): - oddFn = os.path.join(tmpPrefix, tsId + EXT_MRCS_TS_ODD_NAME) - paramsTilt['InputProjections'] = oddFn - oddEvenTmp[0] = os.path.join(tmpPrefix, firstItem.parseFileName(extension="_odd.rec")) + paramsTilt['InputProjections'] = self.getTmpOutFile(tsId, suffix=ODD) + oddEvenTmp[0] = self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT) paramsTilt['OutputFile'] = oddEvenTmp[0] Plugin.runImod(self, 'tilt', argsTilt % paramsTilt) - evenFn = os.path.join(tmpPrefix, tsId + EXT_MRCS_TS_EVEN_NAME) - paramsTilt['InputProjections'] = evenFn - oddEvenTmp[1] = os.path.join(tmpPrefix, firstItem.parseFileName(extension="_even.rec")) + paramsTilt['InputProjections'] = self.getTmpOutFile(tsId, suffix=EVEN) + oddEvenTmp[1] = self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT) paramsTilt['OutputFile'] = oddEvenTmp[1] Plugin.runImod(self, 'tilt', argsTilt % paramsTilt) paramsTrimVol = { - 'input': os.path.join(tmpPrefix, firstItem.parseFileName(extension=".rec")), - 'output': os.path.join(extraPrefix, firstItem.parseFileName(extension=".mrc")), + 'input': self.getTmpOutFile(tsId, ext=MRC_EXT), + 'output': self.getExtraOutFile(tsId, ext=MRC_EXT), 'options': getArgs() } @@ -281,21 +332,16 @@ def getArgs(): if self.applyToOddEven(ts): paramsTrimVol['input'] = oddEvenTmp[0] - paramsTrimVol['output'] = os.path.join(extraPrefix, tsId + EXT_MRC_ODD_NAME) + paramsTrimVol['output'] = self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT) Plugin.runImod(self, 'trimvol', argsTrimvol % paramsTrimVol) paramsTrimVol['input'] = oddEvenTmp[1] - paramsTrimVol['output'] = os.path.join(extraPrefix, tsId + EXT_MRC_EVEN_NAME) + paramsTrimVol['output'] = self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT) Plugin.runImod(self, 'trimvol', argsTrimvol % paramsTrimVol) - def createOutputStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - firstItem = ts.getFirstItem() - - extraPrefix = self._getExtraPath(tsId) - - tomoLocation = os.path.join(extraPrefix, firstItem.parseFileName(extension=".mrc")) + def createOutputStep(self, tsId): + ts = self.tsDict[tsId] + tomoLocation = self.getExtraOutFile(tsId, ext=MRC_EXT) if os.path.exists(tomoLocation): output = self.getOutputSetOfTomograms(self.inputSetOfTiltSeries.get()) @@ -304,8 +350,8 @@ def createOutputStep(self, tsObjId): newTomogram.setLocation(tomoLocation) if self.applyToOddEven(ts): - halfMapsList = [os.path.join(extraPrefix, tsId + EXT_MRC_ODD_NAME), - os.path.join(extraPrefix, tsId + EXT_MRC_EVEN_NAME)] + halfMapsList = [self.getExtraOutFile(tsId, suffix=ODD, ext=MRC_EXT), + self.getExtraOutFile(tsId, suffix=EVEN, ext=MRC_EXT)] newTomogram.setHalfMaps(halfMapsList) newTomogram.setTsId(tsId) @@ -325,6 +371,10 @@ def createOutputStep(self, tsObjId): output.write() self._store() + def createOutputFailedStep(self, tsId): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts) + def closeOutputSetsStep(self): for _, output in self.iterOutputAttributes(): output.setStreamState(Set.STREAM_CLOSED) diff --git a/imod/protocols/protocol_tsNormalization.py b/imod/protocols/protocol_tsPreprocess.py similarity index 69% rename from imod/protocols/protocol_tsNormalization.py rename to imod/protocols/protocol_tsPreprocess.py index cd390a9d..55c8a325 100644 --- a/imod/protocols/protocol_tsNormalization.py +++ b/imod/protocols/protocol_tsPreprocess.py @@ -23,29 +23,44 @@ # * e-mail address 'scipion@cnb.csic.es' # * # ***************************************************************************** - -import os from pyworkflow import BETA from pyworkflow.object import Set import pyworkflow.protocol.params as params from pwem.emlib.image import ImageHandler import tomo.objects as tomoObj - from .. import Plugin -from .protocol_base import ProtImodBase, EXT_MRCS_TS_EVEN_NAME, EXT_MRCS_TS_ODD_NAME, OUTPUT_TILTSERIES_NAME -from ..utils import formatTransformFile +from .protocol_base import ProtImodBase, OUTPUT_TILTSERIES_NAME, XF_EXT, ODD, EVEN +from ..utils import genXfFile -class ProtImodTSNormalization(ProtImodBase): +class ProtImodTsNormalization(ProtImodBase): """ Normalize input tilt-series and change its storing formatting. More info: https://bio3d.colorado.edu/imod/doc/man/newstack.html + + IMOD tilt series preprocess makes use of the Newstack command. + In particular, three functionalities are possible:\n + + _1 Binning_: The protocol also allows to bin tilt series. This + means to reduce the dimensions of the tilt series keeping but + keeping most of the information. The binning factor or simply + binning is an integer number and represent the scaling factor + of the images. Binning 2 means that the original images will + be twice the binned ones. + _2 Normalization_: This protocol allows to scale the gray values + of the images, also called normalization, to a common range or + mean of density. The most used normalization consists in zero + mean and standard deviation one.\n + + _3 storage format_: IMOD is able to modify the number of bit of + the stored data in order to reduce the disc occupancy. + """ _label = 'Tilt-series preprocess' _devStatus = BETA - _possibleOutputs = {OUTPUT_TILTSERIES_NAME:tomoObj.SetOfTiltSeries} + _possibleOutputs = {OUTPUT_TILTSERIES_NAME: tomoObj.SetOfTiltSeries} # -------------------------- DEFINE param functions ----------------------- def _defineParams(self, form): @@ -61,9 +76,12 @@ def _defineParams(self, form): default=1, label='Binning', important=True, - help='Binning to be applied to the normalized tilt-series ' - 'in IMOD convention. Images will be binned by the ' - 'given factor. Must be an integer bigger than 1') + help='Binning is an scaling factor for the output images. ' + 'Must be an integer greater than 1. IMOD uses ordinary' + 'binning to reduce images in size by the given factor. ' + 'The value of a binned pixel is the average of pixel ' + 'values in each block of pixels being binned. Binning ' + 'is applied before all') form.addParam('applyAlignment', params.BooleanParam, @@ -73,38 +91,98 @@ def _defineParams(self, form): form.addParam('floatDensities', params.EnumParam, - choices=['default', '1', '2', '3', '4'], + choices=['No adjust', + 'range between min and max', + 'scaled to common mean and standard deviation', + 'shifted to a common mean without scaling', + 'shifted to mean and rescaled to a min and max'], default=0, label='Adjust densities mode', - display=params.EnumParam.DISPLAY_HLIST, + display=params.EnumParam.DISPLAY_COMBO, help='Adjust densities of sections individually:\n' '-Default: no adjustment performed\n' - '-Mode 1: sections fill the data range\n' - '-Mode 2: sections scaled to common mean and standard deviation.\n' - '-Mode 3: sections shifted to a common mean without scaling\n' - '-Mode 4: sections shifted to a common mean and then ' - 'rescale the resulting minimum and maximum densities ' - 'to the Min and Max values specified') + + '-Range between min and max: This option will scale the gray values' + 'to be in a range given by a minimum and a maximum values.' + 'This is the mode 1 in newstack flag -floatDensities.\n' + + '-Scaled to common mean and standard deviation: This is the most ' + 'common normalization procedure. The new tilt series will have' + 'a meand and a standard deviation introduced by the user. Generaly,' + 'a zero meand and a standard deviation one is a good choice.' + 'This is the mode 2 in newstack flag -floatDensities.\n' + + '-Shifted to a common mean without scaling: This option only' + 'add an offset to the gray values of the images. The offset will' + 'be calculated such as the new images will present a mean gray value' + 'introduced by the user. This is the mode 3 in newstack flag ' + 'floatDensities.\n' + + '-shifted to mean and rescaled to a min and max: In this case, an ' + 'offset is added to the images in order to achieve a mean gray value' + ' then they are rescale the resulting minimum and maximum densities ' + 'to the Min and Max values specified. This is the mode 4 in newstack' + ' flag -floatDensities.\n') + + groupMeanSd = form.addGroup('Mean and SD', + condition='floatDensities==2', + help='Scale all images to the given mean ' + 'and standard deviation.') + + groupMeanSd.addParam('meanSdToggle', + params.EnumParam, + choices=['Yes', 'No'], + default=1, + label='Set mean and SD?', + display=params.EnumParam.DISPLAY_HLIST, + help='Set mean and SD values') + + groupMeanSd.addParam('scaleMean', + params.FloatParam, + default=0, + label='Mean', + help='Mean value for the rescaling') + + groupMeanSd.addParam('scaleSd', + params.FloatParam, + default=1, + label='SD', + help='Standard deviation value for the rescaling') + + groupScale = form.addGroup('Scaling values', + condition='floatDensities==4') + + groupScale.addParam('scaleMax', + params.FloatParam, + default=255, + label='Max.', + help='Maximum value for the rescaling') + + groupScale.addParam('scaleMin', + params.FloatParam, + default=0, + label='Min.', + help='Minimum value for the rescaling') form.addParam('modeToOutput', params.EnumParam, - choices=['default', '4-bit', 'byte', 'signed 16-bit', + expertLevel=params.LEVEL_ADVANCED, + choices=['default', '4-bit', '8-bit', 'signed 16-bit', 'unsigned 16-bit', '32-bit float'], default=0, label='Storage data type', - display=params.EnumParam.DISPLAY_HLIST, - help='Apply one density scaling to all sections to ' - 'map current min and max to the given Min and ' - 'Max. The storage mode of the output file. The ' + display=params.EnumParam.DISPLAY_COMBO, + help='The storage mode of the output file. The ' 'default is the mode of the first input file, ' 'except for a 4-bit input file, where the default ' 'is to output as bytes') form.addParam('scaleRangeToggle', params.EnumParam, + choices=['Yes', 'No'], - condition="floatDensities==0 or floatDensities==1 or floatDensities==3", - default=1, + condition="floatDensities==1 or floatDensities==3", + default=0, label='Set scaling range values?', display=params.EnumParam.DISPLAY_HLIST, help='This option will rescale the densities of all ' @@ -112,30 +190,31 @@ def _defineParams(self, form): 'minimum and maximum density will be mapped ' 'to the Min and Max values that are entered') - form.addParam('scaleRangeMax', - params.FloatParam, - condition="(floatDensities==0 or floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", - default=255, - label='Max.', - help='Maximum value for the rescaling') - form.addParam('scaleRangeMin', params.FloatParam, - condition="(floatDensities==0 or floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", + condition="(floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", default=0, label='Min.', help='Minimum value for the rescaling') + form.addParam('scaleRangeMax', + params.FloatParam, + condition="(floatDensities==1 or floatDensities==3) and scaleRangeToggle==0", + default=255, + label='Max.', + help='Maximum value for the rescaling') + form.addParam('antialias', params.EnumParam, + expertLevel=params.LEVEL_ADVANCED, choices=['None', 'Blackman', 'Triangle', 'Mitchell', - 'Lanczos 2', 'Lanczos 3'], + 'Lanczos 2 lobes', 'Lanczos 3 lobes'], default=5, label='Antialias method:', - display=params.EnumParam.DISPLAY_HLIST, + display=params.EnumParam.DISPLAY_COMBO, help='Type of antialiasing filter to use when reducing images.\n' 'The available types of filters are:\n\n' - 'None\n' + 'None - Antialias will not be applied\n' 'Blackman - fast but not as good at antialiasing as slower filters\n' 'Triangle - fast but smooths more than Blackman\n' 'Mitchell - good at antialiasing, smooths a bit\n' @@ -147,49 +226,6 @@ def _defineParams(self, form): 'based on images of natural scenes where there are ' 'sharp edges.') - groupMeanSd = form.addGroup('Mean and SD', - condition='floatDensities==2', - help='Scale all images to the given mean ' - 'and standard deviation. This option ' - 'implies -float 2 and is incompatible ' - 'with all other scaling options. If no ' - 'values are set, mean=0 and SD=1 by default') - - groupMeanSd.addParam('meanSdToggle', - params.EnumParam, - choices=['Yes', 'No'], - default=1, - label='Set mean and SD?', - display=params.EnumParam.DISPLAY_HLIST, - help='Set mean and SD values') - - groupMeanSd.addParam('scaleMean', - params.FloatParam, - default=0, - label='Mean', - help='Mean value for the rescaling') - - groupMeanSd.addParam('scaleSd', - params.FloatParam, - default=1, - label='SD', - help='Standard deviation value for the rescaling') - - groupScale = form.addGroup('Scaling values', - condition='floatDensities==4') - - groupScale.addParam('scaleMax', - params.FloatParam, - default=255, - label='Max.', - help='Maximum value for the rescaling') - - groupScale.addParam('scaleMin', - params.FloatParam, - default=0, - label='Min.', - help='Minimum value for the rescaling') - form.addParam('processOddEven', params.BooleanParam, expertLevel=params.LEVEL_ADVANCED, @@ -200,42 +236,41 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - self._failedTs = [] - - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.convertInputStep, ts.getObjId()) - self._insertFunctionStep(self.generateOutputStackStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputFailedSet, ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + self._insertFunctionStep(self.generateOutputStackStep, tsId) + self._insertFunctionStep(self.createOutputFailedStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def convertInputStep(self, tsObjId, **kwargs): + def _initialize(self): + self._failedTs = [] + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + + def convertInputStep(self, tsId, **kwargs): oddEvenFlag = self.applyToOddEven(self.inputSetOfTiltSeries.get()) # Interpolation will be done in the generateOutputStep - super().convertInputStep(tsObjId, imodInterpolation=None, - generateAngleFile=False, oddEven=oddEvenFlag) + super().convertInputStep(tsId, + imodInterpolation=None, + generateAngleFile=False, + oddEven=oddEvenFlag) @ProtImodBase.tryExceptDecorator - def generateOutputStackStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + def generateOutputStackStep(self, tsId): + ts = self.tsDict[tsId] firstItem = ts.getFirstItem() - xfFile = None if self.applyAlignment.get() and ts.hasAlignment(): - xfFile = os.path.join(tmpPrefix, firstItem.parseFileName(extension=".xf")) - formatTransformFile(ts, xfFile) + xfFile = self.getExtraOutFile(tsId, ext=XF_EXT) + genXfFile(ts, xfFile) binning = self.binning.get() argsNewstack, paramsNewstack = self.getBasicNewstackParams(ts, - os.path.join(extraPrefix, firstItem.parseFileName()), - inputTsFileName=os.path.join(tmpPrefix, - firstItem.parseFileName()), + self.getExtraOutFile(tsId), + inputTsFileName=self.getTmpOutFile(tsId), xfFile=xfFile, firstItem=firstItem, binning=binning, @@ -275,10 +310,10 @@ def generateOutputStackStep(self, tsObjId): oddFn = firstItem.getOdd().split('@')[1] evenFn = firstItem.getEven().split('@')[1] paramsNewstack['input'] = oddFn - paramsNewstack['output'] = os.path.join(extraPrefix, tsId + EXT_MRCS_TS_ODD_NAME) + paramsNewstack['output'] = self.getExtraOutFile(tsId, suffix=ODD) Plugin.runImod(self, 'newstack', argsNewstack % paramsNewstack) paramsNewstack['input'] = evenFn - paramsNewstack['output'] = os.path.join(extraPrefix, tsId + EXT_MRCS_TS_EVEN_NAME) + paramsNewstack['output'] = self.getExtraOutFile(tsId, suffix=EVEN) Plugin.runImod(self, 'newstack', argsNewstack % paramsNewstack) output = self.getOutputSetOfTiltSeries(self.inputSetOfTiltSeries.get(), self.binning.get()) @@ -304,14 +339,14 @@ def generateOutputStackStep(self, tsObjId): newTi.setAcquisition(tiltImage.getAcquisition()) if self.applyToOddEven(ts): - locationOdd = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_ODD_NAME)) - locationEven = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_EVEN_NAME)) + locationOdd = index + 1, self.getExtraOutFile(tsId, suffix=ODD) + locationEven = index + 1, self.getExtraOutFile(tsId, suffix=EVEN) newTi.setOddEven([ih.locationToXmipp(locationOdd), ih.locationToXmipp(locationEven)]) else: newTi.setOddEven([]) - newTi.setLocation(index + 1, - (os.path.join(extraPrefix, tiltImage.parseFileName()))) + newTi.setLocation(index + 1, self.getExtraOutFile(tsId)) + if binning > 1: newTi.setSamplingRate(tiltImage.getSamplingRate() * binning) newTs.append(newTi) @@ -325,6 +360,10 @@ def generateOutputStackStep(self, tsObjId): output.write() self._store() + def createOutputFailedStep(self, tsId): + ts = self.tsDict[tsId] + super().createOutputFailedSet(ts) + def closeOutputSetsStep(self): for _, output in self.iterOutputAttributes(): output.setStreamState(Set.STREAM_CLOSED) diff --git a/imod/protocols/protocol_xCorrPrealignment.py b/imod/protocols/protocol_xCorrPrealignment.py index 5c29b065..8376c803 100644 --- a/imod/protocols/protocol_xCorrPrealignment.py +++ b/imod/protocols/protocol_xCorrPrealignment.py @@ -24,7 +24,6 @@ # * # ***************************************************************************** -import os import numpy as np from pyworkflow import BETA @@ -35,7 +34,7 @@ import tomo.objects as tomoObj from .. import Plugin, utils -from .protocol_base import ProtImodBase +from .protocol_base import ProtImodBase, TLT_EXT, PREXF_EXT, PREXG_EXT class ProtImodXcorrPrealignment(ProtImodBase): @@ -43,6 +42,16 @@ class ProtImodXcorrPrealignment(ProtImodBase): Tilt-series cross correlation alignment based on the IMOD procedure. More info: https://bio3d.colorado.edu/imod/doc/man/tiltxcorr.html + + Tiltxcorr uses cross-correlation to find an initial translational + alignment between successive images of a tilt series. For a given pair + of images, it stretches the image with the larger tilt angle perpendic- + ular to the tilt axis, by an amount equal to the ratio of the cosines + of the two tilt angles (cosine stretch). The stretched image is corre- + lated with the other image, and the position of the peak of the corre- + lation indicates the relative shift between the images. + + """ _label = 'Coarse prealignment' @@ -57,7 +66,7 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTiltSeries', important=True, - label='Input set of tilt-series') + label='Tilt-series to be prealigned') form.addParam('cumulativeCorr', params.EnumParam, @@ -65,7 +74,7 @@ def _defineParams(self, form): default=1, label='Use cumulative correlation?', display=params.EnumParam.DISPLAY_HLIST, - help='With this option, the program will take the image at zero tilt as the first' + help='The program will take the image at zero tilt as the first' 'reference, and correlate it with the image at the next most negative tilt.' 'It will then add the aligned image to the first reference to make the ' 'reference for the next tilt. At each tilt, the reference will be the sum of ' @@ -79,18 +88,25 @@ def _defineParams(self, form): label='Generate interpolated tilt-series?', important=True, display=params.EnumParam.DISPLAY_HLIST, - help='Generate and save the interpolated tilt-series ' - 'applying the obtained transformation matrices.') + help='Generate and save the interpolated tilt-series applying the obtained transformation ' + 'matrices.\n' + 'By default, the output of this protocol will be a tilseries that will have associated' + 'the alignment information as a transformation matrix. When this option is set as Yes, ' + 'then a second output, called interpolated tilt series, is generated. The interpolated tilt ' + 'series should be used for visualization purpose but not for image processing') form.addParam('binning', params.IntParam, condition='computeAlignment==0', default=1, label='Binning for the interpolated', - help='Binning to be applied to the interpolated ' - 'tilt-series in IMOD convention. Images will be ' - 'binned by the given factor. Must be an integer ' - 'bigger than 1') + help='Binning to be applied to the interpolated tilt-series in IMOD ' + 'convention. \n' + 'Binning is an scaling factor given by an integer greater than 1. ' + 'IMOD uses ordinary binning to reduce images in size by the given factor. ' + 'The value of a binned pixel is the average of pixel values in each block ' + 'of pixels being binned. Binning is applied before all other image ' + 'transformations.') form.addParam('Trimming parameters', params.LabelParam, label='Tilt axis angle detected from import. In case another value is desired please adjust the ' @@ -106,138 +122,36 @@ def _defineParams(self, form): 'Usually, it will be 90 degrees less than the RotationAngle in a ' 'system with no axis inversions') - trimming = form.addGroup('Trimming parameters', - expertLevel=params.LEVEL_ADVANCED) - - xtrimming = trimming.addLine('Horizontal: Number of pixels to avoid from the', - expertLevel=params.LEVEL_ADVANCED, - help="Starting and ending X coordinates of a region to correlate, " - "based on the position of the region at zero tilt.") - - xtrimming.addParam('xmin', - params.IntParam, - label='left', - allowsNull=True, - expertLevel=params.LEVEL_ADVANCED) - - xtrimming.addParam('xmax', - params.IntParam, - label='right', - allowsNull=True, - expertLevel=params.LEVEL_ADVANCED) - - ytrimming = trimming.addLine('Vertical: Number of pixels to avoid from the', - expertLevel=params.LEVEL_ADVANCED, - help="Starting and ending Y coordinates of a region to correlate.") - - ytrimming.addParam('ymin', - params.IntParam, - label='top', - allowsNull=True, - expertLevel=params.LEVEL_ADVANCED) - - ytrimming.addParam('ymax', - params.IntParam, - label='botton', - allowsNull=True, - expertLevel=params.LEVEL_ADVANCED) - - filtering = form.addGroup('Filtering parameters', - expertLevel=params.LEVEL_ADVANCED) - - line1 = filtering.addLine('High pass filter', - expertLevel=params.LEVEL_ADVANCED, - help="Some high pass filtering, using a small value of Sigma1 such " - "as 0.03, may be needed to keep the program from being misled by very " - "large scale features in the images. If the images are noisy, some low " - "pass filtering with Sigma2 and Radius2 is appropriate (e.g. 0.05 for " - " Sigma2, 0.25 for Radius2). If the images are binned, these values " - "specify frequencies in the binned image, so a higher cutoff (less filtering) " - "might be appropriate.\n\n" - "" - "*FilterRadius1*: Low spatial frequencies in the cross-correlation " - "will be attenuated by a Gaussian curve that is 1 " - "at this cutoff radius and falls off below this " - "radius with a standard deviation specified by " - "FilterSigma2. Spatial frequency units range from " - "0 to 0.5.\n" - "*Filter sigma 1*: Sigma value to filter low frequencies in the " - "correlations with a curve that is an inverted " - "Gaussian. This filter is 0 at 0 frequency and " - "decays up to 1 with the given sigma value. " - "However, if a negative value of radius1 is entered, " - "this filter will be zero from 0 to " - "|radius1| then decay up to 1.") - - line1.addParam('filterRadius1', - params.FloatParam, - label='Filter radius 1', - default='0.0', - expertLevel=params.LEVEL_ADVANCED) - - line1.addParam('filterSigma1', - params.FloatParam, - label='Filter sigma 1', - default='0.03', - expertLevel=params.LEVEL_ADVANCED) - - line2 = filtering.addLine('Low pass filter', - expertLevel=params.LEVEL_ADVANCED, - help="If the images are noisy, some low " - "pass filtering with Sigma2 and Radius2 is appropriate (e.g. 0.05 for " - " Sigma2, 0.25 for Radius2). If the images are binned, these values " - "specify frequencies in the binned image, so a higher cutoff (less filtering) " - "might be appropriate.\n\n" - "*Filter radius 2*: High spatial frequencies in the cross-correlation " - "will be attenuated by a Gaussian curve that is 1 " - "at this cutoff radius and falls off above this " - "radius with a standard deviation specified by " - "FilterSigma2.\n" - "*Filter sigma 2*: Sigma value for the Gaussian rolloff below and " - "above the cutoff frequencies specified by " - "FilterRadius1 and FilterRadius2") - - line2.addParam('filterRadius2', - params.FloatParam, - label='Filter radius 2', - default='0.25', - expertLevel=params.LEVEL_ADVANCED) - - line2.addParam('filterSigma2', - params.FloatParam, - label='Filter sigma 2', - default='0.05', - expertLevel=params.LEVEL_ADVANCED) + trimming = form.addGroup('Trimming parameters', expertLevel=params.LEVEL_ADVANCED) + + self.trimimgForm(trimming, pxTrimCondition='False', + correlationCondition='True', levelType=params.LEVEL_ADVANCED) + self.filteringParametersForm(form, condition='True', levelType=params.LEVEL_ADVANCED) # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.convertInputStep, ts.getObjId()) - self._insertFunctionStep(self.computeXcorrStep, ts.getObjId()) - self._insertFunctionStep(self.generateOutputStackStep, - ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + self._insertFunctionStep(self.computeXcorrStep, tsId) + self._insertFunctionStep(self.generateOutputStackStep, tsId) if self.computeAlignment.get() == 0: - self._insertFunctionStep(self.computeInterpolatedStackStep, - ts.getObjId()) + self._insertFunctionStep(self.computeInterpolatedStackStep, tsId) self._insertFunctionStep(self.closeOutputSetsStep) # --------------------------- STEPS functions ----------------------------- - def computeXcorrStep(self, tsObjId): - """Compute transformation matrix for each tilt series""" - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + def _initialize(self): + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} + def computeXcorrStep(self, tsId): + """Compute transformation matrix for each tilt series""" + ts = self.tsDict[tsId] tiltAxisAngle = self.getTiltAxisOrientation(ts) paramsXcorr = { - 'input': os.path.join(tmpPrefix, - ts.getFirstItem().parseFileName()), - 'output': os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".prexf")), - 'tiltfile': os.path.join(tmpPrefix, - ts.getFirstItem().parseFileName(extension=".tlt")), + 'input': self.getTmpOutFile(tsId), + 'output': self.getExtraOutFile(tsId, ext=PREXF_EXT), + 'tiltfile': self.getExtraOutFile(tsId, ext=TLT_EXT), 'rotationAngle': tiltAxisAngle, 'filterSigma1': self.filterSigma1.get(), 'filterSigma2': self.filterSigma2.get(), @@ -284,39 +198,22 @@ def computeXcorrStep(self, tsObjId): Plugin.runImod(self, 'tiltxcorr', argsXcorr % paramsXcorr) paramsXftoxg = { - 'input': os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".prexf")), - 'goutput': os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".prexg")), + 'input': self.getExtraOutFile(tsId, ext=PREXF_EXT), + 'goutput': self.getExtraOutFile(tsId, ext=PREXG_EXT), } argsXftoxg = "-input %(input)s " \ "-NumberToFit 0 " \ "-goutput %(goutput)s " Plugin.runImod(self, 'xftoxg', argsXftoxg % paramsXftoxg) - def getTiltAxisOrientation(self, ts): - if self.tiltAxisAngle.get(): - return self.tiltAxisAngle.get() - else: - return ts.getAcquisition().getTiltAxisAngle() - - def generateOutputStackStep(self, tsObjId): + def generateOutputStackStep(self, tsId): """ Generate tilt-serie with the associated transform matrix """ - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - + ts = self.tsDict[tsId] output = self.getOutputSetOfTiltSeries(self.inputSetOfTiltSeries.get()) - - alignmentMatrix = utils.formatTransformationMatrix( - os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".prexg"))) - + alignmentMatrix = utils.formatTransformationMatrix(self.getExtraOutFile(tsId, ext=PREXG_EXT)) newTs = tomoObj.TiltSeries(tsId=tsId) newTs.copyInfo(ts) newTs.getAcquisition().setTiltAxisAngle(self.getTiltAxisOrientation(ts)) - output.append(newTs) for index, tiltImage in enumerate(ts): @@ -352,20 +249,14 @@ def generateOutputStackStep(self, tsObjId): self._store() - def computeInterpolatedStackStep(self, tsObjId): + def computeInterpolatedStackStep(self, tsId): output = self.getOutputInterpolatedSetOfTiltSeries(self.inputSetOfTiltSeries.get()) - - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + ts = self.tsDict[tsId] paramsAlignment = { - 'input': os.path.join(tmpPrefix, ts.getFirstItem().parseFileName()), - 'output': os.path.join(extraPrefix, ts.getFirstItem().parseFileName()), - 'xform': os.path.join(extraPrefix, - ts.getFirstItem().parseFileName(extension=".prexg")), + 'input': self.getTmpOutFile(tsId), + 'output': self.getExtraOutFile(tsId), + 'xform': self.getExtraOutFile(tsId, ext=PREXG_EXT), 'bin': self.binning.get(), 'imagebinned': 1.0 } @@ -392,9 +283,7 @@ def computeInterpolatedStackStep(self, tsObjId): for index, tiltImage in enumerate(ts): newTi = tomoObj.TiltImage() newTi.copyInfo(tiltImage, copyId=True) - newTi.setLocation(index + 1, - (os.path.join(extraPrefix, - tiltImage.parseFileName()))) + newTi.setLocation(index + 1, self.getExtraOutFile(tsId)) if self.binning > 1: newTi.setSamplingRate(tiltImage.getSamplingRate() * self.binning.get()) newTs.append(newTi) @@ -418,6 +307,13 @@ def closeOutputSetsStep(self): self._store() + # --------------------------- UTILS functions ------------------------------ + def getTiltAxisOrientation(self, ts): + if self.tiltAxisAngle.get(): + return self.tiltAxisAngle.get() + else: + return ts.getAcquisition().getTiltAxisAngle() + # --------------------------- INFO functions ------------------------------ def _summary(self): summary = [] diff --git a/imod/protocols/protocol_xRaysEraser.py b/imod/protocols/protocol_xRaysEraser.py index d8b032e6..59005dc4 100644 --- a/imod/protocols/protocol_xRaysEraser.py +++ b/imod/protocols/protocol_xRaysEraser.py @@ -24,8 +24,6 @@ # * # ***************************************************************************** -import os - from pyworkflow import BETA import pyworkflow.protocol.params as params from pyworkflow.object import Set @@ -33,14 +31,83 @@ import tomo.objects as tomoObj from .. import Plugin -from .protocol_base import ProtImodBase, OUTPUT_TILTSERIES_NAME, EXT_MRCS_TS_ODD_NAME, EXT_MRCS_TS_EVEN_NAME +from .protocol_base import (ProtImodBase, OUTPUT_TILTSERIES_NAME, + ODD, EVEN, MOD_EXT) class ProtImodXraysEraser(ProtImodBase): """ - Erase X-rays from aligned tilt-series based on the IMOD procedure. + Erase Xrays from aligned tilt-series based on the IMOD procedure. More info: https://bio3d.colorado.edu/imod/doc/man/ccderaser.html + + This program replaces deviant pixels with interpolated values from + surrounding pixels. It is designed to correct defects in electron + microscope images from CCD cameras. It can use two algorithms to + automatically remove peaks in intensity caused by X-rays.\n + + The automatic removal of X-rays works by dividing the area of each + image into patches for scanning. The mean and standard deviation (SD) + of the pixels in a patch are computed. The patch is then scanned for + pixels that deviate from the mean by more than a criterion number of + SDs (the scan criterion, a relatively low number to keep from missing + peaks). When such a pixel is found, the program searches neighboring + pixels to find a peak in intensity. It then computes the mean and SD + of pixels in an annulus around the peak and makes sure that the peak + deviates from this local mean by more than a criterion number of SDs + (the peak criterion). Neighboring pixels inside the inner radius of + the annulus are added to the list of pixels to be replaced if they + deviate by a lower criterion (the grow criterion). The patch of pixels + is then replaced by fitting a polynomial to adjacent pixels and inter- + polating from the polynomial. If the peak does not deviate suffi- + ciently from this local mean, but is stronger than the mean of the scan + area by the scan criterion plus 1, then the mean and SD is again com- + puted in a larger annulus. If the peak deviates from this mean by a + number of SDs bigger than another criterion for extra-large peaks, a + patch of pixels is found, but it is replaced only if enough pixels dif- + fer from adjacent ones by large enough amounts (see the -big option + below). The reason for these two stages is that the inner radius for + the first stage must be set safely smaller than the radius of gold + beads to avoid erasing part of the beads, whereas the second stage can + work with larger areas because it has more stringent criteria that + reject gold beads. + + After the peaks are found in a scanning patch, the program next finds + the difference between each pixel and the mean of the eight adjacent + pixels. The mean and SD of this difference is computed, then pixels + are sought that deviate from the mean by yet another criterion, the + difference criterion. When such a pixel is found, neighboring pixels + are examined and added to the patch of pixels to replace if their dif- + ference exceeds the grow criterion. If the number of pixels in the + patch does not exceed a specified maximum, replacement proceeds as + above; otherwise the patch is ignored. + + Two methods are used because the first method is probably more reliable + for dealing with strong peaks that extend over several pixels, while + the second method is definitely better for finding small X-rays. + + After all the patches have been scanned for a section, the program then + searches for single pixels with large interpixel differences at the + edges of the image, over the width set by the -border option. A dif- + ference between a pixel and the mean of whatever adjacent pixels exist + is computed and its deviation from the overall mean interpixel differ- + ence is divided by the maximum SD of interpixel differences over all of + the scans. When this value exceeds the difference criterion and the + interpixel difference is greater than that of its neighbors, the pixel + is replaced with the mean. This procedure is iterated up to 4 times to + catch adjacent extreme pixels. + + Tuning the removal of X-rays would primarily involve adjusting two of + the criteria. The peak and difference criteria would be adjusted down + or up to increase or decrease the number of deviant pixels that are + found. The grow criterion could also be adjusted down or up depending + on whether too few or too many pixels are included in a patch that is + replaced, but this step is not usually done in practice. If there are + strong, large artifacts that are not being removed, the big difference + criterion for extra-large peaks should be lowered first, then if neces- + sary, the maximum radius and criterion strength for extra-large peaks + can be adjusted. + """ _label = 'X-rays eraser' @@ -56,24 +123,24 @@ def _defineParams(self, form): params.PointerParam, pointerClass='SetOfTiltSeries', important=True, - label='Input set of tilt-series') + label='Tilt Series') form.addParam('peakCriterion', params.FloatParam, - default=8.0, - label='Peak criterion', + default=10.0, + label='Peak criterion (in std)', expertLevel=params.LEVEL_ADVANCED, help='Criterion # of SDs above local mean for erasing ' - 'peak based on intensity (the default is 8 SDs)') + 'peak based on intensity (the default is 10 SDs)') form.addParam('diffCriterion', params.FloatParam, - default=6.0, - label='Difference criterion', + default=8.0, + label='Difference criterion (in std)', expertLevel=params.LEVEL_ADVANCED, help='Criterion # of SDs above mean pixel-to-pixel ' 'difference for erasing a peak based on ' - 'differences (the default is 6 SDs).') + 'differences (the default is 8 SDs).') form.addParam('maximumRadius', params.FloatParam, @@ -86,7 +153,7 @@ def _defineParams(self, form): form.addParam('bigDiffCriterion', params.IntParam, default=19, - label='Big difference criterion', + label='Big difference criterion (in std)', expertLevel=params.LEVEL_ADVANCED, help='An extra-large peak will be erased only if the ' 'value for the maximum difference between ' @@ -111,29 +178,30 @@ def _defineParams(self, form): # -------------------------- INSERT steps functions ----------------------- def _insertAllSteps(self): - for ts in self.inputSetOfTiltSeries.get(): - self._insertFunctionStep(self.convertInputStep, ts.getObjId()) - self._insertFunctionStep(self.eraseXraysStep, ts.getObjId()) - self._insertFunctionStep(self.createOutputStep, ts.getObjId()) + self._initialize() + for tsId in self.tsDict.keys(): + self._insertFunctionStep(self.convertInputStep, tsId) + self._insertFunctionStep(self.eraseXraysStep, tsId) + self._insertFunctionStep(self.createOutputStep, tsId) self._insertFunctionStep(self.closeOutputStep) - def convertInputStep(self, tsObjId, **kwargs): - oddEvenFlag = self.applyToOddEven(self.inputSetOfTiltSeries.get()) - super().convertInputStep(tsObjId, imodInterpolation=None, - generateAngleFile=False, oddEven=oddEvenFlag) - - def eraseXraysStep(self, tsObjId): - ts = self.inputSetOfTiltSeries.get()[tsObjId] + def _initialize(self): + self.tsDict = {ts.getTsId(): ts.clone(ignoreAttrs=[]) for ts in self.inputSetOfTiltSeries.get()} - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - tmpPrefix = self._getTmpPath(tsId) + def convertInputStep(self, tsId, **kwargs): + oddEvenFlag = self.applyToOddEven(self.inputSetOfTiltSeries.get()) + super().convertInputStep(tsId, + imodInterpolation=None, + generateAngleFile=False, + oddEven=oddEvenFlag) + def eraseXraysStep(self, tsId): + ts = self.tsDict[tsId] firstItem = ts.getFirstItem() paramsCcderaser = { - 'input': os.path.join(tmpPrefix, firstItem.parseFileName()), - 'output': os.path.join(extraPrefix, firstItem.parseFileName()), + 'input': self.getTmpOutFile(tsId), + 'output': self.getExtraOutFile(tsId), 'findPeaks': 1, 'peakCriterion': self.peakCriterion.get(), 'diffCriterion': self.diffCriterion.get(), @@ -146,9 +214,7 @@ def eraseXraysStep(self, tsObjId): 'annulusWidth': 2.0, 'xyScanSize': 100, 'edgeExclusionWidth': 4, - 'pointModel': os.path.join(extraPrefix, - firstItem.parseFileName(suffix="_fid", - extension=".mod")), + 'pointModel': self.getExtraOutFile(tsId, suffix="fid", ext=MOD_EXT), 'borderSize': 2, 'polynomialOrder': 2, } @@ -176,19 +242,16 @@ def eraseXraysStep(self, tsObjId): oddFn = firstItem.getOdd().split('@')[1] evenFn = firstItem.getEven().split('@')[1] paramsCcderaser['input'] = oddFn - paramsCcderaser['output'] = os.path.join(extraPrefix, tsId + EXT_MRCS_TS_ODD_NAME) + paramsCcderaser['output'] = self.getExtraOutFile(tsId, suffix=ODD) Plugin.runImod(self, 'ccderaser', argsCcderaser % paramsCcderaser) paramsCcderaser['input'] = evenFn - paramsCcderaser['output'] = os.path.join(extraPrefix, tsId + EXT_MRCS_TS_EVEN_NAME) + paramsCcderaser['output'] = self.getExtraOutFile(tsId, suffix=EVEN) Plugin.runImod(self, 'ccderaser', argsCcderaser % paramsCcderaser) - def createOutputStep(self, tsObjId): + def createOutputStep(self, tsId): output = self.getOutputSetOfTiltSeries(self.inputSetOfTiltSeries.get()) - ts = self.inputSetOfTiltSeries.get()[tsObjId] - tsId = ts.getTsId() - extraPrefix = self._getExtraPath(tsId) - + ts = self.tsDict[tsId] newTs = tomoObj.TiltSeries(tsId=tsId) newTs.copyInfo(ts) output.append(newTs) @@ -199,12 +262,11 @@ def createOutputStep(self, tsObjId): newTi = tomoObj.TiltImage() newTi.copyInfo(tiltImage, copyId=True, copyTM=True) newTi.setAcquisition(tiltImage.getAcquisition()) - newTi.setLocation(index + 1, - (os.path.join(extraPrefix, - tiltImage.parseFileName()))) + newTi.setLocation(index + 1, self.getExtraOutFile(tsId)) + if self.applyToOddEven(ts): - locationOdd = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_ODD_NAME)) - locationEven = index + 1, (os.path.join(extraPrefix, tsId + EXT_MRCS_TS_EVEN_NAME)) + locationOdd = index + 1, self.getExtraOutFile(tsId, suffix=ODD) + locationEven = index + 1, self.getExtraOutFile(tsId, suffix=EVEN) newTi.setOddEven([ih.locationToXmipp(locationOdd), ih.locationToXmipp(locationEven)]) else: newTi.setOddEven([]) diff --git a/imod/tests/test_protocols_imod.py b/imod/tests/test_protocols_imod.py index ad94f3e5..943045f6 100644 --- a/imod/tests/test_protocols_imod.py +++ b/imod/tests/test_protocols_imod.py @@ -33,6 +33,9 @@ class TestImodBase(BaseTest): + atsId = 'a' + btsId = 'b' + @classmethod def setUpClass(cls): setupTestProject(cls) @@ -108,7 +111,7 @@ def _runTSNormalization(cls, inputSoTS, binning, floatDensities, modeToOutput, scaleRangeToggle, scaleRangeMax, scaleRangeMin, meanSdToggle, scaleMean, scaleSd, scaleMax, scaleMin): - cls.protTSNormalization = cls.newProtocol(ProtImodTSNormalization, + cls.protTSNormalization = cls.newProtocol(ProtImodTsNormalization, inputSetOfTiltSeries=inputSoTS, binning=binning, floatDensities=floatDensities, @@ -142,6 +145,7 @@ def _runFiducialModels(cls, inputSoTS, twoSurfaces, fiducialRadius, numberFiducial, rotationAngle, shiftsNearZeroFraction) -> ProtImodFiducialModel: cls.protFiducialAlignment = cls.newProtocol(ProtImodFiducialModel, + typeOfModel=0, inputSetOfTiltSeries=inputSoTS, twoSurfaces=twoSurfaces, fiducialRadius=fiducialRadius, @@ -151,6 +155,14 @@ def _runFiducialModels(cls, inputSoTS, twoSurfaces, fiducialRadius, cls.launchProtocol(cls.protFiducialAlignment) return cls.protFiducialAlignment + @classmethod + def _runFiducialModelsPT(cls, inputSoTS) -> ProtImodFiducialModel: + cls.protFiducialAlignment = cls.newProtocol(ProtImodFiducialModel, + typeOfModel=1, + inputSetOfTiltSeries=inputSoTS) + cls.launchProtocol(cls.protFiducialAlignment) + return cls.protFiducialAlignment + @classmethod def _runFiducialAlignemnt(cls, inputSoLM, twoSurfaces, rotationAngle, computeAlignment, binning) -> ProtImodFiducialAlignment: @@ -308,7 +320,7 @@ def setUpClass(cls): cls.inputTMFolder = os.path.split(cls.inputDataSet.getFile('tm1'))[0] - cls.excludeViewsOutputSizes = {'a': 57, 'b': 56} + cls.excludeViewsOutputSizes = {cls.atsId: 57, cls.btsId: 56} cls.binningTsNormalization = 2 @@ -382,6 +394,8 @@ def setUpClass(cls): rotationAngle=-12.5, shiftsNearZeroFraction=0.2) + cls.protFiducialModelsPT = cls._runFiducialModelsPT(inputSoTS=cls.protXcorr.TiltSeries) + cls.protFiducialAlignment = cls._runFiducialAlignemnt(inputSoLM=cls.protFiducialModels.FiducialModelGaps, twoSurfaces=0, rotationAngle=-12.5, @@ -446,7 +460,7 @@ def test_doseFilterOutputTS(self): tsId = ts.getFirstItem().getTsId() self.assertTrue(os.path.exists(os.path.join(self.protDoseFilter._getExtraPath(tsId), - "BB" + tsId + ".st"))) + tsId + ".mrcs"))) def test_xRaysEraserOutputTS(self): @@ -456,7 +470,7 @@ def test_xRaysEraserOutputTS(self): tsId = ts.getFirstItem().getTsId() self.assertTrue(os.path.exists(os.path.join(self.protXRaysEraser._getExtraPath(tsId), - "BB" + tsId + ".st"))) + tsId + ".mrcs"))) def test_excludeViewsOutputTS(self): @@ -466,7 +480,7 @@ def test_excludeViewsOutputTS(self): tsId = ts.getFirstItem().getTsId() self.assertTrue(os.path.exists(os.path.join(self.protExcludeViews._getExtraPath(tsId), - "BB" + tsId + ".st"))) + tsId + ".mrcs"))) for index, tsOut in enumerate(ts): self.assertEqual(tsOut.getSize(), self.excludeViewsOutputSizes[tsOut.getTsId()]) @@ -479,7 +493,7 @@ def test_normalizationOutputTS(self): tsId = ts.getFirstItem().getTsId() self.assertTrue(os.path.exists(os.path.join(self.protTSNormalization._getExtraPath(tsId), - "BB" + tsId + ".st"))) + tsId + ".mrcs"))) inSamplingRate = self.protTSNormalization.inputSetOfTiltSeries.get().getSamplingRate() outSamplingRate = ts.getSamplingRate() @@ -493,7 +507,7 @@ def test_prealignmentOutputTS(self): tsId = ts.getFirstItem().getTsId() outputLocation = os.path.join(self.protXcorr._getExtraPath(tsId), - "BB" + tsId + ".st") + tsId + ".mrcs") self.assertTrue(os.path.exists(outputLocation)) @@ -507,7 +521,7 @@ def test_prealignmentOutputInterpolatedTS(self): tsId = ts.getFirstItem().getTsId() outputLocation = os.path.join(self.protXcorr._getExtraPath(tsId), - "BB" + tsId + ".st") + tsId + ".mrcs") self.assertTrue(os.path.exists(outputLocation)) @@ -523,9 +537,9 @@ def test_fiducialModelstOutputFiducialModelGaps(self): tsId = output.getFirstItem().getTsId() outputLocationImod = os.path.join(self.protFiducialModels._getExtraPath(tsId), - "BB" + tsId + "_gaps.fid") + tsId + "_gaps.fid") outputLocationScipion = os.path.join(self.protFiducialModels._getExtraPath(tsId), - "BB" + tsId + "_gaps.sfid") + tsId + "_gaps.sfid") self.assertTrue(os.path.exists(outputLocationImod)) self.assertTrue(os.path.exists(outputLocationScipion)) @@ -538,7 +552,7 @@ def test_fiducialAlignmentOutputTS(self): tsId = output.getFirstItem().getTsId() outputLocation = os.path.join(self.protFiducialAlignment._getExtraPath(tsId), - "BB" + tsId + ".st") + tsId + ".mrcs") self.assertTrue(os.path.exists(outputLocation)) @@ -551,7 +565,7 @@ def test_fiducialAlignmentOutputTS(self): tsId = output.getFirstItem().getTsId() outputLocation = os.path.join(self.protFiducialAlignment._getExtraPath(tsId), - "BB" + tsId + ".st") + tsId + ".mrcs") self.assertTrue(os.path.exists(outputLocation)) @@ -567,7 +581,7 @@ def test_fiducialAlignmentOutputTS(self): tsId = output.getFirstItem().getTsId() outputLocation = os.path.join(self.protFiducialAlignment._getExtraPath(tsId), - "BB" + tsId + "_noGaps.sfid") + tsId + "_noGaps.sfid") self.assertTrue(os.path.exists(outputLocation)) @@ -577,7 +591,7 @@ def test_fiducialAlignmentOutputTS(self): tsId = output.getFirstItem().getTsId() outputLocation = os.path.join(self.protFiducialAlignment._getExtraPath(tsId), - "BB" + tsId + "_fid.xyz") + tsId + "_fid.xyz") self.assertTrue(os.path.exists(outputLocation)) @@ -594,7 +608,7 @@ def test_applyTransformationMatrixOutputInterpolatedTS(self): tsId = output.getFirstItem().getTsId() outputLocation = os.path.join(self.protApplyTransformationMatrix._getExtraPath(tsId), - "BB" + tsId + ".st") + tsId + ".mrcs") self.assertTrue(os.path.exists(outputLocation)) @@ -610,7 +624,7 @@ def test_tomoReconstructionOutputTomogram(self): tomoId = self.protTomoReconstruction.inputSetOfTiltSeries.get().getFirstItem().getTsId() outputLocation = os.path.join(self.protTomoReconstruction._getExtraPath(tomoId), - "BB" + tomoId + ".mrc") + tomoId + ".mrc") self.assertTrue(os.path.exists(outputLocation)) @@ -626,10 +640,10 @@ def test_goldBeadPeaker3DOutput(self): self.assertSetSize(output, size=52, diffDelta=30) tomoId = self.protGoldBeadPicker3D.inputSetOfTomograms.get().getFirstItem().getTsId() - location = self.protGoldBeadPicker3D._getExtraPath("BB" + tomoId) + location = self.protGoldBeadPicker3D._getExtraPath(tomoId) self.assertTrue(os.path.exists(os.path.join(location, - "BB" + tomoId + ".mod"))) + tomoId + ".mod"))) def test_tomoNormalizationOutput(self): @@ -656,7 +670,7 @@ def test_tomoProjectionOutputTiltSeriesSize(self): ih = ImageHandler() outputDimensions = ih.getDimensions(output.getFirstItem().getFirstItem().getFileName()) - self.assertEqual(outputDimensions, (256, 256, 61, 1)) + self.assertEqual(outputDimensions, (256, 256, 1, 61)) # Sampling rate inSamplingRate = self.protTomoProjection.inputSetOfTomograms.get().getSamplingRate() @@ -678,7 +692,7 @@ def setUpClass(cls): # Create links to the input tilt-series and its associated mdoc file to test the protocols with a set of two # elements to make the tests more robust linkTs = os.path.join(os.path.split(cls.inputSoTS)[0], - "WTI042413_1series4_copy.st") + "WTI042413_1series4_copy.mrc") if not os.path.exists(linkTs): path.createLink(cls.inputSoTS, linkTs) @@ -751,7 +765,7 @@ def test_ctfCorrectionOutput(self): for ts in output: tsId = ts.getTsId() outputLocation = os.path.join(self.protCTFCorrection._getExtraPath(tsId), - '%s.st' % tsId) + '%s.mrcs' % tsId) self.assertTrue(os.path.exists(outputLocation)) diff --git a/imod/utils.py b/imod/utils.py index 0f764ec5..e9266008 100644 --- a/imod/utils.py +++ b/imod/utils.py @@ -28,37 +28,51 @@ """ import logging import os - -logger = logging.getLogger(__name__) import csv import math +from typing import Union import numpy as np - import pyworkflow.object as pwobj import pyworkflow.utils as pwutils from imod import Plugin +from tomo.objects import TiltSeries +logger = logging.getLogger(__name__) -def formatTransformFile(ts, transformFilePath): - """ This method takes a tilt series and the output transformation file path - and creates an IMOD-based transform - file in the location indicated. """ - tsMatrixTransformList = [] +def genXfFile(ts: TiltSeries, outXfName: str, + presentAcqOrders: Union[set, None] = None, + onlyEnabled: bool = False) -> None: + """ This method takes a tilt series and the output transformation file path + and creates an IMOD-based transform file in the location indicated. The transformation matrix + of a tilt-image is only added if its acquisition order is contained in a set composed of the + acquisition orders present in both the given tilt-series and the CTFTomoSeries. If presentAcqOrders + is not None, it is considered before the attribute onlyEnabled, as presentAcqOrders may also have been + generated considering the enabled elements of the intersection. + """ - for ti in ts: - transform = ti.getTransform().getMatrix().flatten() + def formatMatrix(tiltImage): + transform = tiltImage.getTransform().getMatrix().flatten() transformIMOD = ['%.7f' % transform[0], '%.7f' % transform[1], '%.7f' % transform[3], '%.7f' % transform[4], "{:>6}".format('%.3g' % transform[2]), "{:>6}".format('%.3g' % transform[5])] - tsMatrixTransformList.append(transformIMOD) + return transformIMOD - with open(transformFilePath, 'w') as f: + if presentAcqOrders: + tsMatrixList = [formatMatrix(ti) for ti in ts if ti.getAcquisitionOrder() in presentAcqOrders] + else: + tsMatrixList = [] + for ti in ts: + if onlyEnabled and not ti.isEnabled(): + continue + tsMatrixList.append(formatMatrix(ti)) + + with open(outXfName, 'w') as f: csvW = csv.writer(f, delimiter='\t') - csvW.writerows(tsMatrixTransformList) + csvW.writerows(tsMatrixList) def formatTransformFileFromTransformList(transformMatrixList, transformFilePath): @@ -165,7 +179,7 @@ def generateIMODFiducialTextFile(landmarkModel, outputFilePath): for vector in infoTable: outputLines.append("\t%s\t%s\t%s\t%s\n" % (vector[3], vector[0], - vector[1], int(vector[2])-1)) + vector[1], int(vector[2]) - 1)) with open(outputFilePath, 'w') as f: f.writelines(outputLines) @@ -602,7 +616,7 @@ def refactorCTFDefocusAstigmatismPhaseShiftCutOnFreqEstimationInfo(ctfInfoIMODTa def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, - isRelion=False, inputTiltSeries=None): + isRelion=False, inputTiltSeries=None, presentAcqOrders=None): """ This method takes a ctfTomoSeries object a generate a defocus information file in IMOD formatting containing the same information in the specified location. """ @@ -633,8 +647,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, lines = [] pattern = "%d\t%d\t%.2f\t%.2f\t%d\n" - for index in defocusUDict.keys(): - + for index in sorted(defocusUDict.keys()): if index + nEstimationsInRange > len(defocusUDict.keys()): break @@ -647,7 +660,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, round(tiltSeries[index + ctfTomoSeries.getNumberOfEstimationsInRange()].getTiltAngle(), 2), round(tiltSeries[index].getTiltAngle(), 2), # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - int(float(defocusUDict[index][0])/10) + int(float(defocusUDict[index][0]) / 10) )) lines.append(newLine) @@ -672,7 +685,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, # order to match the IMOD defocus file format. lines = ["1\t0\t0.0\t0.0\t0.0\t3\n"] - for index in defocusUDict.keys(): + for index in sorted(defocusUDict.keys()): if index + ctfTomoSeries.getNumberOfEstimationsInRange() > len(defocusUDict.keys()): break @@ -686,9 +699,9 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, round(tiltSeries[index + ctfTomoSeries.getNumberOfEstimationsInRange()].getTiltAngle(), 2), round(tiltSeries[index].getTiltAngle(), 2), # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusUDict[index][0])/10, + float(defocusUDict[index][0]) / 10, # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusVDict[index][0])/10, + float(defocusVDict[index][0]) / 10, float(defocusAngleDict[index][0]), )) @@ -708,7 +721,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, # order to match the IMOD defocus file format. lines = ["4\t0\t0.0\t0.0\t0.0\t3\n"] - for index in defocusUDict.keys(): + for index in sorted(defocusUDict.keys()): if index + ctfTomoSeries.getNumberOfEstimationsInRange() > len(defocusUDict.keys()): break @@ -722,7 +735,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, round(tiltSeries[index + ctfTomoSeries.getNumberOfEstimationsInRange()].getTiltAngle(), 2), round(tiltSeries[index].getTiltAngle(), 2), # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusUDict[index][0])/10, + float(defocusUDict[index][0]) / 10, float(phaseShiftDict[index][0]), )) @@ -743,9 +756,9 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, with open(defocusFilePath, 'w') as f: # This line is added at the beginning of the file in order # to match the IMOD defocus file format - lines = ["5\t0\t0.0\t0.0\t0.0\t3\n"] + lines = ["5\t0\t0.0\t0.0\t0.0\t3\n"] - for index in defocusUDict.keys(): + for index in sorted(defocusUDict.keys()): if index + ctfTomoSeries.getNumberOfEstimationsInRange() > len(defocusUDict.keys()): break @@ -759,9 +772,9 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, round(tiltSeries[index + ctfTomoSeries.getNumberOfEstimationsInRange()].getTiltAngle(), 2), round(tiltSeries[index].getTiltAngle(), 2), # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusUDict[index][0])/10, + float(defocusUDict[index][0]) / 10, # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusVDict[index][0])/10, + float(defocusVDict[index][0]) / 10, float(defocusAngleDict[index][0]), float(phaseShiftDict[index][0]) )) @@ -785,7 +798,7 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, # to match the IMOD defocus file format lines = ["37\t0\t0.0\t0.0\t0.0\t3\n"] - for index in defocusUDict.keys(): + for index in sorted(defocusUDict.keys()): if index + ctfTomoSeries.getNumberOfEstimationsInRange() > len(defocusUDict.keys()): break @@ -799,9 +812,9 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, round(tiltSeries[index + ctfTomoSeries.getNumberOfEstimationsInRange()].getTiltAngle(), 2), round(tiltSeries[index].getTiltAngle(), 2), # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusUDict[index][0])/10, + float(defocusUDict[index][0]) / 10, # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - float(defocusVDict[index][0])/10, + float(defocusVDict[index][0]) / 10, float(defocusAngleDict[index][0]), float(phaseShiftDict[index][0]), float(cutOnFreqDict[index][0]) @@ -818,30 +831,30 @@ def generateDefocusIMODFileFromObject(ctfTomoSeries, defocusFilePath, else: # There is no information available as list (not an IMOD CTF estimation) - logger.debug("Defocus file generated form a defocus attributes.") + logger.info("Defocus file generated from defocus attributes.") with open(defocusFilePath, 'w') as f: lines = ["1\t0\t0.0\t0.0\t0.0\t3\n"] + ind = 1 + for ti in tiltSeries: + ctfTomo = ctfTomoSeries.getCtfTomoFromTi(ti) + if ctfTomo: + if presentAcqOrders and ctfTomo.getAcquisitionOrder() not in presentAcqOrders: + continue + tiltAngle = ti.getTiltAngle() + newLine = ("%d\t%d\t%.2f\t%.2f\t%.1f\t%.1f\t%.2f\n" % ( + ind, + ind, + tiltAngle, + tiltAngle, + # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) + ctfTomo.getDefocusU() / 10, + # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) + ctfTomo.getDefocusV() / 10, + ctfTomo.getDefocusAngle())) - # CtfTomoSeries is iterated inversely because IMOD set indexes - # upside down Scipion (highest index for - # the tilt-image with the highest negative angle) - for ctfTomo in ctfTomoSeries: - index = ctfTomo.getIndex().get() - - newLine = ("%d\t%d\t%.2f\t%.2f\t%.1f\t%.1f\t%.2f\n" % ( - index, - index, - tiltSeries[index].getTiltAngle(), - tiltSeries[index].getTiltAngle(), - # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - ctfTomo.getDefocusU()/10, - # CONVERT DEFOCUS VALUE TO NANOMETERS (IMOD CONVENTION) - ctfTomo.getDefocusV()/10, - ctfTomo.getDefocusAngle()) - ) - - lines.append(newLine) + lines.append(newLine) + ind += 1 f.writelines(lines) @@ -857,7 +870,7 @@ def generateDefocusUDictionary(ctfTomoSeries): else ctfTomo.getDefocusVList() defocusInfoList = defocusInfoList.split(",") - index = ctfTomo.getIndex().get() + index = ctfTomo.getIndex() defocusUDict[index] = defocusInfoList @@ -875,7 +888,7 @@ def generateDefocusVDictionary(ctfTomoSeries): defocusInfoList = ctfTomo.getDefocusVList() defocusInfoList = defocusInfoList.split(",") - index = ctfTomo.getIndex().get() + index = ctfTomo.getIndex() defocusVDict[index] = defocusInfoList @@ -893,7 +906,7 @@ def generateDefocusAngleDictionary(ctfTomoSeries): defocusAngleList = ctfTomo.getDefocusAngleList() defocusAngleList = defocusAngleList.split(",") - index = ctfTomo.getIndex().get() + index = ctfTomo.getIndex() defocusAngleDict[index] = defocusAngleList @@ -911,7 +924,7 @@ def generatePhaseShiftDictionary(ctfTomoSeries): phaseShiftList = ctfTomo.getPhaseShiftList() phaseShiftList = phaseShiftList.split(",") - index = ctfTomo.getIndex().get() + index = ctfTomo.getIndex() phaseShiftDict[index] = phaseShiftList @@ -929,7 +942,7 @@ def generateCutOnFreqDictionary(ctfTomoSeries): cutOnFreqList = ctfTomo.getCutOnFreqList() cutOnFreqList = cutOnFreqList.split(",") - index = ctfTomo.getIndex().get() + index = ctfTomo.getIndex() cutOnFreqDict[index] = cutOnFreqList @@ -966,7 +979,7 @@ def calculateRotationAngleFromTM(ts): tm = ti.getTransform().getMatrix() cosRotationAngle = tm[0][0] sinRotationAngle = tm[1][0] - avgRotationAngle += math.degrees(math.atan(sinRotationAngle/cosRotationAngle)) + avgRotationAngle += math.degrees(math.atan(sinRotationAngle / cosRotationAngle)) avgRotationAngle = avgRotationAngle / ts.getSize() @@ -984,7 +997,7 @@ def generateDoseFileFromDoseTS(ts, doseFileOutputPath): for ti in ts: acq = ti.getAcquisition() - doseInfoList.append((acq.getAccumDose()-acq.getDosePerFrame(), acq.getDosePerFrame())) + doseInfoList.append((acq.getAccumDose() - acq.getDosePerFrame(), acq.getDosePerFrame())) np.savetxt(doseFileOutputPath, np.asarray(doseInfoList), fmt='%f', delimiter=" ") diff --git a/imod/viewers/viewers.py b/imod/viewers/viewers.py index 794b6255..13b97867 100644 --- a/imod/viewers/viewers.py +++ b/imod/viewers/viewers.py @@ -58,6 +58,7 @@ class ImodViewer(pwviewer.Viewer): tomoObj.SetOfLandmarkModels, tomoObj.LandmarkModel ] + _name = 'Imod' def _visualize(self, obj, **kwargs): env = Plugin.getEnviron() diff --git a/requirements.txt b/requirements.txt index d0dcd66b..5ab7c8c3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1 @@ -scipion-em-tomo \ No newline at end of file +scipion-em-tomo>=3.7.0 \ No newline at end of file diff --git a/setup.py b/setup.py index 28064e18..ae719cc7 100644 --- a/setup.py +++ b/setup.py @@ -111,7 +111,7 @@ # 3 - Alpha # 4 - Beta # 5 - Production/Stable - 'Development Status :: 4 - Beta', + 'Development Status :: 5 - Production/Stable', # Indicate who your project is intended for # 'Intended Audience :: Users',