Image Interpolation and Reduction and Linear Transforms

The functions listed here are defined in cubinterp.c, zoomdown.c, reduce_by_binning.c, linearxforms.c, and amat_to_rotmagstr.c. They are declared in cfsemshare.h, which includes mrcslice.h and is included by b3dutil.h.

Function for Image Transformation by Interpolation


void cubinterp(float *array, float *bray, int nxa, int nya, int nxb, int nyb, float amat[2][2], float xc, float yc, float xt, float yt, float scale, float dmean, int linear)

Functions for Image Reduction with Filtering


int selectZoomFilter(int type, double zoom, int *outWidth)
int selectzoomfilter(int *type, float *zoom, int *outWidth)
int selectZoomFilterXY(int type, double xzoom, double yzoom, int *outWidthX, int *outWidthY)
int selectzoomfilterXY(int *type, float *xzoom, float *yzoom, int *outWidthX, int *outWidthY)
void setZoomValueScaling(float factor)
int zoomWithFilter(unsigned char **slines, int aXsize, int aYsize, float aXoff, float aYoff, int bXsize, int bYsize, int bXdim, int bXoff, int dtype, void *outData, b3dUInt32 *cindex, unsigned char *bindex)
int zoomwithfilter(float *array, int *aXsize, int *aYsize, float *aXoff, float *aYoff, int *bXsize, int *bYsize, int *bXdim, int *bXoff, float *outData)
int zoomFiltInterp(float *array, float *bray, int nxa, int nya, int nxb, int nyb, float xc, float yc, float xt, float yt, float dmean)
int zoomfiltinterp(float *array, float *bray, int *nxa, int *nya, int *nxb, int *nyb, float *xc, float *yc, float *xt, float *yt, float *dmean)
double zoomFiltValue(float radius)
double zoomfiltvalue(float *radius)

Binning Function


int extractAndBinIntoArray(void *array, int type, int nxDim, int xStart, int xEnd, int yStart, int yEnd, int nbin, void *brray, int nxBdim, int bxOffset, int byOffset, int keepByte, int *nxr, int *nyr)
int extractWithBinning(void *array, int type, int nxDim, int xStart, int xEnd, int yStart, int yEnd, int nbin, void *brray, int keepByte, int *nxr, int *nyr)
int extractwithbinning(float *array, int *nxDim, int *xStart, int *xEnd, int *yStart, int *yEnd, int *nbin, float *brray, int *keepByte, int *nxr, int *nyr)
void irepak(void *brray, void *array, int *nxin, int *nyin, int *xStart, int *xEnd, int *yStart, int *yEnd)
int reduceByBinning(void *array, int type, int nxin, int nyin, int nbin, void *brray, int keepByte, int *nxr, int *nyr)
void reduce_by_binning(float *array, int *nx, int *ny, int *nbin, float *brray, int *nxr, int *nyr)
void binIntoSlice(float *array, int nxDim, float *brray, int nxBin, int nyBin, int binFacX, int binFacY, float zWeight)
void binintoslice(float *array, int *nxDim, float *brray, int *nxBin, int *nyBin, int *binFacX, int *binFacY, float *zWeight)

Functions for Working with Linear Transformations


void xfUnit(float *f, float val, int rows)
void xfCopy(float *f1, int rows1, float *f2, int rows2)
void xfMult(float *f1, float *f2, float *prod, int rows)
void xfInvert(float *f, float *finv, int rows)
void xfApply(float *f, float xcen, float ycen, float x, float y, float *xp, float *yp, int rows)
void anglesToMatrix(float *angles, float *matrix, int rows)
void icalc_matrix(float *angles, float *matrix)
int matrixToAngles(float *matrix, double *x, double *y, double *z, int rows)
void icalc_angles(float *angles, float *matrix)

Functions for Converting Linear Transformation


