Fourier Correlation, Filtering, and Image Transformation

Functions for Fourier Correlation, Filtering, and Image Transformation

Sobel Filter Function

Functions for Fourier Correlation, Filtering, and Image Transformation

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 dose; 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.

int findManyXCorrPeaks(float *array, int nxdim, int ny, int ixOffset, int iyOffset, float *xpeak, float *ypeak, float *peak, int maxPeaks, int maxGrow, int *numFound)

Finds multiple cross-correlation peaks, such as from correlating with a template, much more efficiently than XcorrPeakFindWidth.  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, the peak values in peak, and the number of peaks in numFound.  For a template correlation where a template has been wrapped to give correlation peaks at the same location as the feature, the offsets ixOffset and iyOffset can be set to the amount of padding on each side in X and Y to ignore the borders where the array was padded.  For a regular cross-correlation to find shifts between images, pass iyOffset as -1 and ixOffset as either -1 or -2. This will allow the routine to find peaks along the borders of the correlation, where ixOffset of -1 would make it loop in X all the way out to nxdim - 1 and wrap X from 0 to nxdim - 1 when testing peaks, whereas -2 handles the more common case of a padded inverse correlation and makes it loop X to nxdim - 3 and wrap X between 0 and nxdim - 3.  In this case, the routine returns positive and negative shifts as XcorrPeakFindWidth does.  
Rather than maintaining peaks in order while finding them, up to maxGrow peaks are collected into arrays allocated for this purpose, then sorted, and the top maxPeaks peaks are retained.  The array is allowed to grow repeatedly.  After the first time, only the newest peaks need to be sorted, but all the data need to be repacked to integrate new high ones.  Values for maxGrow of 1.2 - 1.5 were found to be most efficient in a case where several thousand peaks were being selected out of several hundred thousand.  Returns 1 if maxGrow is not more than 1.05 times maxPeaks and -1 for memory error.

int findSpacedXCorrPeaks(float *array, int nxdim, int ixMin, int ixMax, int iyMin, int iyMax, float *xpeak, float *ypeak, float *peak, int maxPeaks, float minSpacing, int *numPeaks, float minStrength)

Finds multiple cross-correlation peaks with a minimum spacing between them even more efficiently than findManyXCorrPeaks.  Finds the coordinates of up to the maxPeaks highest peaks separated by minSpacing in array, whose X dimension is nxdim, within the * coordinates from ixMin to ixMax and iyMin to iyMax, inclusive.  Returns the positions in xpeak, ypeak, the peak values in peak, and the number of peaks in numPeaks.  If minStrength is non-zero, only those peaks that are greater than that fraction of the highest peak will be returned. The function divides the image into blocks of minSpacing / sqrt(2.) pixels, searches for the highest peak within each block, and stores that peak for each block.  The loop on blocks is paralellized with OpenMP. At the end, data are merged from the different threads, sorted, and weaker peaks too close to other ones eliminated.  The coordinates of the peaks within the image are returned, not offsets from the origin as in XcorrPeakFindWidth. Returns -1 for memory error.

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

double subareaCCCoefficient(float *array, float *brray, int nxdim, int xstart, int xend, int ystart, int yend, int delx, int dely)

