#!/usr/bin/env python # autofidseed - finds fiducial seed model # # Author: David Mastronarde # # $Id: autofidseed,v 50f35d4c13fb 2025/01/09 04:38:56 mast $ # progname = 'autofidseed' prefix = 'ERROR: ' + progname + ' - ' # Clean up all the temp files before exiting def cleanup(exts, pidstr): if leavetmp: leftStr = 'All' else: leftStr = 'Some' for ext in exts: cleanlist = glob.glob(tmproot + '*' + ext + '*') cleanupFiles(cleanlist) if leavetmp or len(exts) == len(cleanExts): leftStr += ' temporary files left as ' + tmpdir + '/afs' + pidInfo + '*' if pidstr != pidInfo and leavetmp: leftStr += ' or ' + 'afs' + pidstr + '*' prnstr(leftStr) # Cleanup files, issue top message if any, and do the Imod error exit def cleanExitError(message = ''): cleanup(cleanExts + resumeExts + infoExts, pid) if message: prnstr(prefix + message) exitFromImodError(progname) # Compute a default border size if needed def setBorderSize(sizeForBorder, entered): if entered > 0: return entered border = (lowBorderPixPerK * sizeForBorder) // 1024 if sizeForBorder > lowHighBreakpoint: border = (highBorderPixPerK * lowHighBreakpoint + lowBorderPixPerK * \ (sizeForBorder - lowHighBreakpoint)) // 1024 return border # Fill in first gap in view list going out from middle in given direction def fillFirstGap(viewList, lookDir, midz): rangeLo = list(range(midz, trackStart - 1, -1)) rangeHi = list(range(midz, trackEnd + 1)) trackRange = rangeLo + rangeHi if lookDir > 0: trackRange = rangeHi + rangeLo for view in trackRange: if view not in viewList: return view else: # ELSE ON FOR exitError('Inconsistency picking seed views') # Select set of seed views given the number desired def selectSeedViews(numNeeded): angCrit = 1.9 # If all views in range are needed, just return that to avoid violating assumptions # of the logic below if numNeeded == trackEnd + 1 - trackStart: return range(trackStart, trackEnd + 1) # Get the 3 basic views at least 2 degrees apart inside the range seedViews = [midView - 1, midView, midView + 1] for ind in (0, 2): while seedViews[ind] > trackStart and seedViews[ind] < trackEnd and \ fabs(tiltIncluded[midView] - tiltIncluded[seedViews[ind]]) < angCrit: seedViews[ind] += ind - 1 # If the min tilt view is not in the list, force it in if minTiltInd not in seedViews: # Find which one is closest to the zero tilt subInd = 0 for ind in (1, 2): if fabs(tiltIncluded[seedViews[subInd]]) > fabs(tiltIncluded[seedViews[ind]]): subInd = ind # put the min tilt view in that spot and set up to walk out from there in steps seedViews[subInd] = minTiltInd if subInd == 1: dirList = [-1, 1] endInds = [trackStart, trackEnd] elif subInd == 0: dirList = [1, 1] endInds = [trackEnd - 1, trackEnd] else: dirList = [-1, -1] endInds = [trackStart + 1, trackStart] # Take two steps in the listed directions and with the proper end point for each # step and find the next view at the correct separation for loop in (0, 1): indDir = dirList[loop] ind = subInd + indDir seedViews[ind] = seedViews[subInd] + indDir while seedViews[ind] != endInds[loop] and \ fabs(tiltIncluded[seedViews[subInd]] - \ tiltIncluded[seedViews[ind]]) < angCrit: seedViews[ind] += indDir if dirList[0] == dirList[1]: subInd += indDir checkDir = -1 endView = seedViews[0] divided = [0, 0] midZ = seedViews[1] for loop in (0, 1): if len(seedViews) >= numNeeded: return seedViews viewTry = endView + checkDir # Try to extend by another 2 degrees while viewTry > trackStart and viewTry < trackEnd and \ fabs(tiltIncluded[endView] - tiltIncluded[viewTry]) < angCrit: viewTry += checkDir # Find insertion between the two with most balanced intervals minView = -1 if endView - checkDir != midZ: checkLo = min(endView - checkDir, midZ + checkDir) checkHi = max(endView - checkDir, midZ + checkDir) + 1 minDiff = 1000. for view in range(checkLo, checkHi): interval1 = fabs(tiltIncluded[endView] - tiltIncluded[view]) interval2 = fabs(tiltIncluded[midZ] - tiltIncluded[view]) diff = fabs(interval1 - interval2) if diff < minDiff: minDiff = diff minView = view minInterval = min(interval1, interval2) # If there was any, and either the extension is invalid or the minimum interval # by dividing this range is larger than the extension interval, divide range if minView >= 0: if viewTry < trackStart or viewTry > trackEnd or minInterval > \ 1.5 * fabs(tiltIncluded[endView] - tiltIncluded[viewTry]): viewTry = minView divided[loop] = 1 # And if there is still nothing valid, just fill first gap from this direction if viewTry < trackStart or viewTry > trackEnd: viewTry = fillFirstGap(seedViews, checkDir, midZ) seedViews.append(viewTry) seedViews.sort() checkDir = -checkDir endView = seedViews[3] # Next select 6th and 7th ones endInd = 0 midInd = 2 for loop in (0, 1): if len(seedViews) >= numNeeded: return seedViews viewTry = -1 endView = seedViews[endInd] # If this side was divided previously, try again to go outside if divided[loop]: viewTry = endView + checkDir while viewTry > trackStart and viewTry < trackEnd and \ fabs(tiltIncluded[endView] - tiltIncluded[viewTry]) < angCrit: viewTry += checkDir # If nothing picked yet, find biggest gap in range if viewTry < trackStart or viewTry > trackEnd: fillInd = (endInd + midInd) // 2 # There must be a gap in range if midZ + 2 * checkDir != endView: # If there is no gap on one side or other of filled value, use the other side if midZ + checkDir == seedViews[fillInd]: viewTry = (seedViews[fillInd] + endView) // 2 elif seedViews[fillInd] + checkDir == endView: viewTry = (midZ + seedViews[fillInd]) // 2 # Otherwise pick the biggest angle gap elif fabs(tiltIncluded[midZ] - tiltIncluded[seedViews[fillInd]]) > \ fabs(tiltIncluded[endView] - tiltIncluded[seedViews[fillInd]]): viewTry = (midZ + seedViews[fillInd]) // 2 else: viewTry = (seedViews[fillInd] + endView) // 2 # And if there is still nothing valid, just fill first gap from this direction if viewTry < trackStart or viewTry > trackEnd: viewTry = fillFirstGap(seedViews, checkDir, midZ) seedViews.append(viewTry) seedViews.sort() checkDir = -checkDir midInd = 3 endInd = 5 return seedViews # Run imodfindbeads with current selection of views def runFindBeads(): global newBeadSize sectOpt = fmtstr('SectionsToDo {}', includedViews[seedViews[0]]) viewStr = fmtstr(' on views {}', includedViews[seedViews[0]] + 1) for i in range(1, numSeedViews): seedVw = includedViews[seedViews[i]] sectOpt += ',' + str(seedVw) if i < numSeedViews - 1: viewStr += ', ' + str(seedVw + 1) else: viewStr += ', and ' + str(seedVw + 1) comlines = ['InputImageFile ' + imageFile, 'OutputModelFile ' + tmppeak, sectOpt, fmtstr('BeadSize {}', beadSize), fmtstr('MinSpacing {}', spacing), fmtstr('LinearInterpolation {}', linear), fmtstr('StorageThreshold {}', -peakFraction), fmtstr('MinGuessNumBeads {}', guess), fmtstr('FallbackThresholds {},{}', averageFallback * numSeedViews, storageFallback * numSeedViews)] if lightBeads: comlines.append('LightBeads 1') if usingSobel: comlines.append(fmtstr('KernelSigma {:.3f}', ksigma)) if boundModel: comlines.append('AreaModel ' + boundModel) if excludeAreas: comlines.append('ExcludeInsideAreas 1') if prexgFile: comlines.append('PrealignTransformFile ' + prexgFile) comlines.append(fmtstr('ImagesAreBinned {}', binning)) if adjustSize: comlines.append('AdjustSizes 1') if findOptions: psplit = findOptions.split(' -') for opt in psplit: comlines.append(opt.lstrip('-')) prnstr('RUNNING IMODFINDBEADS' + viewStr) prnstr(' ') try: findlines = runcmd('imodfindbeads -StandardInput', comlines) except ImodpyError: cleanExitError('Running imodfindbeads' + viewStr) for l in findlines: prnstr(l, end='') if adjustSize and l.startswith('Adjusted parameters'): try: newBeadSize = float(l.split()[-1]) except ValueError: exitError('Converting new bead size to float') if ('peaks are above' not in findlines[-2] and \ 'using fallback' not in findlines[-2]) or \ 'total peaks being' not in findlines[-1]: if 'Failed to find dip' in findlines[-1]: exitError('Cannot proceed; imodfindbeads cannot identify the gold beads') exitError('Output from imodfindbeads does not end in the expected way') lsplit = findlines[-1].split() numPeaksTot = convertToInteger(lsplit[0], 'number of total peaks stored') prnstr(' ') return numPeaksTot # Find out if a higher number of seed views is needed to try to reach target def reviseNumSeedViews(peaksTot, nsViews): numPerView = peaksTot // nsViews targNum = targetNumber needed = 3 if targetDensity > 0.: targNum = targetDensity * nx * ny * 1.e-6 if numPerView < moreViewCrit5 * targNum: needed = min(numViews, 5) if numPerView < moreViewCrit7 * targNum: needed = min(numViews, 7) return needed # Test whether two of the shift-near-zero vectors match def vectorsMatch(cor1, pek1, cor2, pek2): if vecAngles[cor1][pek1] < -900 or vecAngles[cor2][pek2] < -900: return False diff = vecAngles[cor1][pek1] - vecAngles[cor2][pek2] if diff < -180: diff += 360 if diff > 180: diff -= 360 if math.fabs(diff) > maxVecAngleDiff: return False return math.fabs(vecLengths[cor1][pek1] - vecLengths[cor2][pek2]) / maxLength < \ maxVecLenDiff #### MAIN PROGRAM #### # # load System Libraries import os, sys, glob, math, shutil, re from math import fabs # # Setup runtime environment if os.getenv('IMOD_DIR') != None: IMOD_DIR = os.environ['IMOD_DIR'] if sys.platform == 'cygwin' and sys.version_info[0] > 2: IMOD_DIR = IMOD_DIR.replace('\\', '/') if IMOD_DIR[1] == ':' and IMOD_DIR[2] == '/': IMOD_DIR = '/cygdrive/' + IMOD_DIR[0].lower() + IMOD_DIR[2:] sys.path.insert(0, os.path.join(IMOD_DIR, 'pylib')) from imodpy import * addIMODbinIgnoreSIGHUP() else: sys.stdout.write(prefix + " IMOD_DIR is not defined!\n") sys.exit(1) # # load IMOD Libraries from pip import * from pysed import * # Fallbacks from ../manpages/autodoc2man 3 1 autofidseed options = ["track:TrackCommandFile:FN:", "append:AppendToSeedModel:B:", "guess:MinGuessNumBeads:I:", "spacing:MinSpacing:F:", "size:BeadSize:F:", "adjust:AdjustSizes:B:", "peak:PeakStorageFraction:F:", "find:FindBeadOptions:CH:", "views:NumberOfSeedViews:I:", "shifts:ShiftsNearZeroFraction:F:", "justshifts:JustFindShiftsNearZero:I:", "boundary:BoundaryModel:FN:", "exclude:ExcludeInsideAreas:B:", "border:BordersInXandY:IP:", "two:TwoSurfaces:B:", "number:TargetNumberOfBeads:I:", "density:TargetDensityOfBeads:F:", "ratio:MaxMajorToMinorRatio:F:", "elongated:ElongatedPointsAllowed:I:", "cluster:ClusteredPointsAllowed:I:", "lower:LowerTargetForClustered:F:", "subarea:SubareaSize:I:", "sort:SortAreasMinNumAndSize:IP:", "ignore:IgnoreSurfaceData:LI:", "drop:DropTracks:LI:", "pick:PickSeedOptions:CH:", "remove:RemoveTempFiles:I:", "output:OutputSeedModel:FN:", "info:InfoFile:FN:", "tempdir:TemporaryDirectory:FN:", "leave:LeaveTempFiles:B:", "help:usage:B:"] (opts, nonopts) = PipReadOrParseOptions(sys.argv, options, progname, 2, 0, 0) os.environ['PIP_PRINT_ENTRIES'] = '0' spacing = 0.85 optLocalPoints = 20 optLocalSize = 1000 minLocalSize = 500 minLocalPoints = 10 minPointsForLocal = 20 minSubareaNum = 50 minSubareaSize = 2500 lowBorderPixPerK = 32 highBorderPixPerK = 24 lowHighBreakpoint = 2048 minDoTilt = 40. moreViewCrit5 = 1. # Not all points are usable moreViewCrit7 = 0.5 # This is unlikely to help more so make the threshold low leavetmp = 0 targetToGuessFrac = 0.25 guessFracIfJustShifts = 0.5 targetToAverageFB = 0.33 targetToStorageFB = 2.0 minSizeForWeightScaling = 12.5 kernelSigmaDefault = 0.5 # THE DEFAULT IN BEADTRACK usableSecPeakFrac = 0.07 maxVecAngleDiff = 15. maxVecLenDiff = 0.33 newBeadSize = 0 candModel = 'clusterElong.mod' # Extensions for files that are always cleaned up, files that are left for resuming # and files that are left for user edification only. Make all extensions unique cleanExts = ['.seed', '.xyzpt', '.surf', '.mop', '.mopxf', '.moptrim'] resumeExts = ['.track', '.xyzmod', '.elong'] infoExts = ['pkmod', '.sortmod'] # Get track command file trackcom = PipGetString('TrackCommandFile', '') if not trackcom: exitError('You must enter the name of a beadtrack command file') # Compose info file name (troot, ext) = os.path.splitext(trackcom) infoSuffix = '' if troot.endswith('a'): infoSuffix = 'a' if troot.endswith('b'): infoSuffix = 'b' infoName = progname + infoSuffix + '.info' # Set root of names for temp files deftmpdir = progname + infoSuffix + '.dir' tmpdir = PipGetString('TemporaryDirectory', deftmpdir) if not os.path.exists(tmpdir): try: os.mkdir(tmpdir) except OSError: exitError('Making the temporary directory ' + tmpdir) if not (os.path.isdir(tmpdir) and os.access(tmpdir, os.W_OK)): if makeCurrentDirWritable(tmpdir): exitError('Cannot write to temporary directory ' + tmpdir) pid = str(os.getpid()) + '.' tmproot = 'afs' if tmpdir: tmproot = tmpdir + '/' + tmproot # Get alternate name,, keep track if using default infoName = PipGetString('InfoFile', infoName) defaultInfo = PipGetErrNo() # Get info file oldValid = os.path.exists(infoName) if oldValid: infoLines = readTextFile(infoName) pidInfo = optionValue(infoLines, 'PID', 0) # See if cleaning up, so we can exit now clean = PipGetInteger('RemoveTempFiles', 0) if clean != 0 and oldValid and pidInfo: cleanup(resumeExts + infoExts, pidInfo) oldValid = False if clean < 0: if tmpdir == deftmpdir: cleanlist = glob.glob(tmpdir + '/' + '*') cleanupFiles(cleanlist) try: os.rmdir(tmpdir) except OSError: prnstr('WARNING: could not remove temporary directory ' + tmpdir) sys.exit(0) # Get more options justShifts = PipGetInteger('JustFindShiftsNearZero', 0) guess = PipGetInteger('MinGuessNumBeads', 0) ifGuess = 1 - PipGetErrNo() twoSurf = PipGetBoolean('TwoSurfaces', 0) appendToSeed = PipGetBoolean('AppendToSeedModel', 0) spacing = PipGetFloat('MinSpacing', spacing) peakFraction = PipGetFloat('PeakStorageFraction', 1.0) adjustSize = PipGetBoolean('AdjustSizes', 0) findOptions = PipGetString('FindBeadOptions', '') targetNumber = PipGetInteger('TargetNumberOfBeads', 0) targetDensity = PipGetFloat('TargetDensityOfBeads', 0.) if targetNumber < 1 and targetDensity <= 0.: exitError('You must enter a positive number or density of beads as the target for ' + \ 'the seed model') boundModel = PipGetString('BoundaryModel', '') excludeAreas = PipGetBoolean('ExcludeInsideAreas', 0) if boundModel and not os.path.exists(boundModel): exitError('Boundary model ' + boundModel + ' does not exist') (minSubareaNum, minSubareaSize) = PipGetTwoIntegers('SortAreasMinNumAndSize', minSubareaNum, minSubareaSize) ifPick = 1 - PipGetErrNo() subareaSize = PipGetInteger('SubareaSize', 0) if ifPick and subareaSize: exitError('You cannot enter both -sort and -subarea') leavetmp = PipGetBoolean('LeaveTempFiles', 0) shiftsFrac = PipGetFloat('ShiftsNearZeroFraction', 0.2) if justShifts > 0: targetNumber = justShifts if shiftsFrac <= 0 or shiftsFrac >= 10.: exitError('You cannot enter -justshifts and -shifts with a value that disables ' + \ 'finding shifts') if not ifGuess: guess = int(round(guessFracIfJustShifts * justShifts)) ifGuess = 1 # Get elongated and clustered and separate the latter if it is an old type of entry elongated = PipGetInteger("ElongatedPointsAllowed", 0) clustered = PipGetInteger("ClusteredPointsAllowed", 0) clustered = min(4, max(0, clustered)) elongated = min(3, max(0, elongated)) if clustered > 1: if elongated: exitError('You cannot enter a value > 1 for -cluster if you enter -elongated') elongated = clustered - 1 clustered = 1 lowTarget = PipGetFloat("LowerTargetForClustered", 0.) maxMajorMinor = PipGetFloat("MaxMajorToMinorRatio", 0.) pickOptions = PipGetString('PickSeedOptions', '') seedFile = PipGetString('OutputSeedModel', '') (xborder, yborder) = PipGetTwoIntegers('BordersInXandY', 0, 0) numSeedViews = PipGetInteger('NumberOfSeedViews', 3) viewsEntered = 1 - PipGetErrNo() if numSeedViews < 3 or numSeedViews > 7: exitError('The number of views to use as seeds must be between 3 and 7') dropStr = PipGetString('DropTracks', '') dropList = [] if dropStr: dropList = parselist(dropStr) ignoreStr = PipGetString('IgnoreSurfaceData', '') ignoreList = [] if ignoreStr: ignoreList = parselist(ignoreStr) for ind in range(len(dropList)): dropList[ind] -= 1 for ind in range(len(ignoreList)): ignoreList[ind] -= 1 # Read com file and get critical entries rawTracklines = readTextFile(trackcom) # find beadtrack command line then end of input (from transferfid) startline = -1 endline = len(rawTracklines) for i in range(endline): line = rawTracklines[i].strip() if line.startswith('$') and 'beadtrack' in line and '-Standard' in line: startline = i + 1 elif startline > 0 and line.startswith('$'): endline = i if startline < 0: exitError('Old version of ' + trackcom + ' cannot be used; convert it by opening ' +\ 'and closing the fiducial tracking panel in etomo') tracklines = rawTracklines[startline:endline] imageFile = optionValue(tracklines, 'ImageFile', STRING_VALUE) if not seedFile: seedFile = optionValue(tracklines, 'InputSeedModel', STRING_VALUE) lightBeads = optionValue(tracklines, 'LightBeads', BOOL_VALUE) rotationArr = optionValue(tracklines, 'RotationAngle', FLOAT_VALUE); doTiltArr = optionValue(tracklines, 'MinTiltRangeToFindAngles', FLOAT_VALUE) if doTiltArr and doTiltArr[0] > minDoTilt: minDoTilt = doTiltArr[0] if not os.path.exists(imageFile): exitError("The prealigned stack " + imageFile + " does not exist yet") try: (nx, ny, nz) = getmrcsize(imageFile) except ImodpyError: exitFromImodError(progname) # Get the border sizes xborder = setBorderSize((nx + ny) // 2, xborder) yborder = setBorderSize((nx + ny) // 2, yborder) # Get the binning and prealign file binning = 1 binArr = optionValue(tracklines, 'ImagesAreBinned', 1) prexgFile = optionValue(tracklines, 'PrealignTransformFile', STRING_VALUE) if binArr: binning = binArr[0] # Get the box size and determine if finding shifts near zero tilt boxSizeArr = optionValue(tracklines, 'BoxSizeXandY', INT_VALUE, numVal = 2) pixelSize = optionValue(tracklines, 'PixelSize', FLOAT_VALUE, numVal = 1) if pixelSize == None: try: pixSpacing = getmrcpixel(imageFile) if pixSpacing > 0 and pixSpacing != 1.0: pixelSize = pixSpacing * binning / 10. except ImodpyError: pass findShiftsNearZero = False if boxSizeArr != None and pixelSize != None and shiftsFrac > 0. and shiftsFrac < 10.: findShiftsNearZero = True if justShifts > 0 and not findShiftsNearZero: exitError('Cannot find shifts near zero tilt; ' + trackcom + \ ' is missing pixel size or box size') # Get the unbinned bead size from the track com file regardless trackDiameter = optionValue(tracklines, 'BeadDiameter', 2, numVal = 1) if trackDiameter != None and trackDiameter <= 0.: exitError('The BeadDiameter entry in ' + trackcom + ' is not positive') beadSize = PipGetFloat('BeadSize', 0.) # If no -size entered and produce the value needed for findbeads if PipGetErrNo(): if trackDiameter == None: exitError('There is no BeadDiameter entry in ' + trackcom + \ '; fix this or enter a size with -size') beadSize = trackDiameter / binning sizeForWeights = beadSize # Get filtering information if sobel activated usingSobel = optionValue(tracklines, 'Sobel', 3) linear = 0 if usingSobel: kernelArr = optionValue(tracklines, 'Kernel', 2) scalableSigma = optionValue(tracklines, 'ScalableSigma', FLOAT_VALUE, numVal = 1) if kernelArr: ksigma = kernelArr[0] if ksigma >= 1.49: linear = 1 elif scalableSigma: ksigma = scalableSigma * beadSize else: ksigma = kernelSigmaDefault # Get tilt angle options and make list of tilt angles tiltFile = optionValue(tracklines, 'TiltFile', 0) tiltAngles = [] if tiltFile: tiltLines = readTextFile(tiltFile) try: for i in range(len(tiltLines)): if tiltLines[i].strip(): tiltAngles.append(float(tiltLines[i])) except ValueError: exitError('Converting lines in ' + tiltFile + ' to floating point values') else: first = optionValue(tracklines, 'FirstTiltAngle', 2) increment = optionValue(tracklines, 'TiltIncrement', 2) if not (first and increment): exitError('The track command file must have either a tilt angle file or starting '+\ 'and increment tilt angles') for i in range(nz): tiltAngles.append(first[0] + i * increment[0]) # Apply whatever angle offset is in track entry angleOffset = optionValue(tracklines, 'AngleOffset', FLOAT_VALUE, numVal = 1) if angleOffset: for ind in range(len(tiltAngles)): tiltAngles[ind] += angleOffset # Get the exclude list, which is numbered from 1 skipListStr = optionValue(tracklines,'SkipViews', 0) excludeList = [] if skipListStr: excludeList = parselist(skipListStr) # Make a list of included views and get the number of them numViews = len(tiltAngles) includedViews = [] tiltIncluded = [] for iv in range(numViews): if iv + 1 not in excludeList: includedViews.append(iv) tiltIncluded.append(tiltAngles[iv]) numIncluded = len(includedViews) # Find minimum tilt view and highest tilt minTiltInd = 0 highestTilt = 0. if numIncluded < 3: exitError('There must be at least 3 views in the tilt series and not in a skip list') numSeedViews = min(numIncluded, numSeedViews) for iv in range(numViews): highestTilt = max(highestTilt, fabs(tiltAngles[iv])) for incl in range(numIncluded): if fabs(tiltIncluded[minTiltInd]) + 0.1 >= fabs(tiltIncluded[incl]): minTiltInd = incl; # Get view range for tracking: do 11 views unless range is > 20 deg, or 9 views unless # range is still > 20 deg; or 7 views for viewInc in (10, 8, 6): trackStart = max(0, minTiltInd - viewInc // 2) trackEnd = min(numIncluded - 1, trackStart + viewInc) trackStart = max(0, trackEnd - viewInc) if fabs(tiltIncluded[trackStart] - tiltIncluded[trackEnd]) <= 20.1: break # For very fine increments, increase the view range to have at least 7.5 degree track if fabs(tiltIncluded[trackStart] - tiltIncluded[trackEnd]) < 7.4: for viewInc in range(12, 80, 2): trackStart = max(0, minTiltInd - viewInc // 2) trackEnd = min(numIncluded - 1, trackStart + viewInc) trackStart = max(0, trackEnd - viewInc) if fabs(tiltIncluded[trackStart] - tiltIncluded[trackEnd]) >= 7.4: break # midView, trackStart, trackEnd and seed lists will all be included view indexes # so need to take includedViews[view] to get true Z midView = (trackEnd + trackStart) // 2 # Convert a density to a target number # If there is a boundary model, find out the total area from imodfindbeads if targetDensity > 0.: targetDensity *= binning * binning totalArea = nx * ny * 1.e-6 if boundModel: comlines = ['InputImageFile ' + imageFile, 'AreaModel ' + boundModel, 'QueryAreaOnSection ' + str(includedViews[minTiltInd])] if excludeAreas: comlines.append('ExcludeInsideAreas 1') try: findlines = runcmd('imodfindbeads -StandardInput', comlines) for l in findlines: ind = l.find('=') if ind > 0 and 'Area (megapixels)' in l: totalArea = float(l[ind + 1:]) break else: exitError('Cannot find area being analyzed from imodfindbeads output') except ImodpyError: cleanExitError('Running imodfindbeads to determine area being analyzed') except ValueError: exitError('Converting area being analyzed from imodfindbeads output') targetNumber = targetDensity * totalArea # Now it is possible to convert to get fallbacks for the guess and thresholds if not ifGuess: guess = max(1, int(round(targetNumber * targetToGuessFrac))) averageFallback = max(1, int(round(targetNumber * targetToAverageFB))) storageFallback = max(1, int(round(targetNumber * targetToStorageFB))) trackMtime = int(os.stat(trackcom).st_mtime) imageMtime = int(os.stat(imageFile).st_mtime) beadStr = str(beadSize) spacingStr = str(spacing) if boundModel: boundMtime = int(os.stat(boundModel).st_mtime) (comdir, comname) = os.path.split(trackcom) comSaved = tmpdir + '/' + comname # Process info file fully now if oldValid: infoLines = readTextFile(infoName) pidInfo = optionValue(infoLines, 'PID', 0) trackInfo = optionValue(infoLines, 'TrackCom', 0) trackTimeArr = optionValue(infoLines, 'TrackTime', 1) imageTimeArr = optionValue(infoLines, 'ImageTime', 1) beadInfo = optionValue(infoLines, 'BeadSize', 0) guessInfoArr = optionValue(infoLines, 'MinGuess', 1) peakFracInfoArr = optionValue(infoLines, 'PeakFraction', 2) spacingInfo = optionValue(infoLines, 'Spacing', 0) twoSurfInfo = optionValue(infoLines, 'TwoSurf', 3) boundInfo = optionValue(infoLines, 'BoundFile', 0) boundTimeArr = optionValue(infoLines, 'BoundTime', 1) numPeakArr = optionValue(infoLines, 'NumPeaks', 1) numSeedArr = optionValue(infoLines, 'NumViews', 1) viewEnterInfo = optionValue(infoLines, 'ViewsEntered', 3) adjSizeArr = optionValue(infoLines, 'AdjustSizes', 2) zeroFracArr = optionValue(infoLines, 'ZeroShiftsFrac', 2) zeroShifts = optionValue(infoLines, 'ShiftsNearZeroTilt', 0) # Need to check contents of track.com if it doesn't match time trackMatches = False if trackInfo == trackcom and trackTimeArr and trackTimeArr[0] != trackMtime: if os.path.exists(comSaved): savedLines = readTextFile(comSaved) if len(savedLines) == len(rawTracklines): for ind in range(len(savedLines)): if savedLines[ind] != rawTracklines[ind]: break else: trackMatches = True # Now ready to evaluate validity oldValid = justShifts <= 0 and pidInfo and trackInfo == trackcom and \ trackTimeArr and (trackTimeArr[0] == trackMtime or trackMatches) and \ imageTimeArr and imageTimeArr[0] == imageMtime and \ beadInfo == beadStr and \ guessInfoArr and guessInfoArr[0] == guess and \ peakFracInfoArr and peakFracInfoArr[0] == peakFraction and \ spacingInfo == spacingStr and \ twoSurf <= twoSurfInfo and \ (numSeedArr and numPeakArr and \ ((not viewsEntered and not viewEnterInfo) or \ (viewsEntered and numSeedArr[0] == numSeedViews))) and \ adjSizeArr and int(round(adjSizeArr[0])) == adjustSize and \ zeroFracArr and zeroFracArr[0] == shiftsFrac and \ ((boundInfo and boundInfo == boundModel and boundTimeArr and \ boundMtime == boundTimeArr[0]) or (not boundInfo and boundModel == '')) if oldValid and not viewsEntered: needed = reviseNumSeedViews(numPeakArr[0], numSeedArr[0]) if needed > numSeedArr[0]: oldValid = False numSeedViews = needed # Passed that test (!), now make sure all the required files are still there if oldValid: for ind in range(numSeedArr[0]): if not oldValid: break for ext in resumeExts: if ext == resumeExts[1] and not twoSurfInfo: continue if not os.path.exists(tmproot + pidInfo + str(ind) + ext): oldValid = False break # If not resuming and there was an old PID, clean up old files if pidInfo and not oldValid: leaveSave = leavetmp leavetmp = 0 cleanup(resumeExts + infoExts, pidInfo) leavetmp = leaveSave # Adjust number of seed views if resuming if oldValid: numSeedViews = numSeedArr[0] if len(adjSizeArr) > 1 and adjSizeArr[1] > 0: sizeForWeights = adjSizeArr[1] prnstr('Adjusted parameters in previous run for new bead size of ' + str(adjSizeArr[1]) + ' [AFS2]') if zeroShifts != None: prnstr('Adjusted parameters in previous run with shifts per view entry of ' + zeroShifts + ' [AFS4]') # If making a new info file with default name, make sure current directory is writable if defaultInfo and not oldValid and not os.access('.', os.W_OK): exitError("You cannot write to the current directory; use -info to make the info " "file elsewhere") # Check validity of drop and ignore lists if (dropList or ignoreList): if not oldValid: exitError("You cannot list tracks to drop or ignore unless resuming with " + \ "existing tracks") numDrop = 0 numIgnore = 0 for i in range(numSeedViews): if i in dropList: numDrop += 1 if i in ignoreList: numIgnore += 1 if numDrop == numSeedViews: exitError('The list of tracks to drop includes all tracks') if numIgnore == numSeedViews: exitError('The list of tracks to ignore surface data from includes all tracks') # For testing that incredible selection logic #for i in range(3,8): # seedViews = selectSeedViews(i) # print seedViews #sys.exit(0) seedViews = selectSeedViews(numSeedViews) skipList = '' trackVwStart = includedViews[trackStart] trackVwEnd = includedViews[trackEnd] if trackVwStart == 1: skipList = '1' elif trackVwStart > 1: skipList = '1-' + str(trackVwStart) if trackVwEnd < nz - 1 and skipList != '': skipList += ',' if trackVwEnd == nz - 2: skipList += str(nz) elif trackVwEnd < nz - 2: skipList += fmtstr('{}-{}', trackVwEnd + 2, nz) if skipList == '': skipList = '0' # Ready to find the beads if not oldValid: pidInfo = pid tmppeak = tmproot + pid + infoExts[0] numPeaksTot = runFindBeads() # Unless number of views was specified, check if there are not enough and retrack with # more views if not viewsEntered: numSeedViews = reviseNumSeedViews(numPeaksTot, numSeedViews) if numSeedViews > 3: prnstr('REDOING IMODFINDBEADS with more views to try to get more points') seedViews = selectSeedViews(numSeedViews) numPeaksTot = runFindBeads() # Skip views to use views at fixed intervals if there are lots of them # But include the seed views so every track is the same skipInterval = (trackEnd + 1 - trackStart) // 11 if skipInterval > 1: lastIncluded = trackStart for view in range(trackStart + 1, trackEnd): if view - lastIncluded < skipInterval and view not in seedViews: skipList += ',' + str(includedViews[view] + 1) else: lastIncluded = view # Finding shifts near zero: get the range of Z values for imodmop shiftsLine = '' adjustSed = [] if findShiftsNearZero and not oldValid: # Find seed view nearest zero nearZero = -1 minDiff = 100 for ind in range(len(seedViews)): diff = math.fabs(seedViews[ind] - minTiltInd) if nearZero < 0 or diff < minDiff: minDiff = diff nearZero = ind if not nearZero: nearZero += 1 if nearZero >= len(seedViews) - 1: nearZero -= 1 viewNearZero = seedViews[nearZero] viewBelowZero = seedViews[nearZero - 1] viewAboveZero = seedViews[nearZero + 1] minZ = includedViews[viewBelowZero] midZ = includedViews[viewNearZero] maxZ = includedViews[viewAboveZero] diam = int(1.5 * beadSize) mopFile = tmproot + pid + 'mop' (nx,ny,nz,mode,px,py,pz,ox,oy,oz,dmin,dmax,dmean) = getmrc(imageFile, doAll = True) mopcmd = fmtstr('imodmop -tube 1 -diam {} -fv {} -zmin {},{} -planar {} "{}" {}', diam, dmean, minZ, maxZ, tmppeak, imageFile, mopFile) try: runcmd(mopcmd) except ImodpyError: cleanExitError('Running imodmop to isolate bead component of images') mopTrim = tmproot + pid + 'moptrim' newstcmd = fmtstr('newstack -sec {},{},{} {} {}', 0, midZ - minZ, maxZ - minZ, mopFile, mopTrim) try: runcmd(newstcmd) except ImodpyError: cleanExitError('Running Newstack to extract imodmop output for correlation') increment = math.fabs(tiltIncluded[viewAboveZero] - tiltIncluded[viewBelowZero]) / 2. maxSwing = 800. * math.sin(math.radians(increment)) length = int(2 * maxSwing / (binning * pixelSize)) xcorrCom = ['InputFile ' + mopTrim, 'OutputFile ' + tmproot + pid + 'mopxf', fmtstr('TiltAngles {},{},{}', tiltIncluded[viewBelowZero], tiltIncluded[viewNearZero], tiltIncluded[viewAboveZero]), 'FilterSigma1 0.03', 'FilterRadius2 0.25', 'FilterSigma2 0.05', fmtstr('SecondPeakBoxSize {},{}', length, diam)] if rotationArr: xcorrCom.append(fmtstr('RotationAngle {}', rotationArr[0])) try: xcorrLines = runcmd('tiltxcorr -StandardInput', xcorrCom) except ImodpyError: cleanExitError('Running Tiltxcorr to determine bead shifts near zero tilt') shifts = [[0] * 6, [0] * 6] vecAngles = [[-999, -999], [-999, -999]] vecLengths = [[0, 0], [0, 0]] maxLength = 0. foundViews = [] for line in xcorrLines: if line.startswith('View'): nums = re.sub('[a-zA-DF-Z,]', '', line) numSplit = nums.split() foundViews.append(numSplit[0]) # Set up how to deal with the sign of the output based on order of views # Tilt angle is irrelevant, the shifts are per view divByViews = [viewAboveZero - viewNearZero, viewNearZero - viewBelowZero] if foundViews == ['3', '1']: sign = -1. keepSign = False elif foundViews == ['2', '3']: sign = -1. keepSign = True divByViews = [viewNearZero - viewBelowZero, viewAboveZero - viewNearZero] elif foundViews == ['2', '1']: sign = 1. keepSign = True else: foundViews = None prnstr('WARNING: View numbers in output from tiltxcorr finding shifts ' +\ 'near zero do not have expected values'); shInd = 0 for line in xcorrLines: if foundViews and line.startswith('View'): nums = re.sub('[a-zA-DF-Z,]', '', line) numSplit = nums.split() # Convert numbers, apply sign to shift for ind in range(1, min(len(numSplit), 7)): try: shifts[shInd][ind - 1] = float(numSplit[ind]) if ind != 3 and ind != 6: shifts[shInd][ind - 1] *= sign / divByViews[shInd] except ValueError: cleanup(cleanExts + resumeExts + infoExts, pid) exitError('Converting number to float in tiltxcorr output line ', line) # Get angle and length for primary, then for secondary if it is strong enough vecAngles[shInd][0] = math.degrees(math.atan2(shifts[shInd][1],shifts[shInd][0])) vecLengths[shInd][0] = math.sqrt(shifts[shInd][0]**2 + shifts[shInd][1]**2) maxLength = max(maxLength, vecLengths[shInd][0]) if shifts[shInd][5] > usableSecPeakFrac * shifts[shInd][2]: vecAngles[shInd][1] = math.degrees(math.atan2(shifts[shInd][4], shifts[shInd][3])) vecLengths[shInd][1] = math.sqrt(shifts[shInd][3]**2 + shifts[shInd][4]**2) maxLength = max(maxLength, vecLengths[shInd][1]) shInd += 1 if not keepSign: sign = 1. if foundViews != None and (shifts[0][2] == 0 or shifts[1][2] == 0): maxLength = 0. prnstr('WARNING: Output from tiltxcorr finding shifts near zero is not usable') # Proceed only if both primaries are within a reasonable distance of 0 and threshold # met if maxLength < length / 2. and maxLength * binning > shiftsFrac * boxSizeArr[0]: # If primary angles match, set first shift and look for second if vectorsMatch(0, 0, 1, 0) or vectorsMatch(0, 1, 1, 1): if vectorsMatch(0, 0, 1, 0): shiftsLine = fmtstr('{:.1f},{:.1f}', binning * (shifts[0][0] + shifts[1][0]) / 2., binning * (shifts[0][1] + shifts[1][1]) / 2.) if vectorsMatch(0, 1, 1, 1): if shiftsLine: shiftsLine += ',' shiftsLine += fmtstr('{:.1f},{:.1f}', binning * (shifts[0][3] + shifts[1][3]) / 2., binning * (shifts[0][4] + shifts[1][4]) / 2.) # Otherwise look for crossed matches and use them if found elif vectorsMatch(0, 0, 1, 1) or vectorsMatch(1, 0, 0, 1): if vectorsMatch(0, 0, 1, 1): shiftsLine = fmtstr('{:.1f},{:.1f}', binning * (shifts[0][0] + shifts[1][3]) / 2., binning * (shifts[0][1] + shifts[1][4]) / 2.) if vectorsMatch(1, 0, 0, 1): if shiftsLine: shiftsLine += ',' shiftsLine += fmtstr('{:.1f},{:.1f}', binning * (shifts[1][0] + shifts[0][3]) / 2., binning * (shifts[1][1] + shifts[0][4]) / 2.) if shiftsLine: adjustSed += sedDelAndAdd('ShiftsNearZeroTilt', shiftsLine, 'BeadDiameter', delim = '?') prnstr('Tracking parameters adjusted with shifts per view of ' + shiftsLine + ' [AFS3]') if justShifts > 0: newinfo = ['PID ' + pid] writeTextFile(infoName, newinfo) cleanup(cleanExts, pid) sys.exit(0) # Adjust tracking parameters if the bead size changed by 5% as the man page said # existing beadSize value here is binned if not oldValid and adjustSize and newBeadSize and \ math.fabs(newBeadSize - beadSize) > 0.05 * beadSize: unbinnedNewSize = newBeadSize * binning adjustSed.append(fmtstr('?^BeadDiameter?s?[ ].*? {}?', unbinnedNewSize)) prnstr(fmtstr('Tracking parameters adjusted for new unbinned bead size of {:.2f}' + \ ' [AFS1]', unbinnedNewSize)) sizeForWeights = newBeadSize minDiam = optionValue(tracklines, 'MinDiamForParamScaling', 2, numVal = 1) # If above the diameter for scaling in the com file, rescale the parameters if minDiam and trackDiameter: oldScale = 1. newScale = 1. if trackDiameter > minDiam: oldScale = trackDiameter / minDiam if newBeadSize * binning > minDiam: newScale = newBeadSize * binning / minDiam scale = newScale / oldScale if scale != 1.: for option in ('DistanceRescueCriterion', 'PostFitRescueResidual', 'MaxRescueDistance'): crit = optionValue(tracklines, option, 2, numVal = 1) if crit: adjustSed.append(fmtstr('?^{}?s?[ ].*? {:.2f}?', option, crit * scale)) critArr = optionValue(tracklines, 'DeletionCriterionMinAndSD', 2) if critArr and len(critArr) > 1: adjustSed.append(fmtstr( '?^DeletionCriterionMinAndSD?s?[ ].*? {:.3f},{:.2f}?', critArr[0] * scale, critArr[1])) boxArr = optionValue(tracklines, 'BoxSizeXandY', 1) if boxArr and len(boxArr) > 1: newXbox = 2 * int(round(boxArr[0] * scale / 2.)) newYbox = 2 * int(round(boxArr[1] * scale / 2.)) adjustSed.append(fmtstr('?^BoxSizeXandY?s?[ ].*? {},{}?', newXbox, newYbox)) # Put out a track file with adjusted parameters as well as use them here if adjustSed: (adjRoot, adjExt) = os.path.splitext(trackcom) adjFile = adjRoot + '_adjusted' + adjExt pysed(adjustSed, rawTracklines, adjFile, False, '?') # Extract the seed models tmpseed = [] tmptrack = [] tmpsurf = [] tmpxyzpt = [] tmpxyzmod = [] tmpelong = [] tmpsortmod = [] for ind in range(numSeedViews): tmpseed.append(tmproot + pid + str(ind) + cleanExts[0]) tmptrack.append(tmproot + pidInfo + str(ind) + resumeExts[0]) tmpxyzpt.append(tmproot + pid + str(ind) + cleanExts[1]) tmpsurf.append(tmproot + pid + str(ind) + cleanExts[2]) tmpxyzmod.append(tmproot + pidInfo + str(ind) + resumeExts[1]) tmpelong.append(tmproot + pidInfo + str(ind) + resumeExts[2]) tmpsortmod.append(tmproot + pidInfo + str(ind) + infoExts[1]) seedVw = includedViews[seedViews[ind]] viewStr = str(seedVw + 1) if not oldValid: try: cliplines = runcmd(fmtstr('clipmodel -keep -zmin {0},{0} "{1}" "{2}"', seedVw, tmppeak, tmpseed[ind])) except ImodpyError: cleanExitError('Running clipmodel to extract seed for view ' + viewStr) numPeaks = -1 for l in cliplines: if 'Number of points' in l: lsplit = l.split() numPeaks = convertToInteger(lsplit[-1], 'number of points remaining') if numPeaks < (numPeaksTot / numSeedViews) / 3: exitError(fmtstr('Found only {} out of {} points on view {}', numPeaks, numPeaksTot, seedVw + 1)) break else: cleanExitError('Clipmodel did not give expected output when extracting ' +\ 'points for view ' + viewStr) localSize = 1000 localTrack = 0 if numPeaks >= minPointsForLocal: localTrack = 1 density = numPeaks / float(nx * ny) localSize = max(minLocalSize, \ int((optLocalPoints * optLocalSize / density)**0.333)) if density * localSize**2 < minLocalPoints: localSize = int(math.sqrt(minLocalPoints / density)) if localSize > 0.8 * nx and localSize > 0.8 * ny: localTrack = 0 sedcom = ['?^InputSeedModel?s?[ ].*? ' + tmpseed[ind] + '?', '?^OutputModel?s?[ ].*? ' + tmptrack[ind] + '?', '?^SkipViews?d', '?^RoundsOfTracking?s?[ ].*? 2?', fmtstr('?^LocalAreaTracking?s?[ ].*? {}?', localTrack), fmtstr('?^LocalAreaTargetSize?s?[ ].*? {}?', localSize), fmtstr('?^MinTiltRangeToFindAngles?s?[ ].*? {}?', minDoTilt), '?^OutputModel?a?ElongationOutputFile ' + tmpelong[ind] + '?', '?^OutputModel?a?SkipViews ' + skipList + '?'] if twoSurf: sedcom.append('?^OutputModel?a?XYZOutputFile ' + tmpxyzpt[ind] + '?') if adjustSed: sedcom += adjustSed sedlines = pysed(sedcom, tracklines, None, False, '?') prnstr('RUNNING BEADTRACK with seed from view ' + viewStr, flush = True) try: tracklog = runcmd('beadtrack -StandardInput', sedlines) except ImodpyError: cleanExitError('Running beadtrack with seed from view ' + viewStr) # Convert point list to model file if twoSurf: pointcom = fmtstr('point2model -values -1 -sphere {} "{}" "{}"', int((beadSize + 2.) / 2.), tmpxyzpt[ind], tmpxyzmod[ind]) try: runcmd(pointcom) except ImodpyError: cleanExitError('Running point2model with XYZ data from view ' + viewStr) # The info file can now be saved for a new run if not oldValid: newinfo = ['PID ' + pid, 'TrackCom ' + trackcom, 'TrackTime ' + str(trackMtime), 'ImageTime ' + str(imageMtime), 'BeadSize ' + beadStr, 'MinGuess ' + str(guess), 'PeakFraction ' + str(peakFraction), 'Spacing ' + spacingStr, 'TwoSurf ' + str(twoSurf), 'NumPeaks ' + str(numPeaksTot), 'NumViews ' + str(numSeedViews), 'ViewsEntered ' + str(viewsEntered), 'ZeroShiftsFrac ' + str(shiftsFrac), fmtstr('AdjustSizes {} {}', adjustSize, newBeadSize)] if boundModel: newinfo += ['BoundFile ' + boundModel, 'BoundTime ' + str(boundMtime)] if shiftsLine: newinfo += ['ShiftsNearZeroTilt ' + shiftsLine] writeTextFile(infoName, newinfo) cleanupFiles(comSaved) try: shutil.copyfile(trackcom, comSaved) except Exception: prnstr('WARNING: autofidseed - failed to copy ' + trackcom + ' to ' + tmpdir) # Sort the beads onto two surfaces for each model if twoSurf: for ind in range(numSeedViews): if ind in dropList or ind in ignoreList: continue viewStr = str(includedViews[seedViews[ind]] + 1) sortcom = ['TextFileWithSurfaces ' + tmpsurf[ind], 'ValuesToRestrainSorting 1', 'FlipYandZ 0', 'InputFile ' + tmpxyzmod[ind], 'OutputFile ' + tmpsortmod[ind]] if subareaSize: sortcom.append('SubareaSize ' + str(subareaSize)) else: sortcom.append(fmtstr('PickAreasMinNumAndSize {} {}', minSubareaNum, minSubareaSize)) prnstr(' ') prnstr('RUNNING SORTBEADSURFS with XYZ positions from view ' + viewStr) try: runcmd('sortbeadsurfs -Stand', sortcom, 'stdout') except ImodpyError: cleanup(cleanExts, pid) prnstr(prefix + 'Running sortbeadsurfs with XYZ values from view ' + viewStr) exitFromImodError(progname) # Now run pickbestseed; for this we need the seed view with minimum tilt lowestTilt = 1000. for ind in range(numSeedViews): if ind not in dropList and fabs(tiltIncluded[seedViews[ind]]) < fabs(lowestTilt): zeroView = seedViews[ind] lowestTilt = tiltIncluded[zeroView] pickcom = ['OutputSeedModel ' + seedFile, fmtstr('ImageSizeXandY {} {}', nx, ny), fmtstr('BordersInXandY {} {}', xborder, yborder), fmtstr('BeadSize {}', beadSize), fmtstr('MiddleZvalue {}', includedViews[zeroView]), fmtstr('CandidateModel {}/{}', tmpdir, candModel)] for ind in range(numSeedViews): if ind not in dropList: pickcom += ['TrackedModel ' + tmptrack[ind], 'ElongationFile ' + tmpelong[ind], fmtstr('SeedZvalue {}', includedViews[seedViews[ind]])] if ind not in ignoreList: pickcom.append('SurfaceFile ' + tmpsurf[ind]) if twoSurf: pickcom.append('TwoSurfaces') if appendToSeed: pickcom.append('AppendToSeedModel') if boundModel: pickcom.append('BoundaryModel ' + boundModel) if excludeAreas: pickcom.append('ExcludeInsideAreas 1') if clustered: pickcom.append(fmtstr('ClusteredPointsAllowed {}', clustered)) if elongated: pickcom.append(fmtstr('ElongatedPointsAllowed {}', elongated)) if (clustered or elongated) and lowTarget: pickcom.append(fmtstr('LowerTargetForClustered {}', lowTarget)) if rotationArr: pickcom.append(fmtstr('RotationAngle {}', rotationArr[0])) pickcom.append(fmtstr('HighestTiltAngle {}', highestTilt)) # Process the additional options if any weightEntered = False if pickOptions: psplit = pickOptions.split(' -') for opt in psplit: pickcom.append(opt.lstrip('-')) optSplit = opt.lstrip('-').split() if 'WeightsForScore'.startswith(optSplit[0]) or 'weights'.startswith(optSplit[0]): weightEntered = True if targetDensity > 0.: pickcom.append(fmtstr('TargetDensityOfBeads {}', targetDensity)) else: pickcom.append(fmtstr('TargetNumberOfBeads {}', targetNumber)) # Unless weights were entered, adjust the weights for large beads so that distance-based # measures have appropriately less weight if not weightEntered and sizeForWeights > minSizeForWeightScaling: scale = minSizeForWeightScaling / sizeForWeights pickcom.append(fmtstr('WeightsForScore 1.,1.,{:.4f},{:.4f}', scale, scale)) prnstr(' ') prnstr('RUNNING PICKBESTSEED') prnstr(' ') try: pickLines = runcmd('pickbestseed -StandardInput', pickcom) for l in pickLines: prnstr(l, end='') if twoSurf and maxMajorMinor > 0.: for lind in range(len(pickLines) - 1, -1, -1): if pickLines[lind].startswith('Final:'): lsplit = pickLines[lind].split() numBeads = [] for ind in range(len(lsplit) - 2, -1, -1): if lsplit[ind].endswith('='): numBeads.append(convertToInteger(lsplit[ind + 1], 'point number in last line')) if len(numBeads) < 3: exitError('Unable to find enough point numbers in last line to assess ' +\ 'ratio between the surfaces') major = max(numBeads[0], numBeads[1]) minor = min(numBeads[0], numBeads[1]) if major > maxMajorMinor * minor + 1: prnstr(' ') prnstr('The ratio between surfaces is too high; running again with ' +\ 'lower target') prnstr(' ') newTarget = int(2 * maxMajorMinor * minor) + 1 pickcom[-1] = fmtstr('TargetNumberOfBeads {}', newTarget) pickcom.append('LimitMajorityToTarget') runcmd('pickbestseed -StandardInput', pickcom, 'stdout') break else: # ELSE ON FOR exitError('Unable to find line starting with Final: for final counts') except ImodpyError: cleanup(cleanExts, pid) prnstr(prefix + 'Running pickbestseed on tracked models') exitFromImodError(progname) cleanup(cleanExts, pid)