void amatToRotmagstr(float a11, float a12, float a21, float a22, float *theta, float *smag, float *str, float *phi)
void amat_to_rotmagstr(float *amat, float *theta, float *smag, float *str, float *phi)
void rotmagstrToAmat(float theta, float smag, float str, float phi, float *a11, float *a12, float *a21, float *a22)
void rotmagstr_to_amat(float *theta, float *smag, float *str, float *phi, float *amat)
void amatToRotMag(float a11, float a12, float a21, float a22, float *theta, float *ydtheta, float *smag, float *ydmag)
void amat_to_rotmag(float *amat, float *theta, float *ydtheta, float *smag, float *ydmag)
void rotMagToAmat(float theta, float ydtheta, float smag, float ydmag, float *a11, float *a12, float *a21, float *a22)
void rotmag_to_amat(float *theta, float *ydtheta, float *smag, float *ydmag, float *amat)

Function for Image Transformation by Interpolation

void cubinterp(float *array, float *bray, int nxa, int nya, int nxb, int nyb, float amat[2][2], float xc, float yc, float xt, float yt, float scale, float dmean, int linear)

Applies a linear transformation to an image using cubic or linear interpolation.  It eliminates all range tests from the inner loop, but falls back from cubic to quadratic interpolation (which range tests) around the edges of the input image area.  
array - The input image array  
bray - The output image array  
nxa,nya - The dimensions of array  
nxb,nyb - The dimensions of bray  
amat - A 2x2 matrix to specify rotation, scaling, and skewing  
xc,yc - The coordinates of the center of array  
xt,yt - The translation to add to the final image. The center of the output array is taken as nxb / 2., nyb / 2.  
scale - A multiplicative scale factor for the intensities  
dmean - Mean intensity of image or other value with which to fill empty image area  
linear - Set greater than 0 to do linear interpolation or less than zero for nearest neighbor interpolation
The coordinate transformation from (Xi, Yi) in the input image to the (Xo, Yo) in the output image is given by:  
Xo = a11(Xi - Xc) + a12(Yi - Yc) + nxb/2. + xt  
Yo = a21(Xi - Xc) + a22(Yi - Yc) + nyb/2. + yt  
where Xi is a coordinate running from 0 at the left edge of the first pixel to nxa at the right edge of the last pixel in X, and similarly for Y.  
When calling from C, indices in amat are transposed: a11 = amat[0][0], a12 = amat[1][0], a21 = amat[0][1], a22 = amat[1][1]  
To call from Fortran, use arguments of the same type and in the same order. Indices in amat are in logical order: a11 = amat(1,1), a12 = amat(1,2), etc.  
This routine is now parallelized with OpenMP.  The allowed threads is proportional to the square root of output image area, falling to 1 at about 32x32.

Functions for Image Reduction with Filtering

int selectZoomFilter(int type, double zoom, int *outWidth)

Selects which filter to use in image reduction with type and sets the scaling factor with zoom, a value that must be less than 1.  The type can be 0 for a box (equivalent to binning), 1 for a Blackman window, 2 for a triangle filter, 3 for a Mitchell filter, or 4 or 5 for Lanczos 2 or Lanczos 3.  The total width of the filter in the source image is returned in outWidth.  Returns 1 for a filter type out of range or 2 for a zoom out of range.

int selectzoomfilter(int *type, float *zoom, int *outWidth)

Fortran wrapper for selectZoomFilter

int selectZoomFilterXY(int type, double xzoom, double yzoom, int *outWidthX, int *outWidthY)

Like selectZoomFilter but it allows separate zooms in X and Y given by xzoom and yzoom, with respective filter widths returned in outWidthX and outWidthY.

int selectzoomfilterXY(int *type, float *xzoom, float *yzoom, int *outWidthX, int *outWidthY)

Fortran wrapper for selectZoomFilterXY

void setZoomValueScaling(float factor)

Sets a factor for scaling the output values from image reduction of data types other than bytes and RGB.

int zoomWithFilter(unsigned char **slines, int aXsize, int aYsize, float aXoff, float aYoff, int bXsize, int bYsize, int bXdim, int bXoff, int dtype, void *outData, b3dUInt32 *cindex, unsigned char *bindex)

