Fourier Correlation and Filtering

The functions listed here are defined in filtxcorr.c and scaledsobel.c. They are declared in cfsemshare.h, which includes mrcslice.h and is included by b3dutil.h.

Functions for Fourier Correlation and Filtering


int niceFrame(int num, int idnum, int limit)
int niceframe(int *num, int *idnum, int *limit)
void XCorrSetCTF(float sigma1, float sigma2, float radius1, float radius2, float *ctf, int nx, int ny, float *delta)
void setctfwsr(float *sigma1, float *sigma2, float *radius1, float *radius2, float *ctf, int *nx, int *ny, float *delta)
void XCorrSetCTFnoScl(float sigma1, float sigma2, float radius1, float radius2, float *ctf, int nx,int ny, float *delta, int *nsizeOut)
void setctfnoscl(float *sigma1, float *sigma2, float *radius1, float *radius2, float *ctf, int *nx, int *ny, float *delta, int *nsize)
void XCorrFilterPart(float *fft, float *array, int nx, int ny, float *ctf, float delta)
void filterpart(float *fft, float *array, int *nx, int *ny, float *ctf, float *delta)
void doseWeightFilter(float startDose, float endDose, float pixelSize, float afac, float bfac, float cfac, float doseScale, float *ctf, int numVals, float maxFreq, float *delta)
void doseweightfilter(float *startDose, float *endDose, float *pixelSize, float *afac, float *bfac, float *cfac, float *doseScale, float *ctf, int *numVals, float *maxFreq, float *delta)
void doseFilterValue(float startDose, float endDose, float frequency, float afac, float bfac, float cfac, float doseScale, float *atten)
void XCorrMeanZero(float *array, int nxdim, int nx, int ny)
void meanzero(float *array, int *nxdim, int *nx, int *ny)
void XCorrPeakFindWidth(float *array, int nxdim, int ny, float *xpeak, float *ypeak, float *peak, float *width, float *widthMin, int maxPeaks, float minStrength)
void xcorrpeakfindwidth(float *array, int *nxdim, int *ny, float *xpeak, float *ypeak, float *peak, float *width, float *widthMin, int *maxPeaks, float *minStrength)
void XCorrPeakFind(float *array, int nxdim, int ny, float *xpeak, float *ypeak, float *peak, int maxPeaks)
void xcorrpeakfind(float *array, int *nxdim, int *ny, float *xpeak, float *ypeak, float *peak, int *maxPeaks)
void setPeakFindLimits(int limXlo, int limXhi, int limYlo, int limYhi, int useEllipse)
void setpeakfindlimits(int *limXlo, int *limXhi, int *limYlo, int *limYhi, int *useEllipse)
double parabolicFitPosition(float y1, float y2, float y3)
double parabolicfitposition(float *y1, float *y2, float *y3)
void conjugateProduct(float *array, float *brray, int nx, int ny)
void conjugateproduct(float *array, float *brray, int *nx, int *ny)
double CCCoefficientTwoPads(float *array, float *brray, int nxdim, int nx, int ny, float xpeak, float ypeak, int nxpadA, int nypadA, int nxpadB, int nypadB, int minPixels, int *nsum)
double cccoefficienttwopads(float *array, float *brray, int *nxdim, int *nx, int *ny, float *xpeak, float *ypeak, int *nxpadA, int *nypadA, int *nxpadB, int *nypadB, int *minPixels, int *nsum)
double XCorrCCCoefficient(float *array, float *brray, int nxdim, int nx, int ny, float xpeak, float ypeak, int nxpad, int nypad, int *nsum)
double cccoefficient(float *array, float *brray, int *nxdim, int *nx, int *ny, float *xpeak, float *ypeak, int *nxpad, int *nypad, int *nsum)
void sliceGaussianKernel(float *mat, int dim, float sigma)
void slicegaussiankernel(float *mat, int *dim, float *sigma)
void scaledGaussianKernel(float *mat, int *dim, int limit, float sigma)
void scaledgaussiankernel(float *mat, int *dim, int *limit, float *sigma)
void applyKernelFilter(float *array, float *brray, int nxdim, int nx, int ny, float *mat, int kdim)
void applykernelfilter(float *array, float *brray, int *nxdim, int *nx, int *ny, float *mat, int *kdim)
void wrapFFTslice(float *array, float *tmpArray, int nx, int ny, int direction)
void wrapfftslice(float *array, float *tmpArray, int *nx, int *ny, int *direction)
int indicesForFFTwrap(int ny, int direction, int *iyOut, int *iyLow, int *iyHigh)
void fourierShiftImage(float *fft, int nxPad, int nyPad, float dx, float dy, float *temp)
void fouriershiftimage(float *fft, int *nxPad, int *nyPad, float *dx, float *dy, float *temp)
void fourierReduceImage(float *fftIn, int nxrIn, int nyrIn, float *fftOut, int nxrOut, int nyrOut, float dxIn, float dyIn, float *temp)
void fourierreduceimage(float *fftIn, int *nxrIn, int *nyrIn, float *fftOut, int *nxrOut, int *nyrOut, float *dxIn, float *dyIn, float *temp)
void fourierRingCorr(float *ffta, float *fftb, int nxReal, int ny, float *ringCorrs, int maxRings, float deltaR, float *temp)
void fourierringcorr(float *ffta, float *fftb, int *nxReal, int *ny, float *ringCorrs, int *maxRings, float *deltaR, float *temp)
int fourierCropSizes(int size, float factor, float padFrac, int minPad, int niceLimit, int *fullPadSize, int *cropPadSize, float *actualFac)

