#!/progs/bin/python '''Subtract ice rings from set of diffraction images. For usage & options, type: "deice.py -h". Subtracts an estimate of the radially-averaged 2-theta dependent background within user-defined ice rings. (Alternatively sets to zero (unmeasured).) Can process a full batch of images (about 2 minutes each). Requires an accurate direct beam location, more accurate than usually available from data collection or integration. This program helps... Limitations: 1) Assumes that the detector is orthogonal to the beam and that the image is not distored. 2) Supports only ADSC Quantum-4 & like detectors. (Others require program modification.) 3) Ignores (passes through) data in the image header. Suggested strategy: 1) Protect the original image files against accidental over-writing. 2) Use adxv or other display program to find approx. locations of direct beam and most prominent ice ring (pixels or d-spacing). 3) Use deice options -c or -C to optimize the direct beam location for a few representative images. 4) Use options -t or -T to print statistics at single-pixel resolution for a single image. 5) Use statistics or plot thereof to determine the range of each ice ring. The limits should be in the tails of the ice diffraction so that the transition between median and interpolated values at the edge of the range is smooth. 6) Run with the -i or -I options on all images to subtract the ice rings. 7) Proceed to your favorite integration / scaling programs. ''' __author__ = "Michael S. Chapman (OHSU)" __version__ = "1.0 (Crude prototype)" __date__ = "12/09/09" import sys, optparse, math, os, time import numpy as N overflow=-999 def options(): '''Command-line option parser for program control. For usage & options, type: "deice.py -h". Full documentation "pydoc deice". ''' global __author__, __version__, __date__ opts = optparse.OptionParser( usage='%s [options] input-image-file[s]' % sys.argv[0], description="Deice diffraction images; "+" ".join( [__author__, __version__, __date__])) opts.add_option('--ice_pixels', '-I', nargs=2, action='append', type='int', help='ice ring (pixels fr center, low, high; repeatable)') opts.add_option('--ice_ring', '-i', nargs=2, action='append', type='float', help='ice ring (d-spacings, high-res, low-res; repeatable)') opts.add_option('--zero', '-z', action='store_true', help='zero ice rings (set intensity = 0)', default=0) opts.add_option('--subtract', '-s', action='store_true', help='subtract ice rings (linearly interpolate median I)', default=0) opts.add_option('--beam_pixels', '-B', nargs=2, action='store', type='float', help='beam center (x y; pixels)') opts.add_option('--beam_center', '-b', nargs=2, action='store', type='float', help='beam center (x y; mm)', default=(93.96, 94.00)) opts.add_option('--center_pixels','-C', action='store', type='int', help='ice ring (pixels fr center, low, high; for center refinement)') opts.add_option('--center_ring','-c', action='store', type='float', help='center ring (d-spacings, high-res, low-res; '+\ 'for center refinement)') opts.add_option('--detector_size', '-S', nargs=2, action='store', type='float', help='detector size (x y; pixels)', default=(2304, 2304)) opts.add_option('--pixel_size', '-p', action='store', type='float', help='pixel size (mm)', default=0.0816) opts.add_option('--detector_distance', '-D', action='store', type='float', help='crystal-detector distance (mm)', default=123.1354) opts.add_option('--wavelength', '-w', action='store', type='float', help='wavelength (Angstrom)', default=1.0) opts.add_option('--verbose', '-v', action='store_true', help='verbose output', default=0) opts.add_option('--file_qualifier', '-o', action='store', type='string', help='output file name modifier (inserted between current path & file)', default='deiced_') opts.add_option('--stats_pixels', '-T', nargs=3, action='store', type='float', default=None, help='report statistics between start, end, step (pixels, 3 floats)') opts.add_option('--stats_dspacing', '-t', nargs=3, action='store', type='float', default=(50.0, 1.5, 25), help='report statistics between low_res, high_res, '+\ '#bins (A, 3 floats, bins <=0 for 1-pixel binning)') options, arguments = opts.parse_args() # Check / convert input assert options.zero + options.subtract <= 1, \ 'Options zero (-z) and subtract (-s) are mutually exclusive.' if options.beam_pixels == None: options.beam_pixels = ( options.beam_center[0]/options.pixel_size - 0.5, options.detector_size[1] - options.beam_center[1]/options.pixel_size - 0.5) if options.verbose: print "Beam center pixels", options.beam_pixels, "from mm input." if options.ice_ring == None: options.ice_ring = [(3.977, 3.520), (2.285, 2.205), (1.945, 1.884)] if options.verbose: print "Default ice-rings (A) approx cubic ice:", ''.join( ['(%5.1f,%5.1f) '%(tuple[0],tuple[1]) for tuple in options.ice_ring]) if options.ice_pixels == None: options.ice_pixels = [] for ring in options.ice_ring: options.ice_pixels.append(( (options.detector_distance / options.pixel_size) * math.tan(2 * math.asin(options.wavelength / (2 * ring[0]))), (options.detector_distance / options.pixel_size) * math.tan(2 * math.asin(options.wavelength / (2 * ring[1]))))) if options.verbose: print 'Ice ring (pixels) fr d-spacings:', ''.join( ['(%5.1f,%5.1f) '%(tuple[0],tuple[1]) for tuple in options.ice_pixels]) for ring in options.ice_pixels: assert ring[0] < ring[1], \ "For ice rings, inner limit (lower resolution) must precede outer" if options.center_pixels == None and options.center_ring is not None: options.center_pixels = ((options.detector_distance/options.pixel_size) * math.tan(2*math.asin(options.wavelength/(2*options.center_ring)))) if options.verbose: print 'Beam refinement ring (pixels) fr d-spacings: %5.1f'%( options.center_pixels) if options.stats_pixels is None and options.stats_dspacing is not None: options.stats_pixels = [ (options.detector_distance / options.pixel_size) * math.tan(2 * math.asin(options.wavelength \ / (2 * options.stats_dspacing[0]))), (options.detector_distance / options.pixel_size) * math.tan(2 * math.asin(options.wavelength \ / (2 * options.stats_dspacing[1])))] if options.stats_dspacing[2] <= 0.0: options.stats_pixels.append(1.0) else: options.stats_pixels.append( (options.stats_pixels[1] - options.stats_pixels[0])\ /options.stats_dspacing[2]) if options.verbose: print "Range, step for stats (pixels) from d-spacings:", ''.join( ['%6.1f '%(element) for element in options.stats_pixels]) if options.stats_pixels != None: assert options.stats_pixels[0] < options.stats_pixels[1], \ "For statistics, inner limit (lower resolution) must precede outer" return (options, arguments) def conventions(options): '''Detector-specific conversion of control input to image units. Input: options contains all user input. See function options. Limitations: Specific for ADSC Quantum-4 & like detectors. Extendedable to others. ''' # # Programming notes: # Contains not quite all detector-specifics. In function options, conversion of # the beam coordinates from mm to pixel units depends on the direction & # sign of the x & y axes and the location of the origin. # # Detailed explanation: # As printed, numpy array in default 'C' configuration (not 'FORTRAN') is: # print array[row_start:row_end, col_start:col_end] # [row=0 [col=0, col=1, col=2] # row=1 [col=1, col=2, ... # ... or [[col=1, col=2, col=3], [col=1, col=2, ... # so array seems to be stored array(slow, fast) # if A is an m-row X n-column array, A[3,:] is 4th length-n "row". # page 23 of NumPy doc # Careful, because this is the opposite of matrix notation: # a(i, j) where i is the column #, j is the row. # Thus an m-by-n matrix has m rows of n elements; n columns of m elements. # http://en.wikipedia.org/wiki/Matrix_(mathematics) # In ADXV (diffraction images) coordinates are defined as either # (X_pixel, Y_pixel) or (x_mm, y_mm). # These are read into NumPy arrays as array [Y_pixel, X_pixel], i.e. X-fast. # Confusing! # Furthermore, # x_mm = (X_pixel + 0.5) * pixel_size # y_mm = (detector_size_pixels - (Y_pixel + 0.5)) * pixel_size # Noting: (1) y_mm goes in the opposite direction to Y_pixels # (2) mm coordinates are offset 0.5 pixels to center in pixel. # Reverse direction, # X_pixel = (x_mm / pixel_size) - 0.5 # Y_pixel = detector_size_pixels - (y_mm / pixel_size) - 0.5 # global overflow slow_size = options.detector_size[1] slow_center = options.beam_pixels[1] # note conventions imposed in "options" fast_size = options.detector_size[0] fast_center = options.beam_pixels[0] # note conventions imposed in "options" overflow=65535 # maximum intensity allowed / pixel return (slow_size, slow_center, fast_size, fast_center) def sort_distance(options): '''Sort image indices by distance from the direct beam. Input: options contains the command-line input. See function options. Returns: (For a ns by nf image:) sorted = a NumPy array, shape(ns.nf), of distances fr beam (increasing). hash_2d = corresonding array of shape(ns.nf,2) (sorted by increasing distance from the beam) containing the (slow, fast) indices of the original image. Note that needs to be run only once, not for each image. ''' if options.verbose: start = time.clock() print "Creating detector hash table... will take several seconds" ns, cs, nf, cf = conventions(options) # how the image is stored # Following uses (fast) numpy routines to fill elements of an image-sized # array with the pixel indices, then subtract the beam position, then # calculate the distance (pixel units) of each pixel from the direct beam. a = N.indices((ns,nf)) distances = N.hypot(a[0,...]-cs, a[1,...]-cf) hash_1d = distances.flatten().argsort(0) hash_2d = N.empty(shape=(ns*nf,2), dtype='ushort') for isort in range(ns * nf): (hash_2d[isort,0], hash_2d[isort,1]) = divmod(hash_1d[isort], nf) sorted = distances.copy().flatten() sorted.sort(0) if __debug__: # -O command-line option if 0: print "min", distances.min(), "@", distances.argmin(), \ "; max", distances.max(), "@", distances.argmax() print "start/end of hash_1d", hash_1d[:3], "/", hash_1d[-3:] print "start/end of sorted list", sorted[:3], "/", sorted[-3:] print "smallest, largest distances through hash_1d:", \ distances.flatten()[hash_1d[0]], distances.flatten()[hash_1d[1]], \ distances.flatten()[hash_1d[2]], "/", \ distances.flatten()[hash_1d[-3]], distances.flatten()[hash_1d[-2]],\ distances.flatten()[hash_1d[-1]] print "2D indices of smallest, largest distances:", hash_2d[0,:], \ hash_2d[2,:], hash_2d[3,:], "/", hash_2d[-3,:], hash_2d[-2,:], \ hash_2d[-1,:] print "smallest/largest distances from 2D hash:", \ distances[hash_2d[0,0],hash_2d[0,1]], \ distances[hash_2d[1,0],hash_2d[1,1]], \ distances[hash_2d[2,0],hash_2d[2,1]], "/", \ distances[hash_2d[-3,0],hash_2d[-3,1]], \ distances[hash_2d[-2,0],hash_2d[-2,1]], \ distances[hash_2d[-1,0],hash_2d[-1,1]] if options.verbose: print 'Hash table completed in %.0fs' % (time.clock()-start) return (sorted, hash_2d) def stats(data, distances, hash, start, stop, step, verbose, wavelength=None, detector_distance=None, pixel_size=None): '''Intensity statistics binned by distance from the direct beam. Input: data = 2D NumPy array (shape(ns, nf)) containing an image. distances = sorted (increasing) 1D NumPy array (shape(ns * nf) of pixel distances from the direct beam. hash = 2D NumPy array (shape(ns * nf, 2) of array element indices mapping each sorted distance back to the original image (data). (distances and hash are prepared by function sort_distance.) start, stop, step give the limits and increment of the table's bin size in pixel units. verbose = 1 if table is to be printed (else just returns list of median intensities). optional parameters wavelength (Angstrom), detector_distance (mm) and pixel_size (mm) are needed if tabled statistics are to be printed. Returns: a Python list of the median intensities for each bin. Printed output (for intensities in each bin): Mean, median, minimum, maximum, number pixels = 0, number overflow pixels, 1st derivative re: distance from beam (pixels) and 2nd derivative, together with the bin's d-spacing. ''' pixel_bin = N.arange(start, stop+step/2, step) bin = distances.searchsorted(pixel_bin) # bin will contain the index of the 1st element > pixel_bin. # Thus, the slice distances[bin[i]:bin[i+1]] will give all elements w/in bin mean_i = []; median_i = []; min_i = []; max_i = []; n_zero = [] n_limit = []; mean_pixel = [] n_pixels = N.empty(shape=(bin.size-1), dtype='int') for ibin in range(bin.size-1): n_pixels[ibin] = distances[bin[ibin]:bin[ibin+1]].size intensity = N.empty(shape=n_pixels[ibin], dtype='ushort') for pixel in range(n_pixels[ibin]): if 0 and __debug__: print "pixel", pixel, "bin[ibin]", bin[ibin], \ "hash", hash[bin[ibin],:], "data", \ data[hash[bin[ibin]+pixel][0],hash[bin[ibin]+pixel][1]] exit(0) intensity[pixel] \ = data[hash[bin[ibin]+pixel][0],hash[bin[ibin]+pixel][1]] median_i.append(N.median(intensity)) if verbose: mean_pixel.append(distances[bin[ibin]:bin[ibin+1]].mean()) mean_i.append(intensity.mean()) max_i.append(intensity.max()) n_zero.append(len((intensity==0).nonzero()[0])) n_limit.append(len((intensity>=overflow).nonzero()[0])) if verbose: assert detector_distance is not None and pixel_size is not None \ and wavelength is not None,\ 'If verbose - Optional arguments wavelength, detector_distance '\ + '& pixel size needed.' grad = N.gradient(N.array(median_i), step) curv = N.gradient(grad, step) print '%4s %11s %7s %7s %7s %5s %5s %4s %7s %8s'%\ ("Bin", "Pix fr beam", "d-space", "", "Median", "Max", "#zero", "#overfl", "d(medn)", "d2(medn)") for ibin in range(bin.size-1): d = wavelength / (2.0 * math.sin(math.atan2( mean_pixel[ibin] * pixel_size, detector_distance)/2.0)) print '%4d %8.2f %7.4f %7.1f %7.1f %5d %6d %7d %7.2f %8.4f'\ %(ibin, mean_pixel[ibin], d, mean_i[ibin], median_i[ibin], int(max_i[ibin]), int(n_zero[ibin]), int(n_limit[ibin]), grad[ibin], curv[ibin]) return median_i def interpolate(data, distances, hash, start, stop, verbose=0, approx_step=1.0): '''Subtract the effects of a single ice ring from a diffraction image. Method: Subtracts the difference between the radial median density (binned in circles around the direct beam), and the linear interpolations between intensities at the start and end of the ice ring. By default the bin size is ~1 pixel, and the correction is linearly interpolated between the bins for (sub-pixel) precision. Input: data = 2D NumPy array (shape(ns, nf)) containing an image. distances = sorted (increasing) 1D NumPy array (shape(ns * nf) of pixel distances from the direct beam. hash = 2D NumPy array (shape(ns * nf, 2) of array element indices mapping each sorted distance back to the original image (data). (distances and hash are prepared by function sort_distance.) start, stop = limits of ice ring (pixel units from direct beam). verbose = 1 for printed output. optional - approx_step (default 1 pixel) will be adjusted so that the range of the ice ring is a whole number multiple of the step. Returns: data, in-place corrected to remove the ice ring. Limitations: 1) Assumes detector is orthogonal to beam and gives undistorted image. 2) Assumes radial symmetry which will not be true if there are asymmetric shadows or scattering. 3) Leaves intensities of zero and overflows uncorrected so that integration programs continue to treat them as unmeasured. 4) Interpolation is linear over pixel units. Would be more correct in sin(theta)/lambda space, but likely insiginificant. ''' step = (stop-start) / nint((stop-start)/approx_step) # round to even bins if __debug__: print "step changed from", approx_step, "to", step median_i = stats(data,distances,hash,start-step/2.0,stop+step/2.0,step,0) # for debugging, replace "0)" with: __debug__ and 1, 1.0, 120.0, 0.0816) gradient = (float(median_i[-1]) - float(median_i[0])) / (stop - start) if __debug__ and 0: print "gradient=", gradient, "median[-1]=", median_i[-1], \ "median[0]=", median_i[0], "start=", start, "stop=", stop, \ "len=", len(median_i) correction = [] # will contain corrections at list of pixel steps for i in range(len(median_i)): correction.append(float(median_i[i]) - (float(median_i[0]) + i * step * gradient)) if __debug__ and 0: print "i=", i, "pixel=", start+i*step, "median=", \ median_i[i], "correction=", correction[i] if correction[i] < 0: print "***Warning: negative correction %7.1f @ pixel radius=%6.1f"%( correction[i], start + i*step) which_distances = distances.searchsorted(N.array([start, stop])) # which_distances will contain indices of the 1st element > (start or stop) if __debug__: toprint = 5 # number of points to be printed w/ debug statements n_print = (which_distances[1] - which_distances[0] -1.0) / (toprint-1.0) last_print = None n_overflow = 0 for n, i in enumerate(range(which_distances[0], which_distances[1])): # Interpolate the correction for the actual distance (bin, frac) = divmod(distances[i]-start, step); frac /= step try: local_correction = \ correction[nint(bin)]*(1.0-frac) + correction[nint(bin)+1]*frac except IndexError: if abs(frac) < 0.0001: print "Warning: IndexError from rounding corrected", i local_correction = correction[nint(bin)] else: print "i", i, "distance", distances[i], "which_distances:", \ which_distances, "frac", frac, "bin", bin, \ "nint", nint(bin), "len(median_i)", len(median_i) raise original = data[hash[i][0], hash[i][1]] if original >= overflow: n_overflow += 1 if original > 0 and original < overflow: data[hash[i][0], hash[i][1]] = \ max(nint(data[hash[i][0], hash[i][1]] - local_correction), 0) if __debug__ and nint(n//n_print) is not last_print and 0: last_print = nint(n//n_print) print "n=", n, "i=", i, "distance=", distances[i], "bin=", \ nint(bin), "frac=", frac, "neighboring corrections=", \ correction[nint(bin):nint(bin)+2], "old data=", original,\ "local_correction=", local_correction, \ "new data=", data[hash[i][0], hash[i][1]] if (n_overflow): print '***Warning %d overflow pixels'%(n_overflow),\ 'unchanged by background correction' print ' in ring %7.2f-%7.2f'%(start, stop) return data def nint(floating_point): return int(floating_point + 0.5) def zero(data, distances, hash, start, stop, verbose=0): '''Between pixel distances start & stop, zero the elements of data. This is equivalent to the exclude ice ring options of some integration programs. Interpolation should be better. ''' which_distances = distances.searchsorted(N.array([start, stop])) # which_distances will contain indices of the 1st element > (start or stop) for i in range(which_distances[0], which_distances[1]): data[hash[i][0], hash[i][1]] = 0 return data def center(options, imagelist): ''' Refine the direct beam pixel coordinates from an ice ring. Input arguments (all positions / distances in pixel units): options.beam_pixels - starting position; options.center_pixels - distance to ice ring; options.verbose - 1 for output; imagelist is a list of file names of diffraction images to use. (Usually a representative subset, not every image of a large data set.) Returns: (cs, cf) - new beam position (pixel units) with cs = center in the slowly varying dimension of the image array, cf = center in the fast changing dimension of the image array. Algorithm: From the current center searches for maximum intensity along many radii. For each image finds the center of mass of the ice ring's peak pixels. Detracting from an accurate determination are the Bragg spots. These are removed recursively by narrowing the search for the ice ring peak: a) Using maximal values only if they peak within the line search window b) Narrowing the window from an initial (say 20 pixels) to the median radius of peak pixels +/- (say) 3 sigma (+ threshold). c) Accepting only peaks within (say) 5 sigma of the median intensity for the ring. (Assuming Bragg spots cover < 1/2 area of image.) As currently configured, the search should converge if starting within ~10 pixels. Appears to be accurate to within 0.25 pixels. Limitations: 1) Assumes that the detector is orthogonal to the beam and that the image is not distorted. 2) Assumes that most of ice ring is visible (no large 2-theta offset). ''' __author__ = "Michael S. Chapman (OHSU)" __version__ = "1.0" __date = "12/09/09" ns, cs, nf, cf = conventions(options) # how the image is stored # Following uses (fast) numpy routines to fill elements of an image-sized # array with the pixel indices, then subtract the beam position, then # calculate the distance (pixel units) of each pixel from the direct beam. a = N.indices((ns,nf)) cslist = []; cflist = [] radius = 423.8 # current guess as to location of ice-ring radius = options.center_pixels if __debug__: print "Starting radius for beam center = ", radius rot_steps = 1028 #1/2 no. points to be used. Must be >3; should be large circumference = nint(2.0 * math.pi * radius) if rot_steps > circumference: print "Warning: Excessive rotation steps (", rotsteps, \ ") reset to", circumferance num_cyc = 20 search_width = 20 for ifile, file in enumerate(imagelist): if __debug__ and 0: print "Centering: ", file if options.verbose: start = time.clock() image = open(file, 'rb') header = N.fromfile(image, dtype='c', count=512, sep='') data = N.fromfile(image, dtype='H', count = ns*nf).reshape(ns, nf) lolimit = 0; hilimit = overflow for cycle in range(num_cyc): # Just to document how distances of points from center: (slow, fast) = (a[0,...]-cs, a[1,...]-cf) distances = N.hypot(slow, fast) slowindices = []; fastindices = []; distancelist = [] peakintensity = [] n_pt = 0; no_max = 0; n_0 = 0; n_overfl = 0; n_lo = 0; n_hi = 0 for rotation in range(rot_steps): n_pt += 2 angl = math.pi*(1.0/4.0+float(rotation)/float(rot_steps)) # Note rotation offset by 45 deg to lessen chance of whitespace rotated_unit = (math.cos(angl), math.sin(angl)) if __debug__ and 0: print "Angle:", angl, "Unit vector:", rotated_unit indices = [[],[]]; intensities = [[],[]]; which = [None,None] peakval = [None, None] reject = 0 for i, sign in enumerate(range(-1, 2, 2)): # Angles incremented 0 to pi, paired w/ -ve vector for # pi to 2pi, then both rejected if either fails criteria for delta in range(-(search_width/2+1), search_width/2+2): point = (rotated_unit[0]*(radius+delta)*sign, rotated_unit[1]*(radius+delta)*sign) index = (nint(point[0]+cs), nint(point[1]+cf)) # note unconvention slow before fast ready for array val = data[index[0], index[1]] if val <= 0: n_0 += 2; reject = 1; break elif val >= overflow: n_overfl += 2; reject = 1; break intensities[i].append(val) indices[i].append(index) if reject: break sorted_indices = N.array(intensities[i]).argsort().tolist() which[i] = sorted_indices[-1] if __debug__ and 0: print "highest=", which[i], "of", intensities[i] print "highest=", which[i], "of", indices[i] if which[i] == 0 or which[i] == (len(intensities[i]) - 1): no_max += 2; reject = 1; break peakval[i] = intensities[i][which[i]] if peakval[i] < lolimit: n_lo += 2; reject = 1; break if peakval[i] > hilimit: n_hi += 2; reject = 1; break if reject == 0: for i in range(2): # Add pair of points if both satisfied rejection criteria slowindices.append(indices[i][which[i]][0]) fastindices.append(indices[i][which[i]][1]) distancelist.append( distances[indices[i][which[i]][0], indices[i][which[i]][1]]) peakintensity.append(peakval[i]) n_used = n_pt - no_max - n_0 - n_overfl - n_lo - n_hi if n_used <= 0: print "*** Centering - no acceptable points found:" new_cs = cs; new_cf = cf; avdist = None; stdist = None else: new_cs = N.mean(N.array(slowindices)) new_cf = N.mean(N.array(fastindices)) medist = N.median(N.array(distancelist)) avdist = N.mean(N.array(distancelist)) stdist = N.std(N.array(distancelist)) medint = N.median(N.array(peakintensity)) avint = N.mean(N.array(peakintensity)) stdint = N.std(N.array(peakintensity)) if __debug__ and 0: print "slow mean=", new_cs, "; list: ", slowindices print "fast mean=", new_cf, "; list: ", fastindices print "distances: ", distancelist shift_cf = new_cf - cf; cf = new_cf shift_cs = new_cs - cs; cs = new_cs if __debug__ and 0: print "Cyc %d shift %5.3f-->(%7.2f,%7.2f) in %4.1fs"%( cycle, math.hypot(shift_cf,shift_cs), cf, cs, time.clock() - start) print 'Points: %d/%d used, %d/%d under/overflow,'%( n_used, rot_steps*2, n_0, n_overfl), \ '%d no peak, %d/%d too low/high.'%(no_max, n_lo, n_hi) print 'Radius mean, median: %6.1f,%6.1f+/-%3.1f;'%( avdist, medist, stdist),\ 'Peak mean, median: %5.0f,%5.0f+/-%4.0f'%( avint, medint, stdint) # Update the limits re: median & narrow search ready next cycle? radius = medist; search_width = int(min(stdist,1.001) + 1.0) * 2 + 1 lolimit = medint - 5.0 * stdint; hilimit = medint + 5.0 * stdint image.close() if options.verbose: print '%s beam shift (fast,slow) (%5.2f,%5.2f) cyc %d'%( file, shift_cf, shift_cs, cycle+1), 'to (%7.2f,%7.2f)'%(cf, cs) print 'Points: %d/%d used, %d/%d under/overflow,'%( n_used, rot_steps*2, n_0, n_overfl), \ '%d no peak, %d/%d too low/high.'%(no_max, n_lo, n_hi) print 'Ring radius: %7.2f +/-%5.2f (time:%3.0fs)'%( avdist, stdist, time.clock() - start) cslist.append(cs); cflist.append(cf) if ifile > 0: print 'Combined beam center (fast, slow) pixels for %d images:\n'%( ifile+1), \ 'Mean: (%7.2f,%7.2f); StDev: (%7.2f,%7.2f); Median: (%7.2f,%7.2f)'%( N.mean(N.array(cflist)), N.mean(N.array(cslist)), N.std(N.array(cflist)), N.std(N.array(cslist)), N.median(N.array(cflist)), N.median(N.array(cslist)) ) return (N.median(N.array(cflist)), N.median(N.array(cslist))) if __name__=='__main__': # print "Command: ", ' '.join(sys.argv) options, arguments = options() if __debug__: assert 1 < 2, "Debug termination after options read" ns, cs, nf, cf = conventions(options) # how the image is stored # QQQ may not be needed as call in sort_distance if options.center_pixels is not None: (cs, cf) = center(options, arguments) if __debug__: assert 1 < 2, "Debug termination after center" distances, hash = sort_distance(options) # returns... # sorted list of distances from direct beam; corresponding indices of image for file in arguments: outimage = os.path.join(os.path.dirname(file), options.file_qualifier+os.path.basename(file)) if options.verbose: print "Processing", outimage, "..." image = open(file, 'rb') header = N.fromfile(image, dtype='c', count=512, sep='') data = N.fromfile(image, dtype='H', count = ns*nf).reshape(ns, nf) # file size / header corresp to 512-byte header then 2304 x 2304 # unsigned 2-byte integers. # data types - see numpy book page 22, table 2.1 if options.verbose: print "Input image:" if options.stats_pixels is not None: median_i = \ stats(data, distances, hash, options.stats_pixels[0], \ options.stats_pixels[1], options.stats_pixels[2], options.verbose, options.wavelength, options.detector_distance, options.pixel_size) if __debug__: assert 1 < 2, "Debug termination after stats" for ring in options.ice_pixels: if __debug__: print "ice ring", ring if options.zero: data = zero(data, distances, hash, ring[0], ring[1], options.verbose) if options.subtract: data = interpolate(data, distances, hash, ring[0], ring[1], 1) if options.verbose: print "Output de-iced image:" if options.stats_pixels is not None: median_i = \ stats(data, distances, hash, options.stats_pixels[0], \ options.stats_pixels[1], options.stats_pixels[2], options.verbose, options.wavelength, options.detector_distance, options.pixel_size) output = open(outimage, 'wb') header.tofile(output, sep='') data.tofile(output, sep='') # Not able to get fwrite to work with extended datatypes, but tofile OK. image.close(); output.close() print "bye"