Reduces an image using the interpolation filter and zoom specified in selectZoomFilter.
slines - array of line pointers for the input image
aXsize, aYsize - size of input image
aXoff, aYoff - coordinate in the input image at which the lower left edge of the lower left pixel of the output starts, where input coordinates are zero at the lower left edge of the first input pixel
bXsize, bYsize - size of output image to be created, potentially from a subset of the input in X or Y and into a subset of the output array
bXdim - X dimension of the full output image array
bXoff - Index of first pixel in X to fill in output array
dtype - Type of data, a SLICE_MODE_... value.  BYTE, FLOAT, RGB, SHORT, and USHORT are allowed
outData - Output array, or address of first line to fill in a larger output array. The array is the same data type as the input unless mapping is being done, in which case it must be unsigned integers.
cindex - Unsigned integer RGBA values to map byte or short data into, or NULL for no mapping
bindex - Byte values to map RGB data into, or NULL for no mapping.  Each channel is mapped the same and a fourth channel is added with zero.
Returns 1 if no filter has been selected, 2 for an unsupported data type, 3 for an attempt to map float data, 4 if needed input coordinates to compose the output image go out of range, or 5 for a memory allocation error.
This function is parallelized with OpenMP and uses the same formula as cubinterp for reducing the number of threads for small images.

int zoomwithfilter(float *array, int *aXsize, int *aYsize, float *aXoff, float *aYoff, int *bXsize, int *bYsize, int *bXdim, int *bXoff, float *outData)

Fortran wrapper for zoomWithFilter

int zoomFiltInterp(float *array, float *bray, int nxa, int nya, int nxb, int nyb, float xc, float yc, float xt, float yt, float dmean)

Performs image reduction using zoomWithFilter and can be used the same way as cubinterp could be for scaling an image down and shifting it.  The zoom and filter must already be specified with selectZoomFilter .  
array - The input image array  
bray - The output image array  
nxa,nya - The dimensions of array  
nxb,nyb - The dimensions of bray  
xc,yc - The coordinates of the center of array  
xt,yt - The translation to add to the final image. The center of the output array is taken as nxb / 2., nyb / 2.  
dmean - Mean intensity of image or other value with which to fill empty image area  
Returns the same error values as zoomWithFilter

int zoomfiltinterp(float *array, float *bray, int *nxa, int *nya, int *nxb, int *nyb, float *xc, float *yc, float *xt, float *yt, float *dmean)

Fortran wrapper for zoomFiltInterp

double zoomFiltValue(float radius)

Returns the normalized value of the filter selected by selectZoomFilter at the distance radius from the center of the filter.  Returns 0 if no filter was selected.

double zoomfiltvalue(float *radius)

Fortran wrapper for zoomFiltValue

Binning Function

int extractAndBinIntoArray(void *array, int type, int nxDim, int xStart, int xEnd, int yStart, int yEnd, int nbin, void *brray, int nxBdim, int bxOffset, int byOffset, int keepByte, int *nxr, int *nyr)

Extracts a portion of an array into a portion of another array with optional summing of pixels (binning) by the factor nbin.  array has the input, with X dimension nxDim, and the range of coordinates to extract is given in xStart, xEnd, yStart, and yEnd * (inclusive, numbered from 0).  The slice mode is in type, and can be byte, short, unsigned short, float, or RGB. brray receives the output, with the output size returned in nxr and nyr. nxBdim can specify the X dimension of brray or be 0 to have data packed contigously (X dimension equals X size); bxOffset and byOffset can indicate X and Y offsets into brray.  The output will have the same mode as the input and will be the average of binned values, except in two cases.  If the input is bytes or RGB and keepByte is 0 or -1, then the output will be signed short integers and will be the sum or the average, respectively, of the binned values, and for RGB data the output will be an equal sum or average of red, green, and blue values.  brray can be the same as array, provided that its X dimension is not bigger than the input size, unless the binning is 1, the input is bytes and keepByte is not greater than 0. The output size is obtained by integer division of the input size by the binning.  If the remainder of this division is nonzero, the data are centered in the output array as nearly as possible.  Specifically, the coordinates of the lower left corner of the output array are offset by
  ((nx % nbin) / 2, (ny % nbin) / 2)
relative to the input array.  Returns 1 for an unsupported data type, 2 for an input X range bigger than nxDim, or 3 for an output range bigger than the output X dimension plus offset.