Sobel Filter Function


int scaledSobel(float *inImage, int nxin, int nyin, float scaleFac, float minInterp, int linear, float center, float *outImage, int *nxout, int *nyout, float *xOffset, float *yOffset)

Functions for Fourier Correlation and Filtering

int niceFrame(int num, int idnum, int limit)

Returns num if it is even and has no prime factor greater than limit, or makes the number even and adds idnum until it reaches a value with this property.  Use a values of 2 for idnum and call niceFFTlimit to obtain an optimal value of limit for taking the FFT with the current IMOD package. This value is 15 when linked with FFTW (because higher values use slower algorithms) and 19 with the old IMOD FFT routines (an absolute limit in that case).

int niceframe(int *num, int *idnum, int *limit)

Fortran wrapper for niceFrame

void XCorrSetCTF(float sigma1, float sigma2, float radius1, float radius2, float *ctf, int nx, int ny, float *delta)

Takes the filter parameters sigma1, sigma2, radius1, and radius2 and sets up a contrast transfer (filter) function in the array ctf, which should be dimensioned to at least 8193.  The dimensions of the real space image being filtered are specified in nx and ny, and the step size between the values in ctf is returned in delta.  If no filtering is selected (sigma1 and sigma2 both 0) then 0 is returned in delta. The values in ctf are scaled so that mean of all but the zero element is 1.

void setctfwsr(float *sigma1, float *sigma2, float *radius1, float *radius2, float *ctf, int *nx, int *ny, float *delta)

Fortran wrapper for XCorrSetCTF

void XCorrSetCTFnoScl(float sigma1, float sigma2, float radius1, float radius2, float *ctf, int nx,int ny, float *delta, int *nsizeOut)

Sets up a filter function in ctf as described for XCorrSetCTF, but with no scaling, such that the filter value will be 1 for all unattenuated frequencies.  The number of elements computed in ctf is returned in nizeOut.

void setctfnoscl(float *sigma1, float *sigma2, float *radius1, float *radius2, float *ctf, int *nx, int *ny, float *delta, int *nsize)

Fortran wrapper for XCorrSetCTFnoScl

void XCorrFilterPart(float *fft, float *array, int nx, int ny, float *ctf, float delta)

Applies the filter in ctf to the 2D Fourier transform in fft, and puts result into array, which can be the same as fft.  The dimensions of the real space image are given by nx and ny, and the dimension of the real image array are assumed to be nx + 2 by ny.  delta is the interval in reciprocal space (1/pixel) for the function in ctf.

void filterpart(float *fft, float *array, int *nx, int *ny, float *ctf, float *delta)

Fortran wrapper to XCorrFilterPart

void doseWeightFilter(float startDose, float endDose, float pixelSize, float afac, float bfac, float cfac, float doseScale, float *ctf, int numVals, float maxFreq, float *delta)

Computes a dose weighting filter.  Arguments are:
  startDose, endDose  -  Dose at start and end of image acquisition, in electrons per square Angstrom
  pixelSize  -  Pixel size in Angstroms
  afac, bfac, cfac  -  Three factors in critical dose equation
      afac * frequency