Computes a cross-correlation coefficient in a subarea at one shift between two images in array and brray with X dimension nxdim.  The amount that the image in brray is shifted relative to array is in delx, dely, and the correlation is done between xstart and xend in X and ystart and yend in Y (inclusive coordinates in array, numbered from 0.  There are no checks for the validity of these limits.

double weightedCCCoefficient(float *array, float *brray, int nxDim, int ix0, int ix1, int iy0, int iy1, int delX, int delY, float *aWeights, float *bWeights, int nxWgt, int binning, int wgtXoffset, int wgtYoffset)

Returns a weighted cross-correlation between images in array and brray, with X X dimension nxDim, over the range of positions in array given by ix0 to ix1 in X and iy0 to iy1 in Y (inclusive, numbered from 0), at a shift of B relative to A given by delX, delY.  Weights are supplied in aWeights for A and bWeights for B, with X dimension nxWgt, binned by binning and with offsets wgtXoffset, wgtYoffset, where for position ix, iy in an image, the index of the weight is ix / binning + wgtXoffset + nxWgt * (iy / binning + wgtYoffset).  The weight at each pixel is the product of the weights for A and B.

double weightedCorrFromSums(double aSum, double aSumSq, double bSum, double bSumSq, double abSum, double wSum, double *sumArray, const char *descrip)

Returns the weighted cross-correlation from the set of weighted sums: aSum, bSum are the sums of weighted values in A and B, aSumSq and bSumSq are sums of weighted squares, abSum is the sum of weighted cross-products, and sSum is the sum of weights.  If sumArray is non-null, it will be filled with the 6 sum arguments in order.  If descrip is non-zero length,values will be printed after this text; it can be NULL.

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 nxrOut + 2, unless temp is NULL, in which case the FFT is simply cropped, fourierShiftImage is not called and dxIn and dyIn are ignored.

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 fourierExpandImage(float *fftIn, int nxrIn, int nyrIn, float *fftOut, int nxrOut, int nyrOut, float dxIn, float dyIn, float *temp)

Pads the Fourier transform in fftIn into fftOut to achieve an image expansion 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 nxrIn + 2, unless temp is NULL, in which case the FFT is simply padded, fourierShiftImage is not called and dxIn and dyIn are ignored.

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

Fortran wrapper for fourierExpandImage.

void fourierShiftVolume(float *fft, int nxPad, int nyPad, int nzPad, float dx, float dy, float dz, float *temp)

Applies a phase shift to the 3D Fourier transform in fft to produce a shift of dx, dy, and dz in the corresponding real-space volume.  The size of the volume is given in nxPad, nyPad, nzPad 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 fouriershiftvolume(float *fft, int *nxPad, int *nyPad, int *nzPad, float *dx, float *dy, float *dz, float *temp)

Fortran wrapper for fourierShiftVolume.

void fourierReduceVolume(float *fftIn, int nxrIn, int nyrIn, int nzrIn, float *fftOut, int nxrOut, int nyrOut, int nzrOut, float dxIn, float dyIn, float dzIn, float *temp)

Extracts a portion of the 3-D Fourier transform in fftIn into fftOut to achieve an volume reduction from real-space columesize nxrIn, nyrIn, nzrIn to nxrOut, nyrOut, nzrOut after shifting the volume by dxIn, dyIn, dzIn.  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 nxrOut + 2, unless temp is NULL, in which case the FFT is simply cropped, fourierShiftVolume is not called and dxIn, dyIn, and dzIn are ignored.

void fourierreducevolume(float *fftIn, int *nxrIn, int *nyrIn, int *nzrIn, float *fftOut, int *nxrOut, int *nyrOut, int *nzrOut, float *dxIn, float *dyIn, float *dzIn, float *temp)

Fortran wrapper for fourierReduceVolume.

void fourierExpandVolume(float *fftIn, int nxrIn, int nyrIn, int nzrIn, float *fftOut, int nxrOut, int nyrOut, int nzrOut, float dxIn, float dyIn, float dzIn, float *temp)

Pads the 3-D Fourier transform in fftIn into fftOut to achieve an image expansion from real-space image size nxrIn, nyrIn, nzrIn to nxrOut, nyrOut, nzrOut after shifting the image by dxIn, dyIn, dzIn.  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 nxrIn + 2, unless temp is NULL, in which case the FFT is simply padded, fourierShiftVolume is not called and dxIn, dyIn, and dzIn are ignored.

void fourierexpandvolume(float *fftIn, int *nxrIn, int *nyrIn, int *nzrIn, float *fftOut, int *nxrOut, int *nyrOut, int *nzrOut, float *dxIn, float *dyIn, float *dzIn, float *temp)

Fortran wrapper for fourierExpandImage.

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 or expansion, given the original size in size, the reduction factor or inverse of expansion factor in factor, fraction to pad in padFrac, minimum padding in minPad, and the nice or necessary limiting factor for FFTs in niceLimit. The padded size for the original and processed data are returned in fullPadSize and cropPadSize and the actual factor being applied is returned in actualFac. The reduction or expansion 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 equal to 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.