int extractWithBinning(void *array, int type, int nxDim, int xStart, int xEnd, int yStart, int yEnd, int nbin, void *brray, int keepByte, int *nxr, int *nyr)

Extracts a portion of an array into another array with optional summing of pixels (binning) by the factor nbin.  It calls extractAndBinIntoArray with output X dimension and offsets of 0 to produce a contiguously packed output in brray.  Other arguments are as described there.  brray can be the  same as array unless the binning is 1, the input is bytes and keepByte is not greater than 0.  Returns 1 for an unsupported data type or 2 for an input X range bigger than nxDim.

int extractwithbinning(float *array, int *nxDim, int *xStart, int *xEnd, int *yStart, int *yEnd, int *nbin, float *brray, int *keepByte, int *nxr, int *nyr)

Fortran wrapper for extractWithBinning with floating point data

void irepak(void *brray, void *array, int *nxin, int *nyin, int *xStart, int *xEnd, int *yStart, int *yEnd)

Fortran wrapper for call to extractWithBinning to simply repack a portion of a 2D array sequentially into a 1-D array.

int reduceByBinning(void *array, int type, int nxin, int nyin, int nbin, void *brray, int keepByte, int *nxr, int *nyr)

Reduces a full array in size by summing pixels (binning) by the factor nbin. array has the input, with dimensions nxin by nyin.  It simply calls extractWithBinning with ranges of 0 to nxin - 1 and 0 to nyin - 1; other arguments have the same meaning as there.  brray can be the same as array, and nxr, nyr can be the same variables as nxin, nyin.  Returns 1 for an unsupported data type.

void reduce_by_binning(float *array, int *nx, int *ny, int *nbin, float *brray, int *nxr, int *nyr)

Fortran wrapper for reduceByBinning with floating point data, called as reduce_by_binning. Again, the output variables can safely be the same as the input variables.

void binIntoSlice(float *array, int nxDim, float *brray, int nxBin, int nyBin, int binFacX, int binFacY, float zWeight)

Bins a slice of data in array, with X dimension nxDim, by binning factors of binFacX and binFacY in X and Y, and adds it into the slice in brray with a weighting of zWeight. The number of binned pixels to produce in X and Y is given by nxBin and nyBin, and brray is contiguous (has X dimension nxBin).

void binintoslice(float *array, int *nxDim, float *brray, int *nxBin, int *nyBin, int *binFacX, int *binFacY, float *zWeight)

Fortran wrapper for binIntoSlice

Functions for Working with Linear Transformations

void xfUnit(float *f, float val, int rows)

Initializes transform f with a11 and a22 set equal to val; use 1.0 for a unit transform.

void xfCopy(float *f1, int rows1, float *f2, int rows2)

Copies transform f1 to f2, which have rows1 and rows2 rows, respectively.

void xfMult(float *f1, float *f2, float *prod, int rows)

Multiples transform f1 (the one applied first) by f2 (the one applied second) and places the result in prod, which can be the same as f1 or f2.

void xfInvert(float *f, float *finv, int rows)

Takes the inverse of transform f and returns the result in finv, which can be the same as f.

void xfApply(float *f, float xcen, float ycen, float x, float y, float *xp, float *yp, int rows)

Applies transform f to the point x, y, with the center of transformation at xcen, ycen, and returns the result in xp, yp, which can be the same as x, y.

void anglesToMatrix(float *angles, float *matrix, int rows)

Given a set of rotations angles about the X, Y, and Z axes in the first, second, and third elements of angles, finds the 3-D rotation matrix and returns it in matrix. The number of rows in the matrix is specified by rows, which should be 3 when calling with a packed 3x3 matrix, or 4 when calling with an Imat data array.  The order of data elements in the array is r11, r21, r31, etc.  These are the conventions:
R premultiplies column vector of coordinates:
|xnew|       | r11  r12  r13 | |xold|  
|ynew|   =   | r21  r22  r23 | |yold|  
|znew|       | r31  r32  r33 | |zold|  

rotations applied in the order Z first, X last  
ie  R = XYZ  

rotations are right-handed - i.e., positive angle rotates counterclockwise when looking down the respective axis