bfac + cfac
  pass 0 to use the default for any one of these
  doseScale  -  Amount by which to scale the critical and optimal doses; it should include the factor of 0.8 for 200 instead of 300 kV
  ctf  -  Array to fill with filter
  numVals  -  Number of values to compute in the array
  maxFreq  -  Frequency of last value placed in the array in reciprocal pixels
  delta  -  Returned with spacing in reciprocal pixels between array values

void doseweightfilter(float *startDose, float *endDose, float *pixelSize, float *afac, float *bfac, float *cfac, float *doseScale, float *ctf, int *numVals, float *maxFreq, float *delta)

Fortran wrapper to doseWeightFilter

void doseFilterValue(float startDose, float endDose, float frequency, float afac, float bfac, float cfac, float doseScale, float *atten)

Returns the dose weighting filter value for frequency in reciprocal Angstroms in atten; other arguments are as in doseWeightFilter

void XCorrMeanZero(float *array, int nxdim, int nx, int ny)

Shifts the contents of array to have a zero mean.  nxdim specifies the X dimension of array, and nx and ny are the size of the data in array.

void meanzero(float *array, int *nxdim, int *nx, int *ny)

Fortran wrapper to XCorrMeanZero

void XCorrPeakFindWidth(float *array, int nxdim, int ny, float *xpeak, float *ypeak, float *peak, float *width, float *widthMin, int maxPeaks, float minStrength)

Finds the coordinates of up to the maxPeaks highest peaks in array, which is dimensioned to nxdim by ny, and returns the positions in xpeak, ypeak, and the peak values in peak.  If minStrength is greater than 0, then only those peaks that are greater than that fraction of the highest peak will be returned.  In addition, if width and widthMin are not NULL, the distance from the peak to the position at half of the peak height is measured in 8 directions, the overall mean width of the peak is returned in width, and the minimum width along one of the four axes is returned in widthMin.  The X size of the image is assumed to be nxdim - 2.  The sub-pixel position is determined by fitting a parabola separately in X and Y to the peak and 2 adjacent points.  Positions are numbered from zero and coordinates bigger than half the image size are shifted to be negative.  The positions are thus the amount to shift a second image in a correlation to align it to the first.  If fewer than maxPeaks peaks are found, then the remaining values in peaks will be -1.e30.

void xcorrpeakfindwidth(float *array, int *nxdim, int *ny, float *xpeak, float *ypeak, float *peak, float *width, float *widthMin, int *maxPeaks, float *minStrength)

Fortran wrapper to XCorrPeakFindWidth.  If maxPeaks is 1, then xpeak, ypeak, and peak can be single variables instead of arrays.

void XCorrPeakFind(float *array, int nxdim, int ny, float *xpeak, float *ypeak, float *peak, int maxPeaks)

Calls XCorrPeakFindWidth with width and widthMin NULL and minStrength 0.

void xcorrpeakfind(float *array, int *nxdim, int *ny, float *xpeak, float *ypeak, float *peak, int *maxPeaks)

Fortran wrapper to XCorrPeakFind.  If maxPeaks is 1, then xpeak, ypeak, and peak can be single variables instead of arrays.

void setPeakFindLimits(int limXlo, int limXhi, int limYlo, int limYhi, int useEllipse)

Sets limits for shifts in X and in Y for peaks found in the next call to XCorrPeakFindWidth.  The peak position must be between limXlo and limXhi in X and between limYlo and limYhi in Y.  Also, if useEllipse is non-zero, the peaks are constrained to be in the ellipse bounded by these limits (i.e., a circle if the range of limts is the same in X and Y).

void setpeakfindlimits(int *limXlo, int *limXhi, int *limYlo, int *limYhi, int *useEllipse)

Fortran wrapper to setPeakFindLimits.

double parabolicFitPosition(float y1, float y2, float y3)

Given values at three successive positions, y1, y2, and y3, where y2 is the peak value, this fits a parabola to the values and returns the offset of the peak from the center position, a number between -0.5 and 0.5.

double parabolicfitposition(float *y1, float *y2, float *y3)

Fortran wrapper to parabolicFitPosition

void conjugateProduct(float *array, float *brray, int nx, int ny)

Forms the product of complex numbers in array and the complex conjugate of numbers in brray and puts the result back in array, where the arrays are Fourier transforms of images with dimensions nx by ny.

void conjugateproduct(float *array, float *brray, int *nx, int *ny)

Fortran wrapper to conjugateProduct

