# **************************************************************************
# *
# * Authors: Carlos Oscar Sorzano (coss@cnb.csic.es)
# *
# * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 2 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307 USA
# *
# * All comments concerning this program package may be sent to the
# * e-mail address 'scipion@cnb.csic.es'
# *
# **************************************************************************
"""
Protocol to perform high-resolution reconstructions
"""
from glob import glob
import math
from xmipp3.constants import CUDA_ALIGN_SIGNIFICANT
try:
from itertools import izip
except ImportError:
izip = zip
from os.path import join, exists, split
import os
from pyworkflow import VERSION_1_1
from pyworkflow.protocol.constants import LEVEL_ADVANCED
from pyworkflow.protocol.params import (PointerParam, StringParam, FloatParam,
BooleanParam, IntParam, EnumParam,
NumericListParam, USE_GPU, GPU_LIST)
from pyworkflow.utils.path import (cleanPath, makePath, copyFile, moveFile,
createLink, cleanPattern)
from pwem.protocols import ProtRefine3D
from pwem.objects import SetOfVolumes, Volume
from pwem.emlib.metadata import getFirstRow, getSize
from pyworkflow.utils.utils import getFloatListFromValues
from pwem.emlib.image import ImageHandler
import pwem.emlib.metadata as md
from pwem.constants import ALIGN_PROJ
from pwem import emlib
from xmipp3.base import HelicalFinder, isXmippCudaPresent
from xmipp3.convert import createItemMatrix, setXmippAttributes, writeSetOfParticles
from pyworkflow import UPDATED, PROD
[docs]def getPreviousQuality(img, imgRow):
if hasattr(img,"_xmipp_cost"):
imgRow.setValue(md.MDL_COST,img._xmipp_cost.get())
if hasattr(img,"_xmipp_maxCC"):
imgRow.setValue(md.MDL_MAXCC,img._xmipp_maxCC.get())
[docs]class XmippProtReconstructHighRes(ProtRefine3D, HelicalFinder):
"""This is a 3D refinement protocol whose main input is a volume and a set of particles.
The set of particles has to be at full size (the finer sampling rate available), but
the rest of inputs (reference volume and masks) can be at any downsampling factor.
The protocol scales the input images and volumes to a reasonable size depending on
the resolution of the previous iteration.
The protocol works with any input volume, whichever its resolution, as long as it
is a reasonable initial volume for the set of particles. The protocol does not
resolve the heterogeneous problem (it assumes an homogeneous population),
although it is somewhat tolerant through the use of particle weights in the
reconstruction process.
It is recommended to perform several global alignment iterations before entering
into the local iterations. The switch from global to local should be performed when
a substantial percentage of the particles do not move from one iteration to the next.
The algorithm reports the cross correlation (global alignment) or cost (local) function
per defocus group, so that we can see which was the percentile of each particle in its
defocus group. You may want to perform iterations one by one, and remove from one
iteration to the next, those particles that worse fit the model.
AI Generated
## Overview
The Highres protocol performs iterative high-resolution 3D refinement of a
particle set against one or two reference volumes.
This protocol is intended for the later stages of single-particle cryo-EM
processing, when the user already has a reasonable initial volume and wants to
improve the angular assignment, reconstruction, and map quality. It assumes a
mostly homogeneous particle population, although it can reduce the influence of
poorly fitting particles through several weighting schemes.
The protocol works iteratively. In each iteration, it prepares particles and
references at a sampling rate appropriate for the current estimated resolution,
assigns particle orientations, weights and qualifies particles, reconstructs
two half maps, post-processes the result, evaluates the FSC, and prepares the
next reference.
The main outputs are a refined volume with associated half maps and a particle
set containing the final projection-alignment parameters and Xmipp quality
attributes.
## Inputs and General Workflow
The protocol requires full-size particle images and, unless continuing from a
previous run, one or more initial volumes.
The particle images should be provided at the finest available sampling rate.
The reference volumes and masks do not need to be at the same sampling rate as
the particles; the protocol rescales images, volumes, and masks internally
according to the resolution reached in each iteration.
The workflow is organized into iterations. Each iteration typically performs:
- global, local, stochastic, automatic, or no angular assignment;
- particle weighting;
- particle qualification by alignment quality;
- reconstruction of two half maps;
- post-processing of the reconstructed volume;
- FSC-based resolution evaluation;
- cleanup of intermediate files when requested.
This iterative structure allows the refinement to move progressively from
coarser to finer information.
## Continue from a Previous Run
The **Continue from a previous run?** option allows the protocol to resume from
an earlier Highres run.
When this option is enabled, the user selects a previous run, and several
parameters and intermediate files are reused. This is useful when the user
wants to extend the number of iterations, change later-stage settings, or
continue after a partial workflow.
If new input particles are not provided, the particle metadata from the
previous run are reused. If new particles are provided, they are converted and
used in the continuation.
This option is useful for iterative refinement strategies where the user
examines one run and then continues with adjusted parameters.
## Full-Size Images
The **Full-size Images** parameter defines the particle set to be refined.
The particles should be at full resolution, meaning the finest sampling rate
available for the project. The protocol may internally downsample them during
early iterations, but it keeps the full-size images as the original source.
The particle set should be reasonably clean and should correspond to a mostly
homogeneous structural population. The protocol is not designed to solve
strong heterogeneity by itself.
If the input particles already contain useful quality values from a previous
run, these may be carried into the metadata in some continuation or
no-alignment workflows.
## Initial Volumes
The **Initial volumes** parameter provides the starting reference for the
refinement.
The input can be a single volume or a set of two volumes. If two volumes are
provided, they are used as the initial references for the two half-map
branches.
If no initial volume is provided and the input particles already have angular
assignments, the protocol can generate an initial reconstruction from the
particles using Fourier reconstruction.
The initial volume does not need to be high resolution, but it should be a
reasonable representation of the particle. A poor or incorrect starting volume
can lead to incorrect angular assignments and poor convergence.
## Particle Radius
The **Radius of particle** parameter defines the radius, in pixels, of a
spherical mask covering the particle in the input images.
If the value is -1 or otherwise not positive, the protocol uses approximately
half the particle box size.
This radius is used in several internal operations, such as masking particles
and references. It should cover the particle density while excluding as much
empty background as possible.
For elongated or non-spherical particles, the spherical mask is only an
approximation, but it is still useful for limiting the region used during
alignment and reconstruction.
## Symmetry Group
The **Symmetry group** parameter defines the symmetry imposed during angular
assignment and reconstruction.
Use **c1** when no symmetry is present. If the particle has a known symmetry,
the appropriate Xmipp symmetry group can be provided.
Correct symmetry can improve signal and stability. Incorrect symmetry can
produce misleading maps by forcing together non-equivalent density. Symmetry
should therefore be imposed only when justified by the biological structure.
## Correct CTF Envelope
The **Correct CTF envelope** option controls whether the CTF envelope is
corrected during refinement.
Envelope correction can be useful when the envelope has been estimated
reliably. If the envelope is not reliable, applying this correction may amplify
noise or introduce artifacts.
This option should be used cautiously and according to the quality of the CTF
metadata.
## Remove Intermediate Files
The **Remove intermediate files** option helps save disk space by deleting
temporary files that are no longer needed.
High-resolution refinement can generate many intermediate files, especially
projection galleries, angle metadata, temporary stacks, and iteration
directories. Removing them can substantially reduce storage requirements.
Users who need to debug the protocol or inspect intermediate results may want
to disable this option temporarily.
## Next Reference Preparation
At the end of each iteration, the protocol prepares the reference volumes for
the next iteration.
This preparation can include FSC-based filtering, additional low-pass
filtering, spherical or cylindrical masking, positivity enforcement, user
masks, dropout, and optional custom commands.
The goal is to prevent the next angular assignment from using unreliable
high-frequency information or noisy regions of the map.
## Low-Pass Filter for the Next Reference
The **Low pass filter?** option controls whether the reference for the next
iteration is filtered according to the current estimated resolution.
The maximum frequency used for the next reference is based on the current
resolution plus the **Resolution offset**. A positive offset makes the
reference more conservative by using less high-resolution information. A
negative offset makes the process more aggressive by allowing more information.
This mechanism helps reduce overfitting during refinement.
## FSC Criterion and Resolution Offset
The **FSC criterion** defines how the protocol estimates resolution from the
FSC between the two half maps. Common values are 0.143 and 0.5.
The **Resolution offset** modifies the low-pass filtering applied to the next
reference. For example, if the current resolution is 8 Å and the offset is 2 Å,
the next reference may be filtered around 10 Å.
A conservative offset can make refinement more stable, especially in early
iterations. A more aggressive setting may converge faster but can increase the
risk of overfitting or instability.
## Masks for the Next Reference
The protocol can apply several types of masking to the next reference.
The **Spherical mask?** option applies a spherical mask based on the particle
radius. If helical symmetry post-processing is enabled, a cylindrical mask may
be used instead.
The **Mask** parameter allows the user to provide an external volume mask.
Smooth masks are recommended, and values should range from 0 to 1.
Masks should include the relevant particle density while excluding solvent and
irrelevant regions. Overly tight masks can remove real density, whereas overly
loose masks may allow noise to influence alignment.
## Positivity and Dropout
The **Positivity?** option removes negative values from the next reference.
This can help stabilize refinement when negative density is not meaningful for
the map.
The **Dropout** parameter randomly drops voxels inside the binary mask with a
given probability. This is an advanced regularization option that can reduce
dependence on specific reference features.
Most users should keep dropout at 0 unless they have a specific reason to use
it.
## Angular Assignment Modes
The **Image alignment** parameter defines how particle orientations are updated.
The available modes are:
**Global** searches over a projection gallery and is useful when orientations
are still uncertain.
**Local** refines orientations and other parameters around existing angular
assignments.
**Automatic** starts with global alignment and then switches to local alignment
in later iterations.
**Stochastic** performs global alignment using random subsets of images.
**No alignment** reconstructs using existing particle orientations without
changing them.
The appropriate mode depends on the stage of refinement. Early stages usually
benefit from global alignment. Later high-resolution refinement usually
benefits from local alignment.
## Number of Iterations
The **Number of iterations** parameter controls how many refinement cycles are
performed.
For global, local, and stochastic modes, the user specifies the number of
iterations directly. In automatic mode, the protocol uses a predefined
transition strategy. In no-alignment mode, only one iteration is needed.
It is often better to run a moderate number of iterations, inspect the results,
and then continue from the previous run if further refinement is needed.
## Multiresolution Approach
The **Multiresolution approach** adapts the sampling rate of the images to the
current estimated resolution.
In early iterations, when the map is low resolution, the protocol can work with
downsampled images. As the resolution improves, it moves toward finer sampling.
This reduces computational cost and helps the refinement focus on information
that is meaningful at the current stage.
The **Max. Target Resolution** list defines the target resolution schedule in
angstroms. The protocol also limits how quickly the resolution can improve
between iterations.
## Angular Search Limits
The **Max. shift** parameter defines the maximum allowed particle shift as a
percentage of the image size.
The **Tilt angle** parameters define the minimum and maximum tilt angles used
during angular assignment. In this convention, 0 degrees corresponds to top
views and 90 degrees to side views.
Restricting tilt angles can be useful for specimens with known orientation
ranges, but incorrect restrictions can exclude valid views and bias the
reconstruction.
## Local Alignment Options
In local alignment mode, the protocol can optimize several parameters around
the current particle assignment.
These include:
- shifts;
- scale;
- angles;
- gray-value scale and offset;
- defocus.
Shift and angle optimization are commonly useful. Scale, gray-value, and
defocus optimization are more specialized and should be used only when there is
a clear reason.
The local alignment also uses a Fourier padding factor, which affects
projection accuracy and computational cost.
## Restrict Reconstruction Angles
The **Restrict reconstruction angles** option allows reconstruction to use only
particles whose assigned tilt angles fall within a selected range.
This can be useful for special cases such as helices, where side views close to
90 degrees may be preferred for reconstruction.
This option should be used carefully. Excluding orientations can improve a
specific reconstruction strategy, but it also reduces the number of particles
and angular coverage.
## Particle Weights
The protocol can assign weights to particles during reconstruction.
Several weighting options are available:
**Weight by SSNR** uses signal-to-noise information.
**Weight by Continuous cost** uses the local-assignment cost.
**Weight by angular stability** downweights particles whose angular assignment
is unstable between iterations.
**Weight by CC percentile** weights particles according to their cross-
correlation percentile within their defocus group.
The **Minimum CC weight** defines the lower bound for CC-based weights.
These weights help reduce the influence of poorly fitting or unstable
particles without necessarily removing them completely.
## Particle Qualification
The protocol evaluates particle quality within defocus groups.
For global alignment, the relevant score is usually cross-correlation. For
local alignment, the relevant score may be a cost function. The protocol
computes percentiles within defocus groups so that particles are compared with
others acquired under similar CTF conditions.
This information is stored in the output particle metadata and can be useful
for later particle filtering or inspection.
## Reconstruction and Half Maps
Each iteration reconstructs two half maps from two particle subsets.
The two half maps are used to estimate FSC and monitor resolution. The final
average map for the iteration is produced from the two half maps.
The final output volume keeps the half-map file names associated with it, so
they can be used for validation or post-processing.
## Post-Processing Options
The protocol includes several optional post-processing operations.
These include:
- applying an ad hoc mask;
- symmetrizing within a mask;
- applying helical symmetry;
- running a custom post-processing command;
- significant denoising in real space;
- significant denoising in Fourier space;
- Laplacian denoising;
- blind deconvolution;
- attenuation of undershooting;
- difference evaluation.
These options are advanced and should be used only when the user understands
their purpose. Post-processing can improve interpretability, but inappropriate
post-processing can also introduce misleading features.
## Helical Symmetry Post-Processing
The protocol includes specialized options for applying helical symmetry during
post-processing.
The user can define a helical radius, request dihedral symmetry, and specify
search ranges for helical rotation and axial shift.
These options are useful only for helical specimens. They should not be enabled
for non-helical particles.
## Output Volume
The main output is **outputVolume**.
This volume corresponds to the final averaged reconstruction from the last
iteration. It is registered with the sampling rate used in the final iteration
and includes links to the two half maps.
The output should be inspected together with the FSC, filtered maps, and
particle statistics generated during refinement.
## Output Particles
The protocol also produces **outputParticles** when final angle metadata are
available.
This particle set contains the input particles updated with final projection
alignment parameters. The output particles may also include Xmipp attributes
such as shifts, angles, scale, maximum cross-correlation, percentile scores,
weights, angular stability, continuous-assignment cost, and related quality
fields.
This output is useful for later reconstruction, particle filtering, local
defocus refinement, or additional analysis.
## Practical Recommendations
Use this protocol after obtaining a reasonable initial model and a reasonably
clean particle set.
Start with global alignment if orientations are still uncertain. Move to local
alignment when particles stop changing substantially between iterations.
Use automatic alignment when you want the protocol to perform this transition
for you.
Keep the multiresolution approach enabled for most workflows. It reduces cost
and helps avoid using high-frequency information too early.
Use conservative low-pass filtering for the next reference to reduce
overfitting.
Inspect FSC curves, filtered volumes, particle weights, and the number of used
and rejected particles after each iteration.
Do not rely on this protocol to solve strong heterogeneity. If the particle
population is mixed, use classification or subset selection before high-
resolution refinement.
Use post-processing options cautiously and document them clearly if the output
map will be interpreted biologically.
## Final Perspective
Highres is an advanced iterative 3D refinement protocol for producing
high-resolution cryo-EM reconstructions from full-size particles.
For biological users, the protocol is most useful when the project has already
passed the initial-model and cleaning stages and the goal is to refine a mostly
homogeneous particle set toward a reliable final map.
The quality of the result depends on the starting reference, particle
homogeneity, angular-assignment strategy, CTF information, weighting choices,
and careful control of overfitting through half-map FSC evaluation and
reference filtering.
"""
_label = 'highres'
_devStatus = PROD
_lastUpdateVersion = VERSION_1_1
SPLIT_STOCHASTIC = 0
SPLIT_FIXED = 1
GLOBAL_ALIGNMENT = 0
LOCAL_ALIGNMENT = 1
AUTOMATIC_ALIGNMENT = 2
STOCHASTIC_ALIGNMENT = 3
NO_ALIGNMENT = 4
# --------------------------- DEFINE param functions --------------------------------------------
def _defineParams(self, form):
form.addHidden(USE_GPU, BooleanParam, default=True,
label="Use GPU for execution",
help="This protocol has both CPU and GPU implementation.\
Select the one you want to use.")
form.addHidden(GPU_LIST, StringParam, default='0',
expertLevel=LEVEL_ADVANCED,
label="Choose GPU IDs",
help="Add a list of GPU devices that can be used")
form.addSection(label='Input')
form.addParam('doContinue', BooleanParam, default=False,
label='Continue from a previous run?',
help='If you set to *Yes*, you should select a previous'
'run of type *%s* class and some of the input parameters'
'will be taken from it.' % self.getClassName())
form.addParam('inputParticles', PointerParam, label="Full-size Images", important=True,
pointerClass='SetOfParticles', allowsNull=True,
help='Select a set of images at full resolution')
form.addParam('inputVolumes', PointerParam, label="Initial volumes", allowsNull=True,
condition='not doContinue', pointerClass='Volume, SetOfVolumes',
help='Select a set of volumes with 2 volumes or a single volume. '
'If the input particles have an angular assignment, then you may '
'leave empty this field and a 3D reconstruction of the input images is '
'performed using reconstruct_fourier.')
form.addParam('particleRadius', IntParam, default=-1,
condition='not doContinue', label='Radius of particle (px)',
help='This is the radius (in pixels) of the spherical mask covering the particle in the input images')
form.addParam('continueRun', PointerParam, pointerClass=self.getClassName(),
condition='doContinue', allowsNull=True,
label='Select previous run',
help='Select a previous run to continue from.')
form.addParam('symmetryGroup', StringParam, default="c1",
label='Symmetry group',
help='If no symmetry is present, give c1')
form.addParam('correctEnvelope', BooleanParam, default=False,
label='Correct CTF envelope')
form.addParam("saveSpace", BooleanParam, default=True, label="Remove intermediate files", expertLevel=LEVEL_ADVANCED)
form.addSection(label='Next Reference')
form.addParam('nextLowPass', BooleanParam, label="Low pass filter?", default=True, expertLevel=LEVEL_ADVANCED,
help='Apply a low pass filter to the previous iteration whose maximum frequency is '\
'the current resolution(A) + resolutionOffset(A). If resolutionOffset>0, then fewer information' \
'is used (meant to avoid overfitting). If resolutionOffset<0, then more information is allowed '\
'(meant for a greedy convergence).')
form.addParam('nextResolutionCriterion',FloatParam, label="FSC criterion", default=0.5, expertLevel=LEVEL_ADVANCED,
help='The resolution of the reconstruction is defined as the inverse of the frequency at which '\
'the FSC drops below this value. Typical values are 0.143 and 0.5')
form.addParam('nextResolutionOffset', FloatParam, label="Resolution offset (A)", default=2, expertLevel=LEVEL_ADVANCED, condition='nextLowPass')
form.addParam('nextSpherical', BooleanParam, label="Spherical mask?", default=True, expertLevel=LEVEL_ADVANCED,
help='Apply a spherical mask of the size of the particle. If the postprocessing indicates that it has helical symmetry,'
'then a cylindrical mask is applied')
form.addParam('nextPositivity', BooleanParam, label="Positivity?", default=True, expertLevel=LEVEL_ADVANCED,
help='Remove from the next reference all negative values')
form.addParam('nextMask', PointerParam, label="Mask", pointerClass='VolumeMask', allowsNull=True,
help='The mask values must be between 0 (remove these pixels) and 1 (let them pass). Smooth masks are recommended.')
form.addParam('nextDropout', FloatParam, label="Dropout", default=0.0, expertLevel=LEVEL_ADVANCED,
help='This is the probability with which voxels are dropped (set to 0.0) inside the binary mask')
form.addParam('nextReferenceScript', StringParam, label="Next reference command", default="", expertLevel=LEVEL_ADVANCED,
help='A command template that is used to generate next reference. The following variables can be used '
'%(sampling)s %(dim)s %(volume)s %(iterDir)s. The command should read Spider volumes and modify the input volume.'
'the command should be accessible either from the PATH or provide the absolute path.\n'
'Examples: \n'
'xmipp_transform_filter -i %(volume)s --fourier low_pass 15 --sampling %(sampling)s\n'
'/home/joe/myScript %(volume)s sampling=%(sampling)s dim=%(dim)s')
form.addParam('nextRemove', BooleanParam, label="Remove reference to save space?", default=True, expertLevel=LEVEL_ADVANCED,
help='Remove reference volumes once they are not needed any more.')
form.addSection(label='Angular assignment')
form.addHidden('splitMethod', EnumParam, label='Image split method', choices=['Stochastic','Fixed'], default=self.SPLIT_FIXED,
expertLevel=LEVEL_ADVANCED)
form.addParam('multiresolution', BooleanParam, label='Multiresolution approach', default=True, expertLevel=LEVEL_ADVANCED,
help="In the multiresolution approach the sampling rate of the images is adapted to the current resolution")
form.addParam('angularMaxShift', FloatParam, label="Max. shift (%)", default=10, expertLevel=LEVEL_ADVANCED,
help='Maximum shift as a percentage of the image size')
line=form.addLine('Tilt angle:', help='0 degrees represent top views, 90 degrees represent side views', expertLevel=LEVEL_ADVANCED)
line.addParam('angularMinTilt', FloatParam, label="Min.", default=0, expertLevel=LEVEL_ADVANCED,
help="Side views are around 90 degrees, top views around 0")
line.addParam('angularMaxTilt', FloatParam, label="Max.", default=180, expertLevel=LEVEL_ADVANCED,
help="You may generate redudant galleries by setting this angle to 180, this may help if c1 symmetry is considered")
form.addParam('alignmentMethod', EnumParam, label='Image alignment', choices=['Global','Local','Automatic','Stochastic','No alignment'],
default=self.GLOBAL_ALIGNMENT)
form.addParam('numberOfIterations', IntParam, default=3, label='Number of iterations', condition='alignmentMethod!=2 and alignmentMethod!=4')
form.addParam('NimgsSGD', IntParam, default=250, label='Random subset size', condition='alignmentMethod==3',
expertLevel=LEVEL_ADVANCED, help="Stochastic alignment is performed by taking random subsets of images of this size")
form.addParam('alphaSGD', FloatParam, default=0.1, label='Step size', condition='alignmentMethod==3',
expertLevel=LEVEL_ADVANCED, help="The update is performed as V(k+1)=(1-alpha)*V(k)+alpha*R(k+1), that is, the previous "
"volume weights 1-alpha, while the new one weights alpha")
form.addParam("restrictReconstructionAngles", BooleanParam, label="Restrict reconstruction angles", default=False, expertLevel=LEVEL_ADVANCED,
help="You may reconstruct only with those images falling on a certain range. This is particularly useful for helices where "\
"you may want to use projections very close to 90 degrees")
line=form.addLine('Tilt angle restriction:', help='0 degrees represent top views, 90 degrees represent side views',
condition="restrictReconstructionAngles")
line.addParam('angularMinTiltReconstruct', FloatParam, label="Min.", default=88, condition="restrictReconstructionAngles",
help = "Perform an angular assignment and only use those images whose angles are within these limits")
line.addParam('angularMaxTiltReconstruct', FloatParam, label="Max.", default=92, condition="restrictReconstructionAngles",
help = "Perform an angular assignment and only use those images whose angles are within these limits")
form.addParam('maximumTargetResolution', NumericListParam,
label="Max. Target Resolution", default="15 8 4",
condition='multiresolution',
help="In Angstroms. The actual maximum resolution will be the maximum between this number of 0.5 * previousResolution, meaning that"
"in a single step you cannot increase the resolution more than 1/2")
form.addHidden('numberOfPerturbations', IntParam, label="Number of Perturbations", default=1, condition='alignmentMethod!=1',
expertLevel=LEVEL_ADVANCED, help="The gallery of reprojections is randomly perturbed this number of times")
form.addHidden('numberOfReplicates', IntParam, label="Max. Number of Replicates", default=1, condition='alignmentMethod!=1',
expertLevel=LEVEL_ADVANCED, help="Significant alignment is allowed to replicate each image up to this number of times")
form.addParam('contShift', BooleanParam, label="Optimize shifts?", default=True, condition='alignmentMethod==1',
help='Optimize shifts within a limit')
form.addParam('contMaxShiftVariation', FloatParam, label="Max. shift variation", default=2, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED,
help="Percentage of the image size")
form.addParam('contScale', BooleanParam, label="Optimize scale?", default=False, condition='alignmentMethod==1',
help='Optimize scale within a limit')
form.addParam('contMaxScale', FloatParam, label="Max. scale variation", default=0.02, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED)
form.addParam('contAngles', BooleanParam, label="Optimize angles?", default=True, condition='alignmentMethod==1',
help='Optimize angles within a limit')
form.addParam('contGrayValues', BooleanParam, label="Optimize gray values?", default=False, condition='alignmentMethod==1',
help='Optimize gray values. Do not perform this unless the reconstructed volume is gray-compatible with the projections,'\
' i.e., the volumes haven been produced from projections')
form.addParam('contMaxGrayScale', FloatParam, label="Max. gray scale variation", default=0.5, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED)
form.addParam('contMaxGrayShift', FloatParam, label="Max. gray shift variation", default=3, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED,
help='As a factor of the image standard deviation')
form.addParam('contDefocus', BooleanParam, label="Optimize defocus?", condition='alignmentMethod==1', default=False)
form.addParam('contMaxDefocus', FloatParam, label="Max. defocus variation", default=200, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED,
help="In Angstroms")
form.addParam('contPadding', IntParam, label="Fourier padding factor", default=2, condition='alignmentMethod==1', expertLevel=LEVEL_ADVANCED,
help='The volume is zero padded by this factor to produce projections')
form.addSection(label='Weights')
form.addParam('weightSSNR', BooleanParam, label="Weight by SSNR?", default=False, expertLevel=LEVEL_ADVANCED,
help='Weight input images by SSNR')
form.addParam('weightContinuous', BooleanParam, label="Weight by Continuous cost?", default=False, expertLevel=LEVEL_ADVANCED,
condition='alignmentMethod==1', help='Weight input images by angular assignment cost')
form.addParam('weightJumper', BooleanParam, label="Weight by angular stability?", default=False, expertLevel=LEVEL_ADVANCED,
help='Weight input images by angular stability between iterations')
form.addParam('weightCC', BooleanParam, label="Weight by CC percentile?", default=True, expertLevel=LEVEL_ADVANCED,
help='Weight input images by their fitness (cross correlation) percentile in their defocus group')
form.addParam('weightCCmin', FloatParam, label="Minimum CC weight", default=0.1, expertLevel=LEVEL_ADVANCED,
help='Weights are between this value and 1. If most of the particles are good, this value should be high (e.g., 0.9)')
form.addSection(label='Post-processing')
form.addParam('postAdHocMask', PointerParam, label="Mask", pointerClass='VolumeMask', allowsNull=True,
help='The mask values must be between 0 (remove these pixels) and 1 (let them pass). Smooth masks are recommended.')
groupSymmetry = form.addGroup('Symmetry', expertLevel=LEVEL_ADVANCED)
groupSymmetry.addParam('postSymmetryWithinMask', BooleanParam, label="Symmetrize volume within mask?", default=False)
groupSymmetry.addParam('postSymmetryWithinMaskType', StringParam, label="Mask symmetry", default="i1", condition="postSymmetryWithinMask",
help='If no symmetry is present, give c1')
groupSymmetry.addParam('postSymmetryWithinMaskMask', PointerParam, label="Mask", pointerClass='VolumeMask', allowsNull=True, condition="postSymmetryWithinMask",
help='The mask values must be between 0 (remove these pixels) and 1 (let them pass). Smooth masks are recommended.')
groupSymmetry.addParam('postSymmetryHelical', BooleanParam, label="Apply helical symmetry?", default=False)
groupSymmetry.addParam('postSymmetryHelicalRadius', IntParam, label="Radius", default=-1, condition='postSymmetryHelical',
help="In Angstroms")
groupSymmetry.addParam('postSymmetryHelicalDihedral', BooleanParam, label="Dihedral symmetry", default=False,
condition='postSymmetryHelical')
groupSymmetry.addParam('postSymmetryHelicalMinRot', FloatParam, label="Min. Rotation", default=0, condition='postSymmetryHelical',
help="In degrees")
groupSymmetry.addParam('postSymmetryHelicalMaxRot', FloatParam, label="Max. Rotation", default=360, condition='postSymmetryHelical',
help="In degrees")
groupSymmetry.addParam('postSymmetryHelicalMinZ', FloatParam, label="Min. Z shift", default=0, condition='postSymmetryHelical',
help="In angstroms")
groupSymmetry.addParam('postSymmetryHelicalMaxZ', FloatParam, label="Max. Z shift", default=40, condition='postSymmetryHelical',
help="In angstroms")
form.addParam('postScript', StringParam, label="Post-processing command", default="", expertLevel=LEVEL_ADVANCED,
help='A command template that is used to post-process the reconstruction. The following variables can be used '
'%(sampling)s %(dim)s %(volume)s %(iterDir)s. The command should read Spider volumes and modify the input volume.'
'the command should be accessible either from the PATH or provide the absolute path.\n'
'Examples: \n'
'xmipp_transform_filter -i %(volume)s --fourier low_pass 15 --sampling %(sampling)s\n'
'/home/joe/myScript %(volume)s sampling=%(sampling)s dim=%(dim)s')
form.addParam('postSignificantDenoise', BooleanParam, label="Significant denoising Real space", expertLevel=LEVEL_ADVANCED, default=False)
form.addParam('postFilterBank', BooleanParam, label="Significant denoising Fourier space", expertLevel=LEVEL_ADVANCED, default=False)
form.addParam('postLaplacian', BooleanParam, label="Laplacian denoising", expertLevel=LEVEL_ADVANCED, default=False,
help="It can only be used if there is a mask")
form.addParam('postDeconvolve', BooleanParam, label="Blind deconvolution", expertLevel=LEVEL_ADVANCED, default=False)
form.addParam('postSoftNeg', BooleanParam, label="Attenuate undershooting", expertLevel=LEVEL_ADVANCED, default=False)
form.addParam('postSoftNegK', FloatParam, label="Attenuate undershooting (K)", expertLevel=LEVEL_ADVANCED, default=9,
help="Values below avg-K*sigma are attenuated")
form.addParam('postDifference', BooleanParam, label="Evaluate difference", expertLevel=LEVEL_ADVANCED, default=False)
form.addParallelSection(threads=1, mpi=8)
#--------------------------- INSERT steps functions --------------------------------------------
def _insertAllSteps(self):
self.imgsFn=self._getExtraPath('images.xmd')
if self.doContinue:
self.copyAttributes(self.continueRun.get(), 'particleRadius')
self.copyAttributes(self.continueRun.get(), 'inputVolumes')
if not self.inputParticles.hasValue():
self.copyAttributes(self.continueRun.get(), 'inputParticles')
else:
self._insertFunctionStep('convertInputStep', self.inputParticles.getObjId())
self._insertFunctionStep('copyBasicInformation')
self.firstIteration=self.getNumberOfPreviousIterations()+1
else:
self._insertFunctionStep('convertInputStep', self.inputParticles.getObjId())
if self.weightSSNR:
self._insertFunctionStep('doWeightSSNR')
self._insertFunctionStep('doIteration000')
self.firstIteration=1
self.TsOrig=self.inputParticles.get().getSamplingRate()
numberOfIterations = self.numberOfIterations.get() if self.alignmentMethod.get()!=self.AUTOMATIC_ALIGNMENT else 5
if self.alignmentMethod.get()==self.NO_ALIGNMENT:
numberOfIterations = 1
self._maximumTargetResolution = getFloatListFromValues(self.maximumTargetResolution.get(),self.firstIteration+numberOfIterations-1)
for self.iteration in range(self.firstIteration,self.firstIteration+numberOfIterations):
self.insertIteration(self.iteration)
self._insertFunctionStep("createOutput")
[docs] def insertIteration(self,iteration):
if self.alignmentMethod==self.GLOBAL_ALIGNMENT or \
self.alignmentMethod==self.STOCHASTIC_ALIGNMENT or \
(self.alignmentMethod==self.AUTOMATIC_ALIGNMENT and iteration<=3):
self._insertFunctionStep('globalAssignment',iteration)
elif self.alignmentMethod==self.NO_ALIGNMENT:
self._insertFunctionStep('noAssignment', iteration)
else:
self._insertFunctionStep('localAssignment',iteration)
self._insertFunctionStep('weightParticles',iteration)
self._insertFunctionStep('qualifyParticles',iteration)
self._insertFunctionStep('reconstruct',iteration)
self._insertFunctionStep('postProcessing',iteration)
self._insertFunctionStep('evaluateReconstructions',iteration)
self._insertFunctionStep('cleanDirectory',iteration)
#--------------------------- STEPS functions ---------------------------------------------------
[docs] def getNumGpus(self):
return len(self._stepsExecutor.getGpuList())
[docs] def getGpusList(self, separator):
strGpus = ""
for elem in self._stepsExecutor.getGpuList():
strGpus = strGpus + str(elem) + separator
return strGpus[:-1]
[docs] def setGPU(self, oneGPU=False):
if oneGPU:
gpus = self.getGpusList(",")[0]
else:
gpus = self.getGpusList(",")
os.environ["CUDA_VISIBLE_DEVICES"] = gpus
self.info(f'Visible GPUS: {gpus}')
return gpus
[docs] def createOutput(self):
# get last iteration
fnIterDir=glob(self._getExtraPath("Iter*"))
lastIter=len(fnIterDir)-1
# fnFirstDir=self._getExtraPath("Iter000")
fnLastDir=self._getExtraPath("Iter%03d"%lastIter)
# fnFirstVol=join(fnFirstDir,"volume000.mrc")
fnLastVol=join(fnLastDir,"volumeAvg.mrc")
Ts=self.readInfoField(fnLastDir,"sampling",emlib.MDL_SAMPLINGRATE)
if exists(fnLastVol):
volume=Volume()
volume.setFileName(fnLastVol)
volume.setSamplingRate(Ts)
halfMap1=join(fnLastDir,"volume01.vol")
halfMap2=join(fnLastDir,"volume02.vol")
volume.setHalfMaps([halfMap1, halfMap2])
self._defineOutputs(outputVolume=volume)
self._defineSourceRelation(self.inputParticles.get(),volume)
if not self.doContinue and self.inputVolumes.get() is not None:
self._defineSourceRelation(self.inputVolumes.get(),volume)
fnLastAngles=join(fnLastDir,"angles.xmd")
if exists(fnLastAngles):
fnAngles=self._getPath("angles.xmd")
self.runJob('xmipp_metadata_utilities','-i %s -o %s --operate modify_values "image=image1"'%(fnLastAngles,fnAngles),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities','-i %s --operate sort particleId'%fnAngles,numberOfMpi=1)
self.runJob('xmipp_metadata_utilities','-i %s --operate drop_column image1'%fnAngles,numberOfMpi=1)
self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "itemId=particleId"'%fnAngles,numberOfMpi=1)
imgSet = self.inputParticles.get()
self.scaleFactor=Ts/imgSet.getSamplingRate()
imgSetOut = self._createSetOfParticles()
imgSetOut.copyInfo(imgSet)
imgSetOut.setAlignmentProj()
imgSetOut.setIsPhaseFlipped( imgSet.isPhaseFlipped() )
self.iterMd = md.iterRows(fnAngles, md.MDL_PARTICLE_ID)
self.lastRow = next(self.iterMd)
imgSetOut.copyItems(imgSet,
updateItemCallback=self._updateItem)
self._defineOutputs(outputParticles=imgSetOut)
self._defineSourceRelation(self.inputParticles, imgSetOut)
def _updateItem(self, particle, row):
count = 0
while self.lastRow and particle.getObjId() == self.lastRow.getValue(md.MDL_PARTICLE_ID):
count += 1
if count:
self._createItemMatrix(particle, self.lastRow)
try:
self.lastRow = next(self.iterMd)
except StopIteration:
self.lastRow = None
particle._appendItem = count > 0
def _createItemMatrix(self, particle, row):
if row.containsLabel(emlib.MDL_CONTINUOUS_X):
row.setValue(emlib.MDL_SHIFT_X, row.getValue(emlib.MDL_CONTINUOUS_X))
row.setValue(emlib.MDL_SHIFT_Y, row.getValue(emlib.MDL_CONTINUOUS_Y))
row.setValue(emlib.MDL_FLIP, row.getValue(emlib.MDL_CONTINUOUS_FLIP))
row.setValue(emlib.MDL_SHIFT_X, row.getValue(emlib.MDL_SHIFT_X)*self.scaleFactor)
row.setValue(emlib.MDL_SHIFT_Y, row.getValue(emlib.MDL_SHIFT_Y)*self.scaleFactor)
setXmippAttributes(particle, row, emlib.MDL_SHIFT_X, emlib.MDL_SHIFT_Y, emlib.MDL_ANGLE_TILT,
emlib.MDL_ANGLE_ROT,
emlib.MDL_SCALE, emlib.MDL_MAXCC, emlib.MDL_MAXCC_PERCENTILE, emlib.MDL_WEIGHT)
if row.containsLabel(emlib.MDL_ANGLE_DIFF0):
setXmippAttributes(particle, row, emlib.MDL_ANGLE_DIFF0, emlib.MDL_WEIGHT_JUMPER0)
if row.containsLabel(emlib.MDL_CONTINUOUS_X):
setXmippAttributes(particle, row, emlib.MDL_COST, emlib.MDL_WEIGHT_CONTINUOUS2, emlib.MDL_COST_PERCENTILE,
emlib.MDL_CORRELATION_IDX,
emlib.MDL_CORRELATION_MASK,
emlib.MDL_CORRELATION_WEIGHT,
emlib.MDL_IMED)
if row.containsLabel(emlib.MDL_CONTINUOUS_SCALE_X):
setXmippAttributes(emlib.MDL_CONTINUOUS_SCALE_X, emlib.MDL_CONTINUOUS_SCALE_Y)
if row.containsLabel(emlib.MDL_CONTINUOUS_GRAY_A):
setXmippAttributes(emlib.MDL_CONTINUOUS_GRAY_A, emlib.MDL_CONTINUOUS_GRAY_B)
if row.containsLabel(emlib.MDL_WEIGHT_JUMPER):
setXmippAttributes(particle, row, emlib.MDL_WEIGHT_JUMPER)
if row.containsLabel(emlib.MDL_ANGLE_DIFF):
setXmippAttributes(particle, row, emlib.MDL_ANGLE_DIFF)
if row.containsLabel(emlib.MDL_ANGLE_DIFF2):
setXmippAttributes(particle, row, emlib.MDL_ANGLE_DIFF2)
if row.containsLabel(emlib.MDL_ANGLE_TEMPERATURE):
setXmippAttributes(particle, row, emlib.MDL_ANGLE_TEMPERATURE)
if row.containsLabel(emlib.MDL_WEIGHT_SSNR):
setXmippAttributes(particle, row, emlib.MDL_WEIGHT_SSNR)
createItemMatrix(particle, row, align=ALIGN_PROJ)
[docs] def getLastFinishedIter(self):
fnFscs=sorted(glob(self._getExtraPath("Iter???/fsc.xmd")))
lastDir=split(fnFscs[-1])[0]
return int(lastDir[-3:])
[docs] def getNumberOfPreviousIterations(self):
fnDirs=sorted(glob(self.continueRun.get()._getExtraPath("Iter???")))
lastDir=fnDirs[-1]
return int(lastDir[-3:])
[docs] def doWeightSSNR(self):
R=self.particleRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2
self.runJob("xmipp_image_ssnr", "-i %s -R %d --sampling %f --normalizessnr"%\
(self.imgsFn,R,self.inputParticles.get().getSamplingRate()),numberOfMpi=min(self.numberOfMpi.get(),24))
self.runJob('xmipp_metadata_utilities','-i %s -o %s --operate keep_column "particleId weightSSNR" '%\
(self.imgsFn,self._getExtraPath("ssnrWeights.xmd")),numberOfMpi=1)
[docs] def doIteration000(self):
fnDirCurrent=self._getExtraPath('Iter000')
makePath(fnDirCurrent)
# Split data
if self.splitMethod == self.SPLIT_FIXED:
self.runJob("xmipp_metadata_split","-i %s --oroot %s/images -n 2"%(self.imgsFn,fnDirCurrent),numberOfMpi=1)
for i in range(1,3):
moveFile("%s/images%06d.xmd"%(fnDirCurrent,i),"%s/images%02d.xmd"%(fnDirCurrent,i))
# Get volume sampling rate
if self.inputVolumes.get() is None:
TsCurrent=self.inputParticles.get().getSamplingRate()
else:
TsCurrent=self.inputVolumes.get().getSamplingRate()
self.writeInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE,TsCurrent)
# Copy reference volumes and window if necessary
Xdim=self.inputParticles.get().getDimensions()[0]
newXdim=int(round(Xdim*self.TsOrig/TsCurrent))
self.writeInfoField(fnDirCurrent,"size",emlib.MDL_XSIZE,newXdim)
img = ImageHandler()
if isinstance(self.inputVolumes.get(),SetOfVolumes):
i=1
for vol in self.inputVolumes.get():
fnVol=join(fnDirCurrent,"volume%02d.vol"%i)
if i==1:
fnVol1 = fnVol
else:
fnVol2 = fnVol
img.convert(vol, fnVol)
if newXdim!=vol.getDim()[0]:
self.runJob('xmipp_transform_window',"-i %s --size %d"%(fnVol,newXdim),numberOfMpi=1)
i+=1
else:
fnVol1=join(fnDirCurrent,"volume%02d.vol"%1)
fnVol2=join(fnDirCurrent,"volume%02d.vol"%2)
if self.inputVolumes.get() is None:
args = "-i %s -o %s --max_resolution 0.3 --sampling %f --sym %s" % (
self.imgsFn, fnVol1, TsCurrent, self.symmetryGroup.get())
if self.useGpu.get():
#AJ to make it work with and without queue system
if self.numberOfMpi.get()>1:
gpuId = self.setGPU(oneGPU=False)
args += ' -gpusPerNode %d' % self.getNumGpus()
args += ' -threadsPerGPU %d' % max(self.numberOfThreads.get(),4)
if self.numberOfMpi.get()==1:
gpuId = self.setGPU(oneGPU=True)
args += " --device %s" % (gpuId.replace(',', ' '))
args += ' --thr %s' % self.numberOfThreads.get()
if self.numberOfMpi.get()>1:
self.runJob('xmipp_cuda_reconstruct_fourier', args, numberOfMpi=self.getNumGpus()+1)
else:
self.runJob('xmipp_cuda_reconstruct_fourier', args)
else:
self.runJob("xmipp_reconstruct_fourier_accel", args,
numberOfMpi=self.numberOfMpi.get())
volXdim = Xdim
else:
vol=self.inputVolumes.get()
img.convert(vol, fnVol1)
volXdim = vol.getDim()[0]
if newXdim!=volXdim:
self.runJob('xmipp_transform_window',"-i %s --size %d"%(fnVol1,newXdim),numberOfMpi=1)
if self.alignmentMethod != self.LOCAL_ALIGNMENT:
maxFreq=0.25
else:
maxFreq=0.45
self.runJob('xmipp_transform_randomize_phases',"-i %s -o %s --freq discrete %f"%(fnVol1,fnVol2,maxFreq),numberOfMpi=1)
# Compare both reconstructions
self.evaluateReconstructions(0)
# Get the angles and shifts from images into this directory as if this directory
# had been the result of a global alignment
if self.alignmentMethod.get()==self.LOCAL_ALIGNMENT:
copyFile(self.imgsFn,join(fnDirCurrent,"angles.xmd"))
[docs] def evaluateReconstructions(self,iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
fnVol1=join(fnDirCurrent,"volume%02d.vol"%1)
fnVol2=join(fnDirCurrent,"volume%02d.vol"%2)
fnVolFSC1=join(fnDirCurrent,"volumeFSC%02d.vol"%1)
fnVolFSC2=join(fnDirCurrent,"volumeFSC%02d.vol"%2)
TsCurrent=self.readInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE)
if not exists(fnVolFSC1):
copyFile(fnVol1,fnVolFSC1)
copyFile(fnVol2,fnVolFSC2)
# Apply mask if available
fnMask=""
volXdim = self.readInfoField(fnDirCurrent, "size", emlib.MDL_XSIZE)
if self.postAdHocMask.hasValue():
fnMask=join(fnDirCurrent,"mask.vol")
if not exists(fnMask):
self.prepareMask(self.postAdHocMask.get(), fnMask, TsCurrent, volXdim)
self.runJob("xmipp_image_operate","-i %s --mult %s"%(fnVolFSC1,fnMask),numberOfMpi=1)
self.runJob("xmipp_image_operate","-i %s --mult %s"%(fnVolFSC2,fnMask),numberOfMpi=1)
# Threshold
self.runJob('xmipp_transform_threshold','-i %s --select below 0 --substitute value 0 '%fnVolFSC1,numberOfMpi=1)
self.runJob('xmipp_transform_threshold','-i %s --select below 0 --substitute value 0 '%fnVolFSC2,numberOfMpi=1)
# Estimate resolution
fnFsc=join(fnDirCurrent,"fsc.xmd")
self.runJob('xmipp_resolution_fsc','--ref %s -i %s -o %s --sampling_rate %f'\
%(fnVolFSC1,fnVolFSC2,fnFsc,TsCurrent),numberOfMpi=1)
cleanPath(fnVolFSC1)
cleanPath(fnVolFSC2)
if fnMask!="":
cleanPath(fnMask)
# Estimate resolution before postprocessing
fnBeforeVol1=join(fnDirCurrent,"volumeBeforePostProcessing%02d.vol"%1)
fnBeforeVol2=join(fnDirCurrent,"volumeBeforePostProcessing%02d.vol"%2)
if exists(fnBeforeVol1) and exists(fnBeforeVol2):
fnBeforeFsc=join(fnDirCurrent,"fscBeforePostProcessing.xmd")
self.runJob('xmipp_resolution_fsc','--ref %s -i %s -o %s --sampling_rate %f'%(fnBeforeVol1,fnBeforeVol2,fnBeforeFsc,TsCurrent),
numberOfMpi=1)
mdFSC = emlib.MetaData(fnFsc)
resolution=2*TsCurrent
for objId in mdFSC:
fsc = mdFSC.getValue(emlib.MDL_RESOLUTION_FRC,objId)
if fsc<self.nextResolutionCriterion.get():
resolution=mdFSC.getValue(emlib.MDL_RESOLUTION_FREQREAL,objId)
break
if iteration==0:
if self.alignmentMethod != self.LOCAL_ALIGNMENT:
resolution = TsCurrent * 1/0.25
else:
resolution = TsCurrent * 1/0.45
self.writeInfoField(fnDirCurrent,"resolution",emlib.MDL_RESOLUTION_FREQREAL,resolution)
# Produce a filtered volume
if iteration>0:
self.runJob('xmipp_transform_filter','-i %s -o %s --fourier low_pass %f --sampling %f'%\
(join(fnDirCurrent,"volumeAvg.mrc"),join(fnDirCurrent,"volumeAvgFiltered.mrc"),resolution,TsCurrent),numberOfMpi=1)
# A little bit of statistics (accepted and rejected particles, number of directions, ...)
if iteration>0:
for i in range(1,3):
fnAnglesi = join(fnDirCurrent,"angles%02d.xmd"%i)
mdAngles = emlib.MetaData(fnAnglesi)
mdUnique = emlib.MetaData()
mdUnique.aggregateMdGroupBy(mdAngles, emlib.AGGR_MAX, [emlib.MDL_PARTICLE_ID], emlib.MDL_WEIGHT, emlib.MDL_WEIGHT)
mdUnique.sort(emlib.MDL_PARTICLE_ID)
fnAnglesUnique = join(fnDirCurrent,"imagesUsed%02d.xmd"%i)
mdUnique.write(fnAnglesUnique)
fnUsed=join(fnDirCurrent,"imagesUsed.xmd")
fnUsed1=join(fnDirCurrent,"imagesUsed01.xmd")
fnUsed2=join(fnDirCurrent,"imagesUsed02.xmd")
self.runJob('xmipp_metadata_utilities',"-i %s --set union_all %s -o %s"%(fnUsed1,fnUsed2,fnUsed),numberOfMpi=1)
cleanPath(fnUsed1)
cleanPath(fnUsed2)
fnAngles=join(fnDirCurrent,"angles.xmd")
fnUsedId=join(fnDirCurrent,"imagesUsedId.xmd")
self.runJob('xmipp_metadata_utilities',"-i %s --operate keep_column particleId -o %s"%(fnUsed,fnUsedId),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities',"-i %s --set natural_join %s"%(fnUsedId,fnAngles),numberOfMpi=1)
fnImages=self._getExtraPath("images.xmd")
fnImagesId=self._getExtraPath('imagesId.xmd')
fnImagesRejected=join(fnDirCurrent,"imagesRejected.xmd")
self.runJob('xmipp_metadata_utilities',"-i %s --set subtraction %s particleId -o %s"%(fnImagesId,fnUsedId,fnImagesRejected),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities',"-i %s --set natural_join %s"%(fnImagesRejected,fnImages),numberOfMpi=1)
cleanPath(fnUsedId)
Nimages=getSize(fnImages)
Nrepeated=getSize(join(fnDirCurrent,"angles.xmd"))
Nunique=getSize(fnUsed)
Nrejected=getSize(fnImagesRejected)
fh=open(join(fnDirCurrent,"statistics.txt"),'w')
fh.write("Number of input images: %d\n"%Nimages)
fh.write("Number of used images: %d\n"%Nunique)
fh.write("Number of rejected images: %d\n"%Nrejected)
if Nunique>0:
fh.write("Average number of directions per used image: %f\n"%(float(Nrepeated)/Nunique))
fh.close()
[docs] def checkInfoField(self,fnDir,block):
fnInfo = join(fnDir,"iterInfo.xmd")
if not exists(fnInfo):
return False
blocks = emlib.getBlocksInMetaDataFile(fnInfo)
return block in blocks
[docs] def readInfoField(self,fnDir,block,label):
mdInfo = emlib.MetaData("%s@%s"%(block,join(fnDir,"iterInfo.xmd")))
return mdInfo.getValue(label,mdInfo.firstObject())
[docs] def writeInfoField(self,fnDir,block,label, value):
mdInfo = emlib.MetaData()
objId=mdInfo.addObject()
mdInfo.setValue(label,value,objId)
mdInfo.write("%s@%s"%(block,join(fnDir,"iterInfo.xmd")),emlib.MD_APPEND)
[docs] def prepareImages(self,fnDirPrevious,fnDir,TsCurrent,getShiftsFrom=''):
if self.checkInfoField(fnDir,"count"):
state = self.readInfoField(fnDir, "count", emlib.MDL_COUNT)
if state>=1:
return
print("Preparing images to sampling rate=",TsCurrent)
Xdim=self.inputParticles.get().getDimensions()[0]
newXdim=int(round(Xdim*self.TsOrig/TsCurrent))
if newXdim<40:
newXdim=int(40)
TsCurrent=Xdim*(self.TsOrig/newXdim)
elif newXdim%2==1:
newXdim+=1
TsCurrent=Xdim*(self.TsOrig/newXdim)
self.writeInfoField(fnDir,"sampling",emlib.MDL_SAMPLINGRATE,TsCurrent)
self.writeInfoField(fnDir,"size",emlib.MDL_XSIZE,newXdim)
# Prepare particles
fnDir0=self._getExtraPath("Iter000")
fnNewParticles=join(fnDir,"images.stk")
if newXdim!=Xdim:
self.runJob("xmipp_image_resize","-i %s -o %s --fourier %d"%(self.imgsFn,fnNewParticles,newXdim),numberOfMpi=min(self.numberOfMpi.get(),24))
else:
self.runJob("xmipp_image_convert","-i %s -o %s --save_metadata_stack %s"%(self.imgsFn,fnNewParticles,join(fnDir,"images.xmd")),
numberOfMpi=1)
R=self.particleRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2
R=min(round(R*self.TsOrig/TsCurrent*(1+self.angularMaxShift.get()*0.01)),newXdim/2)
self.runJob("xmipp_transform_mask","-i %s --mask circular -%d"%(fnNewParticles,R),numberOfMpi=min(self.numberOfMpi.get(),24))
fnSource=join(fnDir,"images.xmd")
if not self.inputParticles.get().isPhaseFlipped():
self.runJob("xmipp_ctf_correct_phase", "-i %s --sampling_rate %f" % (fnSource, TsCurrent),
numberOfMpi=min(self.numberOfMpi.get(), 24))
if self.splitMethod==self.SPLIT_STOCHASTIC:
self.runJob('xmipp_metadata_utilities','-i %s --set intersection %s particleId particleId -o %s/all_images.xmd'%\
(fnSource,self._getExtraPath('images.xmd'),fnDir),numberOfMpi=1)
self.runJob("xmipp_metadata_split","-i %s/all_images.xmd --oroot %s/images -n 2"%(fnDir,fnDir),numberOfMpi=1)
cleanPath("%s/all_images.xmd"%fnDir)
for i in range(1,3):
moveFile("%s/images%06d.xmd"%(fnDir,i),"%s/images%02d.xmd"%(fnDir,i))
else:
for i in range(1,3):
fnImagesi=join(fnDir,"images%02d.xmd"%i)
self.runJob('xmipp_metadata_utilities','-i %s --set intersection %s/images%02d.xmd particleId particleId -o %s'%\
(fnSource,fnDir0,i,fnImagesi),numberOfMpi=1)
cleanPath(fnSource)
if self.alignmentMethod==self.STOCHASTIC_ALIGNMENT:
for i in range(1,3):
fnImagesi=join(fnDir,"images%02d.xmd"%i)
self.runJob('xmipp_metadata_utilities','-i %s --operate random_subset %d'%\
(fnImagesi,self.NimgsSGD),numberOfMpi=1)
if getShiftsFrom!="":
fnPreviousAngles=join(getShiftsFrom,"angles.xmd")
TsPrevious=self.readInfoField(getShiftsFrom,"sampling",emlib.MDL_SAMPLINGRATE)
fnAux=join(fnDir,"aux.xmd")
for i in range(1,3):
fnImagesi=join(fnDir,"images%02d.xmd"%i)
self.runJob('xmipp_metadata_utilities','-i %s --set join %s particleId particleId -o %s'%\
(fnImagesi,fnPreviousAngles,fnAux),numberOfMpi=1)
self.adaptShifts(fnAux, TsPrevious, fnImagesi, TsCurrent)
cleanPath(fnAux)
self.writeInfoField(fnDir,"count",emlib.MDL_COUNT,int(1))
[docs] def prepareReferences(self,fnDirPrevious,fnDir,TsCurrent,targetResolution):
if self.checkInfoField(fnDir,"count"):
state = self.readInfoField(fnDir, "count", emlib.MDL_COUNT)
if state>=2:
return
print("Preparing references to sampling rate=",TsCurrent)
fnMask=''
newXdim=self.readInfoField(fnDir,"size",emlib.MDL_XSIZE)
if self.nextMask.hasValue():
fnMask=join(fnDir,"mask.vol")
self.prepareMask(self.nextMask.get(), fnMask, TsCurrent, newXdim)
oldXdim=self.readInfoField(fnDirPrevious,"size",emlib.MDL_XSIZE)
for i in range(1,3):
fnPreviousVol=join(fnDirPrevious,"volume%02d.vol"%i)
fnReferenceVol=join(fnDir,"volumeRef%02d.vol"%i)
if oldXdim!=newXdim:
self.runJob("xmipp_image_resize","-i %s -o %s --dim %d"%(fnPreviousVol,fnReferenceVol,newXdim),numberOfMpi=1)
else:
copyFile(fnPreviousVol, fnReferenceVol)
self.runJob('xmipp_transform_filter','-i %s --fourier fsc %s --sampling %f'%(fnReferenceVol,join(fnDirPrevious,"fsc.xmd"),TsCurrent),numberOfMpi=1)
if self.nextLowPass:
self.runJob('xmipp_transform_filter','-i %s --fourier low_pass %f --sampling %f'%\
(fnReferenceVol,targetResolution+self.nextResolutionOffset.get(),TsCurrent),numberOfMpi=1)
if self.nextSpherical:
if self.postSymmetryHelical:
R=self.postSymmetryHelicalRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2*self.TsOrig
self.runJob('xmipp_transform_mask','-i %s --mask cylinder -%d -%d'%\
(fnReferenceVol,round(R/TsCurrent),newXdim),numberOfMpi=1)
else:
R=self.particleRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2
self.runJob('xmipp_transform_mask','-i %s --mask circular -%d'%\
(fnReferenceVol,round(R*self.TsOrig/TsCurrent)),numberOfMpi=1)
if self.nextPositivity:
self.runJob('xmipp_transform_threshold','-i %s --select below 0 --substitute value 0'%fnReferenceVol,numberOfMpi=1)
if fnMask!='':
self.runJob('xmipp_image_operate','-i %s --mult %s'%(fnReferenceVol,fnMask),numberOfMpi=1)
if self.nextDropout.get()>0.0:
self.runJob('xmipp_image_operate','-i %s --dropout %f'%(fnReferenceVol,self.nextDropout),numberOfMpi=1)
if self.nextReferenceScript!="":
scriptArgs = {'volume': fnReferenceVol,
'sampling': TsCurrent,
'dim': newXdim,
'iterDir': fnDir}
cmd = self.nextReferenceScript % scriptArgs
self.runJob(cmd, '', numberOfMpi=1)
if fnMask!='':
cleanPath(fnMask)
self.writeInfoField(fnDir,"count",emlib.MDL_COUNT,int(2))
[docs] def prepareMask(self,maskObject,fnMask,TsMaskOut,XdimOut):
img=ImageHandler()
img.convert(maskObject, fnMask)
self.runJob('xmipp_image_resize',"-i %s --factor %f"%(fnMask,maskObject.getSamplingRate()/TsMaskOut),numberOfMpi=1)
maskXdim, _, _, _ =img.getDimensions((1,fnMask))
if XdimOut!=maskXdim:
self.runJob('xmipp_transform_window',"-i %s --size %d"%(fnMask,XdimOut),numberOfMpi=1)
[docs] def calculateAngStep(self,newXdim,TsCurrent,ResolutionAlignment):
k=newXdim*TsCurrent/ResolutionAlignment # Freq. index
return math.atan2(1,k)*180.0/math.pi # Corresponding angular step
[docs] def globalAssignment(self,iteration):
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
makePath(fnDirCurrent)
previousResolution=self.readInfoField(fnDirPrevious,"resolution",emlib.MDL_RESOLUTION_FREQREAL)
if self.alignmentMethod==self.GLOBAL_ALIGNMENT or self.alignmentMethod==self.AUTOMATIC_ALIGNMENT or \
self.alignmentMethod==self.STOCHASTIC_ALIGNMENT:
fnGlobal=join(fnDirCurrent,"globalAssignment")
makePath(fnGlobal)
targetResolution=max(previousResolution*0.5,self._maximumTargetResolution[iteration-1])
if self.multiresolution:
TsCurrent=max(self.TsOrig,targetResolution/3)
else:
TsCurrent=self.TsOrig
getShiftsFrom=''
# if iteration>1: # This causes images to be replicated
# getShiftsFrom=fnDirPrevious
self.prepareImages(fnDirPrevious,fnGlobal,TsCurrent,getShiftsFrom)
TsCurrent=self.readInfoField(fnGlobal,"sampling",emlib.MDL_SAMPLINGRATE) # Prepare images may have changed it
self.prepareReferences(fnDirPrevious,fnGlobal,TsCurrent,targetResolution)
# Calculate angular step at this resolution
ResolutionAlignment=previousResolution
if self.nextLowPass:
ResolutionAlignment+=self.nextResolutionOffset.get()
newXdim=self.readInfoField(fnGlobal,"size",emlib.MDL_XSIZE)
angleStep=self.calculateAngStep(newXdim, TsCurrent, ResolutionAlignment)
angleStep=max(angleStep,3.0)
self.writeInfoField(fnGlobal,"angleStep",emlib.MDL_ANGLE_DIFF,float(angleStep))
# Global alignment
perturbationList = [chr(x) for x in range(ord('a'),ord('a')+self.numberOfPerturbations.get())]
for i in range(1,3):
fnDirSignificant=join(fnGlobal,"significant%02d"%i)
fnImgs=join(fnGlobal,"images%02d.xmd"%i)
makePath(fnDirSignificant)
# Create defocus groups
row=getFirstRow(fnImgs)
if row.containsLabel(emlib.MDL_CTF_MODEL) or row.containsLabel(emlib.MDL_CTF_DEFOCUSU):
self.runJob("xmipp_ctf_group","--ctfdat %s -o %s/ctf:stk --pad 1.0 --sampling_rate %f --phase_flipped --error 0.1 --resol %f"%\
(fnImgs,fnDirSignificant,TsCurrent,targetResolution),numberOfMpi=1)
moveFile("%s/ctf_images.sel"%fnDirSignificant,"%s/ctf_groups.xmd"%fnDirSignificant)
cleanPath("%s/ctf_split.doc"%fnDirSignificant)
mdInfo = emlib.MetaData("numberGroups@%s"%join(fnDirSignificant,"ctfInfo.xmd"))
fnCTFs="%s/ctf_ctf.stk"%fnDirSignificant
numberGroups=mdInfo.getValue(emlib.MDL_COUNT,mdInfo.firstObject())
ctfPresent=True
else:
numberGroups=1
ctfPresent=False
fnCTFs=""
# Generate projections
fnReferenceVol=join(fnGlobal,"volumeRef%02d.vol"%i)
for subset in perturbationList:
fnGallery = join(fnDirSignificant,"gallery%02d%s.stk" % (i, subset))
fnGalleryMd = join(fnDirSignificant,"gallery%02d%s.xmd" % (i, subset))
args = "-i %s -o %s --sampling_rate %f --sym %s --min_tilt_angle %f --max_tilt_angle %f --perturb %f " % \
(fnReferenceVol, fnGallery, angleStep,self.symmetryGroup, self.angularMinTilt.get(),self.angularMaxTilt.get(),math.sin(angleStep * math.pi / 180.0) / 4)
args += " --compute_neighbors --angular_distance -1 --experimental_images %s" % self._getExtraPath("images.xmd")
self.runJob("xmipp_angular_project_library", args,numberOfMpi=min(self.numberOfMpi.get(), 24))
cleanPath(join(fnDirSignificant, "gallery_angles%02d%s.doc" % (i, subset)))
moveFile(join(fnDirSignificant,"gallery%02d%s.doc" % (i, subset)), fnGalleryMd)
fnAngles = join(fnGlobal, "anglesDisc%02d%s.xmd" % (i, subset))
for j in range(1, numberGroups + 1):
fnAnglesGroup = join(fnDirSignificant,"angles_group%03d%s.xmd" % (j, subset))
if not exists(fnAnglesGroup):
if ctfPresent:
fnGroup="ctfGroup%06d@%s/ctf_groups.xmd"%(j,fnDirSignificant)
fnCTF ="%d@%s/ctf_ctf.stk"%(j,fnDirSignificant)
fnGalleryGroup=fnGallery=join(fnDirSignificant,"gallery%02d%s_%06d.stk"%(i,subset,j))
fnGalleryGroupMd=fnGallery=join(fnDirSignificant,"gallery%02d%s_%06d.xmd"%(i,subset,j))
self.runJob("xmipp_transform_filter","-i %s -o %s --fourier binary_file %s --save_metadata_stack %s --keep_input_columns"%\
(fnGalleryMd,fnGalleryGroup,fnCTF,fnGalleryGroupMd),
numberOfMpi=min(self.numberOfMpi.get(),24))
else:
fnGroup=fnImgs
fnGalleryGroup=fnGallery
fnGalleryGroupMd=fnGalleryMd
if getSize(fnGroup)==0: # If the group is empty
continue
maxShift = round(self.angularMaxShift.get() * newXdim / 100)
R = self.particleRadius.get()
if R <= 0:
R = self.inputParticles.get().getDimensions()[0] / 2
R = R * self.TsOrig / TsCurrent
if not self.useGpu.get():
args = '-i %s --initgallery %s --maxShift %d --odir %s --dontReconstruct --useForValidation %d --dontCheckMirrors ' % \
(fnGroup, fnGalleryGroupMd, maxShift,fnDirSignificant,self.numberOfReplicates.get() - 1)
self.runJob('xmipp_reconstruct_significant',args,numberOfMpi=self.numberOfMpi.get())
# moveFile(join(fnDirSignificant,"images_significant_iter001_00.xmd"),join(fnDirSignificant,"angles_group%03d%s.xmd"%(j,subset)))
fnAnglesSignificant = join(fnDirSignificant,"angles_iter001_00.xmd")
if exists(fnAnglesSignificant):
moveFile(fnAnglesSignificant, fnAnglesGroup)
cleanPath(join(fnDirSignificant,"images_iter001_00.xmd"))
# cleanPath(join(fnDirSignificant,"angles_iter001_00.xmd"))
cleanPath(join(fnDirSignificant,"images_significant_iter001_00.xmd"))
else:
gpuID = self.setGPU(oneGPU=False)
args = '-i %s -r %s -o %s --keepBestN %f --dev %s ' % \
(fnGroup, fnGalleryGroupMd, fnAnglesGroup, self.numberOfReplicates.get(), gpuID)
self.runJob(CUDA_ALIGN_SIGNIFICANT,args, numberOfMpi=1)
if exists(fnAnglesGroup):
if not exists(fnAngles) and exists(fnAnglesGroup):
copyFile(fnAnglesGroup, fnAngles)
else:
if exists(fnAngles) and exists(fnAnglesGroup):
self.runJob("xmipp_metadata_utilities","-i %s --set union_all %s"%(fnAngles,fnAnglesGroup),numberOfMpi=1)
if exists(fnAngles) and exists(fnImgs):
self.runJob("xmipp_metadata_utilities","-i %s --set join %s image"%(fnAngles,fnImgs),numberOfMpi=1)
if self.saveSpace and ctfPresent:
self.runJob("rm -f",fnDirSignificant+"/gallery*",numberOfMpi=1)
# Evaluate the stability of the alignment
fnOut=join(fnGlobal,"anglesDisc%02d"%i)
for subset1 in perturbationList:
fnOut1=join(fnGlobal,"anglesDisc%02d%s"%(i,subset1))
fnAngles1=fnOut1+".xmd"
counter2 = 0
for subset2 in perturbationList:
if subset1==subset2:
continue
fnAngles2=join(fnGlobal,"anglesDisc%02d%s.xmd"%(i,subset2))
fnOut12=join(fnGlobal,"anglesDisc%02d%s%s"%(i,subset1,subset2))
self.runJob("xmipp_angular_distance","--ang1 %s --ang2 %s --oroot %s --sym %s --compute_weights 1 particleId 0.5 --check_mirrors --set 0"%(fnAngles2,fnAngles1,fnOut12,self.symmetryGroup),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "angleDiff0 shiftDiff0 weightJumper0"'%(fnOut12+"_weights.xmd"),numberOfMpi=1)
if counter2 == 0:
mdWeightsAll = emlib.MetaData(fnOut12+"_weights.xmd")
counter2=1
else:
mdWeights = emlib.MetaData(fnOut12+"_weights.xmd")
if mdWeights.size()==mdWeightsAll.size():
counter2 += 1
for id1, id2 in izip(mdWeights,mdWeightsAll):
angleDiff0 = mdWeights.getValue(emlib.MDL_ANGLE_DIFF0, id1)
shiftDiff0 = mdWeights.getValue(emlib.MDL_SHIFT_DIFF0, id1)
weightJumper0 = mdWeights.getValue(emlib.MDL_WEIGHT_JUMPER0, id1)
angleDiff0All = mdWeightsAll.getValue(emlib.MDL_ANGLE_DIFF0, id2)
shiftDiff0All = mdWeightsAll.getValue(emlib.MDL_SHIFT_DIFF0, id2)
weightJumper0All = mdWeightsAll.getValue(emlib.MDL_WEIGHT_JUMPER0, id2)
mdWeightsAll.setValue(emlib.MDL_ANGLE_DIFF0, angleDiff0+angleDiff0All, id2)
mdWeightsAll.setValue(emlib.MDL_SHIFT_DIFF0, shiftDiff0+shiftDiff0All, id2)
mdWeightsAll.setValue(emlib.MDL_WEIGHT_JUMPER0, weightJumper0+weightJumper0All, id2)
if counter2>1:
iCounter2 = 1.0/counter2
for id in mdWeightsAll:
angleDiff0All = mdWeightsAll.getValue(emlib.MDL_ANGLE_DIFF0, id)
shiftDiff0All = mdWeightsAll.getValue(emlib.MDL_SHIFT_DIFF0, id)
weightJumper0All = mdWeightsAll.getValue(emlib.MDL_WEIGHT_JUMPER0, id)
mdWeightsAll.setValue(emlib.MDL_ANGLE_DIFF0, angleDiff0All*iCounter2, id)
mdWeightsAll.setValue(emlib.MDL_SHIFT_DIFF0, shiftDiff0All*iCounter2, id)
mdWeightsAll.setValue(emlib.MDL_WEIGHT_JUMPER0, weightJumper0All*iCounter2, id)
if counter2>0:
mdWeightsAll.write(fnOut1+"_weights.xmd")
self.runJob("xmipp_metadata_utilities",'-i %s --set merge %s'%(fnAngles1,fnOut1+"_weights.xmd"),numberOfMpi=1)
if not exists(fnOut+".xmd") and exists(fnAngles1):
copyFile(fnAngles1,fnOut+".xmd")
else:
if exists(fnAngles1) and exists(fnOut+".xmd"):
self.runJob("xmipp_metadata_utilities",'-i %s --set union_all %s'%(fnOut+".xmd",fnAngles1),numberOfMpi=1)
cleanPath(join(fnGlobal,"anglesDisc*_weights.xmd"))
[docs] def adaptShifts(self, fnSource, TsSource, fnDest, TsDest):
K=TsSource/TsDest
copyFile(fnSource,fnDest)
row=getFirstRow(fnDest)
if row.containsLabel(emlib.MDL_SHIFT_X):
self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "shiftX=%f*shiftX"'%(fnDest,K),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "shiftY=%f*shiftY"'%(fnDest,K),numberOfMpi=1)
if row.containsLabel(emlib.MDL_CONTINUOUS_X):
self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "continuousX=%f*continuousX"'%(fnDest,K),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities','-i %s --operate modify_values "continuousY=%f*continuousY"'%(fnDest,K),numberOfMpi=1)
[docs] def localAssignment(self,iteration):
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
if self.alignmentMethod==self.LOCAL_ALIGNMENT or \
(self.alignmentMethod==self.AUTOMATIC_ALIGNMENT and iteration>=4):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
fnDirLocal=join(fnDirCurrent,"localAssignment")
makePath(fnDirLocal)
previousResolution=self.readInfoField(fnDirPrevious,"resolution",emlib.MDL_RESOLUTION_FREQREAL)
targetResolution=max(previousResolution*0.8,self._maximumTargetResolution[iteration-1])
if self.multiresolution:
TsCurrent=max(self.TsOrig,targetResolution/3)
else:
TsCurrent=self.TsOrig
self.writeInfoField(fnDirLocal,"sampling",emlib.MDL_SAMPLINGRATE,TsCurrent)
TsCurrent=self.readInfoField(fnDirLocal,"sampling",emlib.MDL_SAMPLINGRATE) # Write and read to guarantee consistency with previous directories
# Prepare images and references
produceNewReferences=True
fnDirGlobal=join(fnDirCurrent,"globalAssignment")
if exists(fnDirGlobal):
TsGlobal=self.readInfoField(fnDirGlobal,"sampling",emlib.MDL_SAMPLINGRATE)
if TsGlobal==TsCurrent:
produceNewReferences=False
if produceNewReferences:
self.prepareImages(fnDirPrevious,fnDirLocal,TsCurrent,fnDirPrevious)
TsCurrent=self.readInfoField(fnDirLocal,"sampling",emlib.MDL_SAMPLINGRATE) # Prepare images may have changed it
self.prepareReferences(fnDirPrevious,fnDirLocal,TsCurrent,targetResolution)
else:
newXdim=self.readInfoField(fnDirGlobal,"size",emlib.MDL_XSIZE)
self.writeInfoField(fnDirLocal,"size",emlib.MDL_XSIZE,newXdim)
for i in range(1,3):
createLink(join(fnDirGlobal,"images%02d.xmd"%i),join(fnDirLocal,"images%02d.xmd"%i))
createLink(join(fnDirGlobal,"volumeRef%02d.vol"%i),join(fnDirLocal,"volumeRef%02d.vol"%i))
# Compute maximum angular deviation
ResolutionAlignment=previousResolution
if self.nextLowPass:
ResolutionAlignment+=self.nextResolutionOffset.get()
newXdim=self.readInfoField(fnDirLocal,"size",emlib.MDL_XSIZE)
maxAngle=3*self.calculateAngStep(newXdim, TsCurrent, ResolutionAlignment)
for i in range(1,3):
state = self.readInfoField(fnDirLocal, "count", emlib.MDL_COUNT)
if state>=2+i:
continue
fnLocalImages=join(fnDirLocal,"images%02d.xmd"%i)
fnLocalImagesIdx=join(fnDirLocal,"images%02d_idx.xmd"%i)
# Starting angles
fnLocalAssignment=join(fnDirLocal,"anglesDisc%02d.xmd"%i)
if exists(fnDirGlobal):
fnGlobalAssignment=join(fnDirGlobal,"anglesDisc%02d.xmd"%i)
TsGlobal=self.readInfoField(fnDirGlobal,"sampling",emlib.MDL_SAMPLINGRATE)
if TsGlobal==TsCurrent:
copyFile(fnGlobalAssignment,fnLocalAssignment)
else:
self.adaptShifts(fnGlobalAssignment,TsGlobal,fnLocalAssignment,TsCurrent)
else:
TsPrevious=self.readInfoField(fnDirPrevious,"sampling",emlib.MDL_SAMPLINGRATE)
fnAux=join(fnDirLocal,"aux.xmd")
self.runJob("xmipp_metadata_utilities","-i %s --set intersection %s particleId particleId -o %s"%\
(join(fnDirPrevious,"angles.xmd"),fnLocalImages,fnAux),numberOfMpi=1)
self.adaptShifts(fnAux,TsPrevious,fnLocalAssignment,TsCurrent)
cleanPath(fnAux)
self.runJob("xmipp_metadata_utilities","-i %s --operate drop_column image"%fnLocalAssignment,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "particleId image" -o %s'%(fnLocalImages,fnLocalImagesIdx),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate remove_duplicates particleId'%(fnLocalImagesIdx),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities","-i %s --set join %s particleId"%(fnLocalAssignment,fnLocalImagesIdx),numberOfMpi=1)
cleanPath(fnLocalImagesIdx)
fnVol=join(fnDirLocal,"volumeRef%02d.vol"%i)
fnLocalStk=join(fnDirLocal,"anglesCont%02d.stk"%i)
R=self.particleRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2
R=round(R*self.TsOrig/TsCurrent)
args="-i %s -o %s --sampling %f --Rmax %d --padding %d --ref %s --max_resolution %f --applyTo image1"%\
(fnLocalAssignment,fnLocalStk,TsCurrent,R,self.contPadding.get(),fnVol,previousResolution)
if self.contShift or self.alignmentMethod.get()==self.AUTOMATIC_ALIGNMENT:
args+=" --optimizeShift --max_shift %f"%(self.contMaxShiftVariation.get()*newXdim*0.01)
if self.contScale or (self.alignmentMethod.get()==self.AUTOMATIC_ALIGNMENT and iteration>=5):
args+=" --optimizeScale --max_scale %f"%self.contMaxScale.get()
if self.contAngles or self.alignmentMethod.get()==self.AUTOMATIC_ALIGNMENT:
args+=" --optimizeAngles --max_angular_change %f"%maxAngle
if self.contDefocus or (self.alignmentMethod.get()==self.AUTOMATIC_ALIGNMENT and iteration>=5):
args+=" --optimizeDefocus --max_defocus_change %f"%self.contMaxDefocus.get()
if self.inputParticles.get().isPhaseFlipped():
args+=" --phaseFlipped"
#if self.weightResiduals:
# args+=" --oresiduals %s"%join(fnDirLocal,"residuals%02i.stk"%i)
#
if self.useGpu:
args+=" --nThreads %d"%self.numberOfMpi.get()
self.runJob('xmipp_cuda_angular_continuous_assign2', args, numberOfMpi=1)
else:
self.runJob("xmipp_angular_continuous_assign2", args, numberOfMpi=self.numberOfMpi.get())
self.runJob("xmipp_transform_mask","-i %s --mask circular -%d"%(fnLocalStk,R),numberOfMpi=min(self.numberOfMpi.get(),24))
self.writeInfoField(fnDirLocal,"count",emlib.MDL_COUNT,int(2+i))
[docs] def noAssignment(self, iteration):
fnDirPrevious = self._getExtraPath("Iter%03d" % (iteration - 1))
fnDirCurrent = self._getExtraPath("Iter%03d" % iteration)
fnDirLocal = join(fnDirCurrent, "noAssignment")
makePath(fnDirLocal)
previousResolution = self.readInfoField(fnDirPrevious, "resolution", emlib.MDL_RESOLUTION_FREQREAL)
targetResolution = max(previousResolution * 0.8, self._maximumTargetResolution[iteration - 1])
if self.multiresolution:
TsCurrent = max(self.TsOrig, targetResolution / 3)
else:
TsCurrent = self.TsOrig
self.writeInfoField(fnDirLocal, "sampling", emlib.MDL_SAMPLINGRATE, TsCurrent)
TsCurrent = self.readInfoField(fnDirLocal, "sampling",
emlib.MDL_SAMPLINGRATE) # Write and read to guarantee consistency with previous directories
# Prepare images and references
self.prepareImages(fnDirPrevious, fnDirLocal, TsCurrent)
for i in range(1, 3):
fnLocalImages = join(fnDirLocal, "images%02d.xmd" % i)
fnAssignment = join(fnDirLocal, "anglesNoAssignment%02d.xmd" % i)
TsPrevious = self.readInfoField(fnDirPrevious, "sampling", emlib.MDL_SAMPLINGRATE)
self.adaptShifts(fnLocalImages, TsPrevious, fnAssignment, TsCurrent)
row = md.getFirstRow(fnAssignment)
if not row.hasLabel(md.MDL_MAXCC) and not row.hasLabel(md.MDL_COST):
self.runJob("xmipp_metadata_utilities", "-i %s --fill maxCC constant 1"% fnAssignment, numberOfMpi=1)
[docs] def weightParticles(self, iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
from math import exp
for i in range(1,3):
# Grab file
fnDirGlobal=join(fnDirCurrent,"globalAssignment")
fnDirLocal=join(fnDirCurrent,"localAssignment")
fnDirNo=join(fnDirCurrent,"noAssignment")
fnAnglesCont=join(fnDirLocal,"anglesCont%02d.xmd"%i)
fnAnglesDisc=join(fnDirGlobal,"anglesDisc%02d.xmd"%i)
fnAnglesNo=join(fnDirNo,"anglesNoAssignment%02d.xmd"%i)
fnAngles=join(fnDirCurrent,"angles%02d.xmd"%i)
if exists(fnAnglesCont):
copyFile(fnAnglesCont, fnAngles)
TsCurrent=self.readInfoField(fnDirLocal,"sampling",emlib.MDL_SAMPLINGRATE)
Xdim=self.readInfoField(fnDirLocal,"size",emlib.MDL_XSIZE)
elif exists(fnAnglesNo):
copyFile(fnAnglesNo, fnAngles)
TsCurrent = self.readInfoField(fnDirNo, "sampling", emlib.MDL_SAMPLINGRATE)
Xdim = self.readInfoField(fnDirNo, "size", emlib.MDL_XSIZE)
else:
if exists(fnAnglesDisc):
copyFile(fnAnglesDisc, fnAngles)
TsCurrent=self.readInfoField(fnDirGlobal,"sampling",emlib.MDL_SAMPLINGRATE)
Xdim=self.readInfoField(fnDirGlobal,"size",emlib.MDL_XSIZE)
else:
raise Exception("Angles for iteration "+str(iteration)+" not found")
self.writeInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE,TsCurrent)
self.writeInfoField(fnDirCurrent,"size",emlib.MDL_XSIZE,Xdim)
if self.weightSSNR:
row=getFirstRow(fnAngles)
if row.containsLabel(emlib.MDL_WEIGHT_SSNR):
self.runJob("xmipp_metadata_utilities","-i %s --operate drop_column weightSSNR"%fnAngles,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities","-i %s --set join %s particleId"%\
(fnAngles,self._getExtraPath("ssnrWeights.xmd")),numberOfMpi=1)
if iteration>1 and self.alignmentMethod!=self.STOCHASTIC_ALIGNMENT:
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
if self.splitMethod == self.SPLIT_FIXED:
fnPreviousAngles=join(fnDirPrevious,"angles%02d.xmd"%i)
else:
fnPreviousAngles=join(fnDirCurrent,"aux.xmd")
self.runJob("xmipp_metadata_utilities","-i %s --set intersection %s particleId particleId -o %s"%\
(join(fnDirPrevious,"angles.xmd"),fnAngles,fnPreviousAngles),numberOfMpi=1)
self.runJob("xmipp_angular_distance","--ang1 %s --ang2 %s --compute_weights --check_mirrors --oroot %s --sym %s"%\
(fnPreviousAngles,fnAngles,fnDirCurrent+"/jumper",self.symmetryGroup),numberOfMpi=1)
moveFile(fnDirCurrent+"/jumper_weights.xmd", fnAngles)
if self.splitMethod == self.SPLIT_STOCHASTIC:
cleanPath(fnPreviousAngles)
if iteration>2:
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-2))
if self.splitMethod == self.SPLIT_FIXED:
fnPreviousAngles=join(fnDirPrevious,"angles%02d.xmd"%i)
else:
fnPreviousAngles=join(fnDirCurrent,"aux.xmd")
self.runJob("xmipp_metadata_utilities","-i %s --set intersection %s particleId particleId -o %s"%\
(join(fnDirPrevious,"angles.xmd"),fnAngles,fnPreviousAngles),numberOfMpi=1)
self.runJob("xmipp_angular_distance","--ang1 %s --ang2 %s --compute_weights --check_mirrors --oroot %s --set 2 --sym %s"%\
(fnPreviousAngles,fnAngles,fnDirCurrent+"/jumper",self.symmetryGroup),numberOfMpi=1)
moveFile(fnDirCurrent+"/jumper_weights.xmd", fnAngles)
if self.splitMethod == self.SPLIT_STOCHASTIC:
cleanPath(fnPreviousAngles)
#if self.weightResiduals and exists(fnAnglesCont):
# fnCovariance=join(fnDirLocal,"covariance%02d.stk"%i)
# self.runJob("xmipp_image_residuals","-i %s -o %s --normalizeDivergence"%(fnAngles,fnCovariance),numberOfMpi=1)
# moveFile(join(fnDirLocal,"covariance%02d.xmd"%i),fnAngles)
mdAngles=emlib.MetaData(fnAngles)
weightCCmin=float(self.weightCCmin.get())
for objId in mdAngles:
weight=1.0
if self.weightJumper and self.alignmentMethod==self.GLOBAL_ALIGNMENT and self.numberOfPerturbations.get()>1:
aux=mdAngles.getValue(emlib.MDL_WEIGHT_JUMPER0,objId)
weight*=aux
if self.weightSSNR:
aux=mdAngles.getValue(emlib.MDL_WEIGHT_SSNR,objId)
weight*=aux
if self.weightContinuous and exists(fnAnglesCont) and self.alignmentMethod==self.LOCAL_ALIGNMENT:
aux=mdAngles.getValue(emlib.MDL_WEIGHT_CONTINUOUS2,objId)
weight*=aux
#if self.weightResiduals and exists(fnAnglesCont):
# aux=mdAngles.getValue(emlib.MDL_ZSCORE_RESCOV,objId)
# aux/=3
# weight*=exp(-0.5*aux*aux)
# aux=mdAngles.getValue(emlib.MDL_ZSCORE_RESMEAN,objId)
# aux/=3
# weight*=exp(-0.5*aux*aux)
# aux=mdAngles.getValue(emlib.MDL_ZSCORE_RESVAR,objId)
# aux/=3
# weight*=exp(-0.5*aux*aux)
if self.weightJumper and iteration>1:
w1=mdAngles.getValue(emlib.MDL_WEIGHT_JUMPER,objId)
w2=1.0
if iteration>2:
w2=mdAngles.getValue(emlib.MDL_WEIGHT_JUMPER2,objId)
weight*=w1*w2
mdAngles.setValue(emlib.MDL_WEIGHT,weight,objId)
mdAngles.write(fnAngles)
fnAngles=join(fnDirCurrent,"angles.xmd")
fnAngles1=join(fnDirCurrent,"angles01.xmd")
fnAngles2=join(fnDirCurrent,"angles02.xmd")
self.runJob('xmipp_metadata_utilities',"-i %s --set union %s -o %s"%(fnAngles1,fnAngles2,fnAngles),numberOfMpi=1)
[docs] def qualifyParticles(self, iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
fnAngles=join(fnDirCurrent,"angles.xmd")
fnAnglesQualified=join(fnDirCurrent,"angles_qualified.xmd")
# Qualify according to CC and COST by defocus groups
row=getFirstRow(fnAngles)
if row.containsLabel(emlib.MDL_CTF_MODEL) or row.containsLabel(emlib.MDL_CTF_DEFOCUSU):
previousResolution=self.readInfoField(fnDirPrevious,"resolution",emlib.MDL_RESOLUTION_FREQREAL)
TsCurrent=self.readInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE)
numberGroups=50
self.runJob("xmipp_ctf_group","--ctfdat %s -o %s/ctf:stk --simple %d"%\
(fnAngles,fnDirCurrent,numberGroups),numberOfMpi=1)
moveFile("%s/ctf_images.sel"%fnDirCurrent,"%s/ctf_groups.xmd"%fnDirCurrent)
ctfPresent=True
else:
numberGroups=1
ctfPresent=False
for j in range(1,numberGroups+2):
fnAnglesGroup=join(fnDirCurrent,"angles_group%03d.xmd"%j)
if ctfPresent:
fnGroup="ctfGroup%06d@%s/ctf_groups.xmd"%(j,fnDirCurrent)
else:
fnGroup=fnAngles
if getSize(fnGroup)>0:
if row.containsLabel(emlib.MDL_MAXCC):
self.runJob("xmipp_metadata_utilities","-i %s --operate percentile maxCC maxCCPerc -o %s"%(fnGroup,fnAnglesGroup),numberOfMpi=1)
fnGroup=fnAnglesGroup
if row.containsLabel(emlib.MDL_COST):
self.runJob("xmipp_metadata_utilities","-i %s --operate percentile cost costPerc -o %s"%(fnGroup,fnAnglesGroup),numberOfMpi=1)
if not exists(fnAnglesQualified):
copyFile(fnAnglesGroup, fnAnglesQualified)
else:
self.runJob("xmipp_metadata_utilities","-i %s --set union %s"%(fnAnglesQualified,fnAnglesGroup),numberOfMpi=1)
cleanPath(fnAnglesGroup)
if ctfPresent:
cleanPath("%s/ctf_groups.xmd"%fnDirCurrent)
moveFile(fnAnglesQualified, fnAngles)
if self.weightCC:
mdAngles=emlib.MetaData(fnAngles)
weightCCmin=float(self.weightCCmin.get())
hasCost = mdAngles.containsLabel(md.MDL_COST_PERCENTILE)
hasCC = mdAngles.containsLabel(md.MDL_MAXCC_PERCENTILE)
for objId in mdAngles:
if self.alignmentMethod==self.LOCAL_ALIGNMENT or \
(self.alignmentMethod==self.NO_ALIGNMENT and hasCost):
w=mdAngles.getValue(emlib.MDL_COST_PERCENTILE,objId)
elif hasCC:
w=mdAngles.getValue(emlib.MDL_MAXCC_PERCENTILE,objId)
else:
w=1
weight=mdAngles.getValue(emlib.MDL_WEIGHT,objId)
weight*=weightCCmin+w*(1-weightCCmin)
mdAngles.setValue(emlib.MDL_WEIGHT,weight,objId)
mdAngles.write(fnAngles)
# Qualify according to angular temperature
if iteration==1:
self.runJob("xmipp_metadata_utilities","-i %s --fill angleTemp constant 0"%fnAngles,numberOfMpi=1)
else:
fnAngles_1=join(fnDirPrevious,"angles.xmd")
fnTemp=join(fnDirCurrent,"previousAngleTemp.xmd")
self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "particleId angleTemp" -o %s'%(fnAngles_1,fnTemp),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate remove_duplicates particleId'%fnTemp,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --set join %s particleId'%(fnAngles,fnTemp),numberOfMpi=1)
cleanPath(fnTemp)
hasDiff0 = row.containsLabel(emlib.MDL_ANGLE_DIFF0)
hasDiff1 = row.containsLabel(emlib.MDL_ANGLE_DIFF)
hasDiff2 = row.containsLabel(emlib.MDL_ANGLE_DIFF2)
if hasDiff0 or hasDiff1 or hasDiff2:
mdAngles=emlib.MetaData(fnAngles)
K=1.0/3.0
iNorm = 1.0/(hasDiff0 + hasDiff1 + hasDiff2)
for objId in mdAngles:
perturbation=0.
if hasDiff0:
perturbation+=mdAngles.getValue(emlib.MDL_ANGLE_DIFF0,objId)
if hasDiff1:
perturbation+=mdAngles.getValue(emlib.MDL_ANGLE_DIFF,objId)
if hasDiff2:
perturbation+=mdAngles.getValue(emlib.MDL_ANGLE_DIFF2,objId)
perturbation*=iNorm
previousTemp = mdAngles.getValue(emlib.MDL_ANGLE_TEMPERATURE,objId)
currentTemp = (1-K)*previousTemp + K*perturbation
mdAngles.setValue(emlib.MDL_ANGLE_TEMPERATURE,currentTemp,objId)
mdAngles.write(fnAngles)
[docs] def reconstruct(self, iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
TsCurrent=self.readInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE)
# Delete previous image files, they exist in case that the last iteration
# was performed as a single iteration
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
fnCorrectedImages1=join(fnDirPrevious,"images_corrected01.stk")
if exists(fnCorrectedImages1) and self.saveSpace.get():
cleanPath(fnCorrectedImages1)
fnCorrectedImages2=join(fnDirPrevious,"images_corrected02.stk")
if exists(fnCorrectedImages2) and self.saveSpace.get():
cleanPath(fnCorrectedImages2)
grayAdjusted=False
for i in range(1,3):
fnAngles=join(fnDirCurrent,"angles%02d.xmd"%i)
fnVol=join(fnDirCurrent,"volume%02d.vol"%i)
if not exists(fnVol):
# Correct for the CTF
fnAnglesToUse = fnAngles
row=getFirstRow(fnAngles)
hasCTF = row.containsLabel(emlib.MDL_CTF_DEFOCUSU) or row.containsLabel(emlib.MDL_CTF_MODEL)
fnCorrectedImagesRoot=join(fnDirCurrent,"images_corrected%02d"%i)
deleteStack = False
if hasCTF:
args="-i %s -o %s.stk --save_metadata_stack %s.xmd --keep_input_columns"%(fnAngles,fnCorrectedImagesRoot,fnCorrectedImagesRoot)
args+=" --sampling_rate %f"%TsCurrent
if self.correctEnvelope:
args+=" --correct_envelope"
if self.inputParticles.get().isPhaseFlipped():
args+=" --phase_flipped"
self.runJob("xmipp_ctf_correct_wiener2d",args,numberOfMpi=min(self.numberOfMpi.get(),24))
self.runJob("xmipp_image_eliminate_byEnergy","-i %s.xmd --sigma2 9 --minSigma2 0.01"%\
fnCorrectedImagesRoot,numberOfMpi=min(self.numberOfMpi.get(),12))
fnAnglesToUse = fnCorrectedImagesRoot+".xmd"
deleteStack = True
deletePattern = fnCorrectedImagesRoot+".*"
if self.alignmentMethod!=self.STOCHASTIC_ALIGNMENT:
self.runJob('xmipp_metadata_utilities','-i %s --set intersection %s particleId particleId'%(fnAngles,fnAnglesToUse),numberOfMpi=1)
# This is because eliminate_byEnergy may have reduced the number of images in fnAngles
if (self.contGrayValues and self.alignmentMethod.get()==self.LOCAL_ALIGNMENT) or \
(self.alignmentMethod.get()==self.AUTOMATIC_ALIGNMENT and iteration>=5):
grayAdjusted=True
R=self.particleRadius.get()
if R<=0:
R=self.inputParticles.get().getDimensions()[0]/2
fnGrayRoot = join(fnDirCurrent,"images_gray%02d"%i)
fnRefVol=join(fnDirCurrent,"localAssignment","volumeRef%02d.vol"%i)
fnDirPrevious=self._getExtraPath("Iter%03d"%(iteration-1))
previousResolution=self.readInfoField(fnDirPrevious,"resolution",emlib.MDL_RESOLUTION_FREQREAL)
args="-i %s -o %s.stk --sampling %f --Rmax %d --padding %d --ref %s --max_resolution %f --save_metadata_stack %s.xmd"%\
(fnAnglesToUse,fnGrayRoot,TsCurrent,R,self.contPadding.get(),fnRefVol,previousResolution,fnGrayRoot)
args+=" --max_gray_scale %f --max_gray_shift %f"%\
(self.contMaxGrayScale.get(),self.contMaxGrayShift.get())
self.runJob("xmipp_transform_adjust_image_grey_levels",args,numberOfMpi=self.numberOfMpi.get())
fnAnglesToUse = fnGrayRoot+".xmd"
if deleteStack:
cleanPattern(deletePattern)
deleteStack = True
deletePattern = fnGrayRoot+".*"
# Save the gray transformation
self.runJob("xmipp_metadata_utilities",'-i %s --operate drop_column "continuousA continuousB"'%fnAngles,numberOfMpi=1)
fnAux = join(fnDirCurrent,"gray_transformation%02d.xmd"%i)
self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "continuousA continuousB" -o %s'%(fnAnglesToUse,fnAux),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --set merge %s'%(fnAngles,fnAux),numberOfMpi=1)
cleanPath(fnAux)
# Restrict the angles
if self.restrictReconstructionAngles:
fnRestricted = join(fnDirCurrent,"images_restricted%02d.xmd"%i)
args = '-i %s --query select "angleRot > %f AND anglePsi < %f" -o %s'%\
(fnAnglesToUse,self.angularMinTiltReconstruct.get(),self.angularMaxTiltReconstruct.get(),fnRestricted)
self.runJob("xmipp_metadata_utilities",args,numberOfMpi=1)
fnAnglesToUse = fnRestricted
# Reconstruct Fourier
args="-i %s -o %s --sym %s --weight"%(fnAnglesToUse,fnVol,self.symmetryGroup)
if self.useGpu.get():
#AJ to make it work with and without queue system
if self.numberOfMpi.get()>1:
gpuId = self.setGPU(oneGPU=False)
args += ' -gpusPerNode %d' % self.getNumGpus()
args += ' -threadsPerGPU %d' % max(self.numberOfThreads.get(),4)
if self.numberOfMpi.get()==1:
gpuId = self.setGPU(oneGPU=True)
args += " --device %s" % (gpuId.replace(',', ' '))
args += ' --thr %s' % self.numberOfThreads.get()
if self.numberOfMpi.get()>1:
self.runJob('xmipp_cuda_reconstruct_fourier', args, numberOfMpi=self.getNumGpus()+1)
else:
self.runJob('xmipp_cuda_reconstruct_fourier', args)
else:
self.runJob("xmipp_reconstruct_fourier_accel", args, numberOfMpi=self.numberOfMpi.get())
# If stochastic gradient descent
if self.alignmentMethod==self.STOCHASTIC_ALIGNMENT:
newXdim = self.readInfoField(fnDirCurrent, "size", emlib.MDL_XSIZE)
fnAuxVol=join(fnDirCurrent,"volume%02d_aux.vol"%i)
fnPreviousVol=join(fnDirPrevious,"volume%02d.vol"%i)
self.runJob("xmipp_image_resize","-i %s -o %s --dim %d"%(fnPreviousVol,fnAuxVol,newXdim))
self.runJob("xmipp_image_operate","-i %s --mult %f"%(fnVol,self.alphaSGD.get()))
self.runJob("xmipp_image_operate","-i %s --mult %f"%(fnAuxVol,1-self.alphaSGD.get()))
self.runJob("xmipp_image_operate","-i %s --plus %s"%(fnVol,fnAuxVol))
cleanPath(fnAuxVol)
if deleteStack:
cleanPattern(deletePattern)
if grayAdjusted:
fnAngles=join(fnDirCurrent,"angles.xmd")
fnAnglesAux=join(fnDirCurrent,"anglesAux.xmd")
fnAnglesAuxId=join(fnDirCurrent,"anglesAuxId.xmd")
fnAngles1=join(fnDirCurrent,"angles01.xmd")
fnAngles2=join(fnDirCurrent,"angles02.xmd")
self.runJob("xmipp_metadata_utilities",'-i %s --operate drop_column "continuousA continuousB"'%fnAngles,numberOfMpi=1)
self.runJob('xmipp_metadata_utilities',"-i %s --set union %s -o %s"%(fnAngles1,fnAngles2,fnAnglesAux),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities',"-i %s --operate keep_column itemId -o %s"%\
(fnAnglesAux,fnAnglesAuxId),numberOfMpi=1)
self.runJob('xmipp_metadata_utilities',"-i %s --set intersection %s itemId itemId"%\
(fnAngles,fnAnglesAuxId),numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate sort itemId'%fnAngles,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate sort itemId'%fnAnglesAux,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --operate keep_column "continuousA continuousB"'%fnAnglesAux,numberOfMpi=1)
self.runJob("xmipp_metadata_utilities",'-i %s --set merge %s'%(fnAngles,fnAnglesAux),numberOfMpi=1)
cleanPath(fnAnglesAux)
cleanPath(fnAnglesAuxId)
[docs] def postProcessing(self, iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
TsCurrent=self.readInfoField(fnDirCurrent,"sampling",emlib.MDL_SAMPLINGRATE)
for i in range(1,3):
fnVol=join(fnDirCurrent,"volume%02d.vol"%i)
fnBeforeVol=join(fnDirCurrent,"volumeBeforePostProcessing%02d.vol"%i)
volXdim = self.readInfoField(fnDirCurrent, "size", emlib.MDL_XSIZE)
if self.postSymmetryWithinMask:
if self.postMaskSymmetry!="c1":
fnMask=join(fnDirCurrent,"mask.vol")
self.prepareMask(self.postSymmetryWithinMaskMask.get(),fnMask,TsCurrent,volXdim)
self.runJob("xmipp_transform_symmetrize","-i %s --sym %s --mask_in %s --dont_wrap"%\
(fnVol,self.postSymmetryWithinMaskType.get(),fnMask),numberOfMpi=1)
cleanPath(fnMask)
if self.postSymmetryHelical:
z0=float(self.postSymmetryHelicalMinZ.get())
zF=float(self.postSymmetryHelicalMaxZ.get())
zStep=(zF-z0)/10
rot0=float(self.postSymmetryHelicalMinRot.get())
rotF=float(self.postSymmetryHelicalMaxRot.get())
rotStep=(rotF-rot0)/10
fnCoarse=join(fnDirCurrent,"coarseHelical%02d.xmd"%i)
fnFine=join(fnDirCurrent,"fineHelical%02d.xmd"%i)
radius=int(self.postSymmetryHelicalRadius.get()/TsCurrent)
height=int(volXdim)
self.runCoarseSearch(fnVol, self.postSymmetryHelicalDihedral, 0.9, z0, zF, zStep, rot0, rotF, rotStep, 1, fnCoarse, 0, radius, height, TsCurrent)
self.runFineSearch(fnVol, self.postSymmetryHelicalDihedral, fnCoarse, fnFine, 0.9, z0, zF, rot0, rotF, 0, radius, height, TsCurrent)
cleanPath(fnCoarse)
self.runSymmetrize(fnVol, self.postSymmetryHelicalDihedral, fnFine, fnVol, 0.9, 0, radius, height, TsCurrent)
if self.postScript!="":
img = ImageHandler()
volXdim, _, _, _ =img.getDimensions((1,fnVol))
scriptArgs = {'volume': fnVol,
'sampling': TsCurrent,
'dim': volXdim,
'iterDir': fnDirCurrent}
cmd = self.postScript % scriptArgs
self.runJob(cmd, '', numberOfMpi=1)
# Align volumes
fnVol1=join(fnDirCurrent,"volume%02d.vol"%1)
fnVol2=join(fnDirCurrent,"volume%02d.vol"%2)
fnVolAvg=join(fnDirCurrent,"volumeAvg.mrc")
self.runJob('xmipp_image_operate','-i %s --plus %s -o %s'%(fnVol1,fnVol2,fnVolAvg),numberOfMpi=1)
self.runJob('xmipp_image_operate','-i %s --mult 0.5'%fnVolAvg,numberOfMpi=1)
self.runJob('xmipp_volume_align','--i1 %s --i2 %s --local --apply'%(fnVolAvg,fnVol1),numberOfMpi=1)
self.runJob('xmipp_volume_align','--i1 %s --i2 %s --local --apply'%(fnVolAvg,fnVol2),numberOfMpi=1)
# Generate mask if available
if self.postAdHocMask.hasValue():
fnMask=join(fnDirCurrent,"mask.vol")
if not exists(fnMask):
volXdim = self.readInfoField(fnDirCurrent, "size", emlib.MDL_XSIZE)
self.prepareMask(self.postAdHocMask.get(), fnMask, TsCurrent, volXdim)
self.runJob('xmipp_transform_threshold',"-i %s --select below 0.5 --substitute binarize"%fnMask,numberOfMpi=1)
else:
fnMask=""
# Remove untrusted background voxels
if self.postSignificantDenoise:
fnRootRestored=join(fnDirCurrent,"volumeRestored")
args='--i1 %s --i2 %s --oroot %s --denoising 1'%(fnVol1,fnVol2,fnRootRestored)
if fnMask!="":
args+=" --mask binary_file %s"%fnMask
if self.useGpu:
self.runJob('xmipp_cuda_volume_halves_restoration', args, numberOfMpi=1)
else:
self.runJob('xmipp_volume_halves_restoration',args,numberOfMpi=1)
moveFile("%s_restored1.vol"%fnRootRestored,fnVol1)
moveFile("%s_restored2.vol"%fnRootRestored,fnVol2)
# Filter bank denoising
if self.postFilterBank:
fnRootRestored=join(fnDirCurrent,"volumeRestored")
args='--i1 %s --i2 %s --oroot %s --filterBank 0.01'%(fnVol1,fnVol2,fnRootRestored)
if fnMask!="":
args+=" --mask binary_file %s"%fnMask
if self.useGpu:
self.runJob('xmipp_cuda_volume_halves_restoration', args, numberOfMpi=1)
else:
self.runJob('xmipp_volume_halves_restoration',args,numberOfMpi=1)
moveFile("%s_restored1.vol"%fnRootRestored,fnVol1)
moveFile("%s_restored2.vol"%fnRootRestored,fnVol2)
cleanPath("%s_filterBank.vol"%fnRootRestored)
# Laplacian Denoising
if self.postLaplacian:
fnRootRestored=join(fnDirCurrent,"volumeRestored")
args = "-i %s --retinex 0.95 "
if fnMask!="":
args+=fnMask
self.runJob('xmipp_transform_filter',args%fnVol1,numberOfMpi=1)
self.runJob('xmipp_transform_filter',args%fnVol2,numberOfMpi=1)
# Blind deconvolution
if self.postDeconvolve:
fnRootRestored=join(fnDirCurrent,"volumeRestored")
args='--i1 %s --i2 %s --oroot %s --deconvolution 1'%(fnVol1,fnVol2,fnRootRestored)
if fnMask!="":
args+=" --mask binary_file %s"%fnMask
if self.useGpu:
self.runJob('xmipp_cuda_volume_halves_restoration', args, numberOfMpi=1)
else:
self.runJob('xmipp_volume_halves_restoration',args,numberOfMpi=1)
moveFile("%s_restored1.vol"%fnRootRestored,fnVol1)
moveFile("%s_restored2.vol"%fnRootRestored,fnVol2)
self.runJob("xmipp_image_convert","-i %s_convolved.vol -o %s -t vol"%(fnRootRestored,fnVolAvg),numberOfMpi=1)
cleanPath("%s_convolved.vol"%fnRootRestored)
cleanPath("%s_deconvolved.vol"%fnRootRestored)
fnForFSC=join(fnDirCurrent,"volumeFSC%02d.vol"%i)
copyFile(fnVol,fnForFSC) # From this point, the two half volumes may be modified
# Attenuate undershooting
if self.postSoftNeg:
removeMask=False
if fnMask=="":
fnMask=join(fnDirCurrent,"mask.vol")
self.runJob("xmipp_transform_mask","-i %s --mask circular %d --create_mask %s"%(fnVol1,-volXdim/2,fnMask),numberOfMpi=1)
removeMask=True
fnFsc=join(fnDirCurrent,"fscSoft.xmd")
self.runJob('xmipp_resolution_fsc','--ref %s -i %s -o %s --sampling_rate %f'%(fnVol1,fnVol2,fnFsc,TsCurrent),numberOfMpi=1)
self.runJob("xmipp_transform_filter","-i %s --softnegative %s %s %f %f"%(fnVol1,fnMask,fnFsc,TsCurrent,self.postSoftNegK),numberOfMpi=1)
self.runJob("xmipp_transform_filter","-i %s --softnegative %s %s %f %f"%(fnVol2,fnMask,fnFsc,TsCurrent,self.postSoftNegK),numberOfMpi=1)
cleanPath(fnFsc)
if removeMask:
cleanPath(fnMask)
fnMask=""
# Difference evaluation and production of a consensus average
if self.postDifference:
fnRootRestored=join(fnDirCurrent,"volumeRestored")
args='--i1 %s --i2 %s --oroot %s --difference 2 2'%(fnVol1,fnVol2,fnRootRestored)
if fnMask!="":
args+=" --mask binary_file %s"%fnMask
if self.useGpu:
self.runJob('xmipp_cuda_volume_halves_restoration', args, numberOfMpi=1)
else:
self.runJob('xmipp_volume_halves_restoration',args,numberOfMpi=1)
self.runJob("xmipp_image_convert","-i %s_avgDiff.vol -o %s -t vol"%(fnRootRestored,fnVolAvg),numberOfMpi=1)
cleanPath("%s_avgDiff.vol"%fnRootRestored)
cleanPath("%s_restored1.vol"%fnRootRestored)
cleanPath("%s_restored2.vol"%fnRootRestored)
# Recalculate the average after alignment and denoising
if not exists(fnVolAvg):
self.runJob('xmipp_image_operate','-i %s --plus %s -o %s'%(fnVol1,fnVol2,fnVolAvg),numberOfMpi=1)
self.runJob('xmipp_image_operate','-i %s --mult 0.5'%fnVolAvg,numberOfMpi=1)
# if fnMask!="":
# self.runJob("xmipp_image_operate","-i %s --mult %s"%(fnVolAvg,fnMask),numberOfMpi=1)
self.runJob('xmipp_image_header','-i %s --sampling_rate %f'%(fnVolAvg,TsCurrent),numberOfMpi=1)
[docs] def cleanDirectory(self, iteration):
fnDirCurrent=self._getExtraPath("Iter%03d"%iteration)
if self.saveSpace:
fnGlobal=join(fnDirCurrent,"globalAssignment")
fnLocal=join(fnDirCurrent,"localAssignment")
if exists(fnGlobal):
cleanPath(join(fnGlobal,"images.stk"))
for i in range(1,3):
if exists(fnGlobal):
cleanPath(join(fnGlobal,"images%02d.xmd"%i))
cleanPath(join(fnGlobal,"volumeRef%02d.vol"%i))
if exists(fnLocal):
cleanPath(join(fnLocal,"images%02d.xmd"%i))
cleanPath(join(fnLocal,"images.stk"))
cleanPath(join(fnLocal,"anglesCont%02d.stk"%i))
cleanPath(join(fnLocal,"anglesDisc%02d.xmd"%i))
cleanPath(join(fnLocal,"volumeRef%02d.vol"%i))
fnCorrectedImages=join(fnDirCurrent,"images_corrected%02d.stk"%i)
if exists(fnCorrectedImages) and iteration!=self.firstIteration+self.numberOfIterations.get()-1:
cleanPath(fnCorrectedImages) # Delete corrected images except for the last iteration
#if self.weightResiduals:
# cleanPath(join(fnLocal,"covariance%02d.stk"%i))
# cleanPath(join(fnLocal,"residuals%02i.stk"%i))
#--------------------------- INFO functions --------------------------------------------
def _validate(self):
errors = []
if isinstance(self.inputVolumes.get(),SetOfVolumes) and self.inputVolumes.get().getSize()!=2:
errors.append("The set of input volumes should have exactly 2 volumes")
if self.postSymmetryWithinMask and not self.postSymmetryWithinMaskMask.hasValue():
errors.append("Symmetrize within mask requires a mask")
if not self.doContinue and not self.inputParticles.hasValue():
errors.append("You must provide input particles")
if not self.doContinue and self.inputParticles.hasValue() and \
self.alignmentMethod.get()==self.LOCAL_ALIGNMENT and not self.inputParticles.get().hasAlignmentProj():
errors.append("If the first iteration is local, then the input particles must have an alignment")
if self.useGpu.get() and not isXmippCudaPresent():
errors.append("You have asked to use GPU, but I cannot find the Xmipp GPU programs in the path")
return errors
def _warnings(self):
warnings = []
return warnings
def _summary(self):
summary = []
summary.append("Symmetry: %s" % self.symmetryGroup.get())
summary.append("Number of iterations: "+str(self.numberOfIterations))
if self.alignmentMethod==self.GLOBAL_ALIGNMENT:
summary.append("Global alignment, max shift=%f"%self.angularMaxShift.get())
else:
auxStr="Local alignment, refining: "
if self.contShift:
auxStr+="shifts "
if self.contScale:
auxStr+="scale "
if self.contAngles:
auxStr+="angles "
if self.contGrayValues:
auxStr+="gray "
if self.contDefocus:
auxStr+="defocus"
summary.append(auxStr)
auxStr="Weights: "
if self.weightSSNR:
auxStr+="SSNR "
if self.weightContinuous and self.alignmentMethod==self.LOCAL_ALIGNMENT:
auxStr+="Continuous "
if self.weightJumper:
auxStr+="Jumper"
summary.append(auxStr)
if self.postSymmetryWithinMask:
summary.append("Symmetrizing within mask: "+self.postMaskSymmetry)
if self.postSymmetryHelical:
summary.append("Looking for helical symmetry")
return summary
def _methods(self):
strline = ''
if hasattr(self, 'outputVolume') or True:
strline += 'We processed %d particles from %s ' % (self.inputParticles.get().getSize(),
self.getObjectTag('inputParticles'))
if self.inputVolumes.get() is not None:
strline += 'using %s as reference and Xmipp highres procedure. ' % (self.getObjectTag('inputVolumes'))
if self.symmetryGroup!="c1":
strline+="We imposed %s symmetry. "%self.symmetryGroup
strline += "We performed %d iterations of "%self.numberOfIterations.get()
if self.alignmentMethod==self.GLOBAL_ALIGNMENT:
strline+=" global alignment (max. shift=%f)"%self.angularMaxShift
else:
strline+=" local alignment, refining "
if self.contShift:
strline+="shifts "
if self.contScale:
strline+="scale "
if self.contAngles:
strline+="angles "
if self.contGrayValues:
strline+="gray "
if self.contDefocus:
strline+="defocus"
strline+=". "
if self.weightSSNR or (self.weightContinuous and self.alignmentMethod==self.LOCAL_ALIGNMENT) or self.weightJumper:
strline+="For reconstruction, we weighted the images according to "
if self.weightSSNR:
strline+="their SSNR "
if self.weightContinuous and self.alignmentMethod==self.LOCAL_ALIGNMENT:
strline+=", their correlation in the continuous alignment "
if self.weightJumper:
strline+=", and their angular stability"
strline+=". "
if self.postAdHocMask.hasValue():
strline+="We masked the reconstruction with %s. "%self.getObjectTag('postAdHocMask')
if self.postSymmetryWithinMask:
strline+="We imposed %s symmetry within the mask %s. "%(self.postSymmetryWithinMaskType.get(),self.getObjectTag('postSymmetryWithinMaskMask'))
if self.postSymmetryHelical:
strline+="Finally, we imposed helical symmetry. "
return [strline]