X =
|   1    0     0   |  
|   0   cosa -sina |  
|   0   sina  cosa |  

Y =
|  cosb  0   sinb  |  
|   0    1     0   |  
| -sinb  0   cosb  |  

Z =
|  cosg -sinb  0   |  
|  sing  cosg  0   |  
|   0    0     1   |  

void icalc_matrix(float *angles, float *matrix)

Fortran wrapper to anglesToMatrix assuming matrix has 3 rows

int matrixToAngles(float *matrix, double *x, double *y, double *z, int rows)

Given a 3D rotation matrix in matrix, whose number of rows is rows, it finds the angles of rotation about the three axes, in the order Z, Y, X, and returns them in x, y, and z.  Returns 1 if the determinant of the matrix is not near 1. Angles are in degrees.  The conventions of anglesToMatrix are followed.

void icalc_angles(float *angles, float *matrix)

Fortran wrapper to matrixToAngles that assumes matrix has three rows and that places X, Y, and Z rotations into angles.  If the determinant is not 1, it exits with matrix and determinant output and an error message.  This function is callable from C with this name.

Functions for Converting Linear Transformation

void amatToRotmagstr(float a11, float a12, float a21, float a22, float *theta, float *smag, float *str, float *phi)

Converts a 2 by 2 transformation matrix into four "natural" parameters of image transformation.  The transformation is specified by a11, a12, a21, and a22, where
   x' = a11 * x + a12 * y
   y' = a21 * x + a22 * y   
In the converted transformation, theta is overall rotation, smag is overall magnification, str is a unidirectional stretch, and phi is the angle of the stretch axis.  Two equivalent solutions are possible, with the stretch axis in the first or fourth quadrant.  The function returns the solution that makes the magnification smag nearer to 1.0.

void amat_to_rotmagstr(float *amat, float *theta, float *smag, float *str, float *phi)

Fortran wrapper to amatToRotmagstr, where amat is dimensioned (2,*) or otherwise has elements in the order a11, a21, a12, a22.

void rotmagstrToAmat(float theta, float smag, float str, float phi, float *a11, float *a12, float *a21, float *a22)

Obtains a 2 by 2 transformation from four parameters of image transformation: theta is overall rotation, smag is overall magnification, str is a unidirectional stretch, and phi is the angle of the stretch axis, where angles are in degrees.  The transformation is returned in a11, a12, a21, and a22.

void rotmagstr_to_amat(float *theta, float *smag, float *str, float *phi, float *amat)

Fortran wrapper to rotmagstrToAmat, where amat is dimensioned (2,*) or otherwise has elements in the order a11, a21, a12, a22.

void amatToRotMag(float a11, float a12, float a21, float a22, float *theta, float *ydtheta, float *smag, float *ydmag)

Converts a 2 by 2 transformation matrix into four "semi-natural" parameters of image transformation.  The transformation is specified by a11, a12, a21, and a22, where
   x' = a11 * x + a12 * y
   y' = a21 * x + a22 * y   
In the converted transformation, theta is overall rotation, smag is overall magnification, ydtheta is the rotation of the Y axis minus the rotation of the X axis, and ydmag is the stretch along the Y axis minus the stretch along the X axis.

void amat_to_rotmag(float *amat, float *theta, float *ydtheta, float *smag, float *ydmag)

Fortran wrapper to amatToRotMag, where amat is dimensioned (2,*) or otherwise has elements in the order a11, a21, a12, a22.

void rotMagToAmat(float theta, float ydtheta, float smag, float ydmag, float *a11, float *a12, float *a21, float *a22)

Converts axial rotation and stretch specifications into a transformation matrix in a11, a12, a21, and a22.  theta is overall rotation, smag is overall magnification, ydtheta is the rotation of the Y axis minus the rotation of the X axis, and ydmag is the stretch along the Y axis minus the stretch along the X axis.

void rotmag_to_amat(float *theta, float *ydtheta, float *smag, float *ydmag, float *amat)

Fortran wrapper to rotMagToAmat, where amat is dimensioned (2,*) or otherwise has elements in the order a11, a21, a12, a22.