double CCCoefficientTwoPads(float *array, float *brray, int nxdim, int nx, int ny, float xpeak, float ypeak, int nxpadA, int nypadA, int nxpadB, int nypadB, int minPixels, int *nsum)

Returns the cross-correlation coefficient between the images in array and brray at a shift between the images given by xpeak, ypeak.  The images have sizes nx by ny and the arrays have X dimension nxdim.  Areas on each side of nxpadA in X and nypadA in Y in array, and of nxpadB and nypadB in brray, are excluded from the correlation.  Otherwise, the correlation will include all the pixels in the overlap between the images at the given shift.  The number of pixels will be returned in nsum, but no computation will be done and 0. will be returned if this number is less than minPixels.

double cccoefficienttwopads(float *array, float *brray, int *nxdim, int *nx, int *ny, float *xpeak, float *ypeak, int *nxpadA, int *nypadA, int *nxpadB, int *nypadB, int *minPixels, int *nsum)

Fortran wrapper to CCCoefficientTwoPads

double XCorrCCCoefficient(float *array, float *brray, int nxdim, int nx, int ny, float xpeak, float ypeak, int nxpad, int nypad, int *nsum)

Calls CCCoefficientTwoPads with both excluded areas specified by nxpad, nypad and with a value of 25 for minPixels.

double cccoefficient(float *array, float *brray, int *nxdim, int *nx, int *ny, float *xpeak, float *ypeak, int *nxpad, int *nypad, int *nsum)

Fortran wrapper to XCorrCCCoefficient

void sliceGaussianKernel(float *mat, int dim, float sigma)

Fills mat with the coefficients of a dim x dim Gaussian kernel with standard deviation sigma.  The coefficients are scaled to sum to 1.

void slicegaussiankernel(float *mat, int *dim, float *sigma)

Fortran wrapper for sliceGaussianKernel

void scaledGaussianKernel(float *mat, int *dim, int limit, float sigma)

Fills mat with the coefficients of a Gaussian kernel with standard deviation sigma.  The size of the kernel is set to 3 for sigma up to 1., 5 for sigma up to 2., etc., up to the size given by limit. The size is returned in dim. The coefficients are scaled to sum to 1.  Values are placed sequentially in mat; i.e., it is not treated as a 2D array of size limit.

void scaledgaussiankernel(float *mat, int *dim, int *limit, float *sigma)

Fortran wrapper for scaledGaussianKernel

void applyKernelFilter(float *array, float *brray, int nxdim, int nx, int ny, float *mat, int kdim)

Applies the kernel in mat, with the size given in kdim, to a whole floating point image in array with X and Y sizes nx amd ny and X dimension nxdim, and places the output into brray with the same X dimension.  Scaling coefficients should be sequential in mat with the X dimension advancing fastest. The filtering is parallelized with OpenMP with the same size-dependent limitation on the number of threads as in cubinterp

void applykernelfilter(float *array, float *brray, int *nxdim, int *nx, int *ny, float *mat, int *kdim)

Fortran wrapper for applyKernelFilter

void wrapFFTslice(float *array, float *tmpArray, int nx, int ny, int direction)

Converts a 2D FFT or Z slice of a 3D FFT  between the organization used by FFT routines, with the origin at the lower left, and the organization used in MRC files, with the origin at the middle left.  The complex data are packed in array with nx and ny complex elements in X and Y, direction is 0 to convert from FFT routine organization to file organization, and non-zero for the reverse.  tmpArray is a temporary array that must be dimensioned to at least 2 * nx.

void wrapfftslice(float *array, float *tmpArray, int *nx, int *ny, int *direction)

Fortran wrapper for wrapFFTslice

int indicesForFFTwrap(int ny, int direction, int *iyOut, int *iyLow, int *iyHigh)

Produces indices needed for converting between FFT routine organization and MRC file organization.  The size of the dimension is in ny and the direction is 0 for FFT routine to file, or non-zero otherwise.  Returns the starting index to copy data from on the low and high ends of the dimension in iyLow and iyHigh, respectively, and the starting index to copy data from the high end to in iyOut.  The return value is the increment for these indices.  In all cases, a loop is run ny / 2 times.  For ny even, the loop simply swaps each element or line between iyHigh to iyLow.  For ny odd, first the element at iyOut is saved, then the loop is run to copy from iyHigh to iyOut and from iyLow to iyHigh, then the saved element is copied to iyOut.

void fourierShiftImage(float *fft, int nxPad, int nyPad, float dx, float dy, float *temp)

Applies a phase shift to the Fourier transform in fft to produce a shift of dx and dy in the corresponding real-space image.  The size of the image is given in nxPad, nyPad and the X-dimension of fft is assumed to be nxPad + 2.  temp is a temporary array that must be dimensioned to at least nxPad + 2.

void fouriershiftimage(float *fft, int *nxPad, int *nyPad, float *dx, float *dy, float *temp)

Fortran wrapper for fourierShiftImage.

void fourierReduceImage(float *fftIn, int nxrIn, int nyrIn, float *fftOut, int nxrOut, int nyrOut, float dxIn, float dyIn, float *temp)

Extracts a portion of the Fourier transform in fftIn into fftOut to achieve an image reduction from real-space image size nxrIn, nyrIn to nxrOut, nyrOut, after shifting the image by dxIn, dyIn.  A Fourier shift is required to keep the image centered anyway, so it is efficient to combine this with another shift if one is needed. The X-dimension of fftIn is assumed to be nxrIn + 2 and that of fftOut is assumed to be nxrOut + 2.  temp is a temporary array that must be dimensioned to at least nxPad + 2.

void fourierreduceimage(float *fftIn, int *nxrIn, int *nyrIn, float *fftOut, int *nxrOut, int *nyrOut, float *dxIn, float *dyIn, float *temp)

Fortran wrapper for fourierReduceImage.

void fourierRingCorr(float *ffta, float *fftb, int nxReal, int ny, float *ringCorrs, int maxRings, float deltaR, float *temp)

Computes the Fourier Ring Correlation between two FFTs in ffta and fftb with real-space image size of nxReal by ny, where each line has nxReal + 2 floating point numbers.  The ring spacing is specified by deltaR and maxRings sets the maximum number of rings to compute.  The correlation coefficients are returned into ringCorrs.  temp is a temporary array that must have at least 4 * maxRings elements; it can be the same as ringCorrs.

void fourierringcorr(float *ffta, float *fftb, int *nxReal, int *ny, float *ringCorrs, int *maxRings, float *deltaR, float *temp)

Fortran wrapper for fourierRingCorr

int fourierCropSizes(int size, float factor, float padFrac, int minPad, int niceLimit, int *fullPadSize, int *cropPadSize, float *actualFac)

Computes padded sizes in one dimension for Fourier reduction, give the original size in size, the reduction factor in factor, fraction to pad in oadFrac, minimum padding in minPad, and the nice or necessary limiting factor for FFTs in niceLimit. The padded size for the original and cropped data are returned in fullPadSize and cropPadSize and the actual factor being applied is returned in actualFac. The reduction factor or its multiple by 2, 3, 4, 5, 6, 8, or 10 is required to be an integer, where the factor must be accurate to the third decimal place for a multiple to be acceptable.  The return value is 1 for a factor not greater than 1 or 2 if it is not an acceptable fractional value.

Sobel Filter Function

int scaledSobel(float *inImage, int nxin, int nyin, float scaleFac, float minInterp, int linear, float center, float *outImage, int *nxout, int *nyout, float *xOffset, float *yOffset)

Applies a Sobel-type gradient filter to an image after scaling it by a specified amount.  When an image is being scaled down, some of the scaling can be done with binning and the remainder is done with interpolation.
  inImage - input image
  nxin, nyin - size of image, which is assumed to be contiguous in the image array (X dimension equals nxin.)
  scaleFac - overall amount by which image is being scaled down (a value < 1 will scale an image up).
  minInterp - minimum amount of this scaling to be done by interpolation.
  linear - set to 1 to use linear rather than cubic interpolation, or to -1 to use zoomFiltInterp instead of binning followed by cubinterp
  center - weighting of center pixel, 2 for Sobel or 1 for Prewitt filter
  outImage - output image
  nxout, nyout - returned with size of output image
  xOffset, yOffset - coordinate offset in output image:
     input_coord = output_coord * scalingFactor + xOffset
The return value is 1 for a failure to allocate a temporary array, or the return value from zoomdown if it gives an error.  
If inImage is NULL or center < 0, the function will compute the output sizes and offsets and return.  
If center is 0, the scaled image is computed and returned without Sobel filtering.
The filtering operation is parallelized with OpenMP with the same limitation on number of threads as used in cubinterp.  
The call from Fortran is the same as that from C.