Utilities in libcfshr

The utility functions listed here are defined in various separate files. They are declared in cfsemshare.h, which includes mrcslice.h and is included by b3dutil.h.

Function to Test for Point Inside Boundary


int InsideContour(float *ptX, float *ptY, int np, float x, float y)
int insidecontour(float *ptX, float *ptY, int *np, float *x, float *y)

Function for Finding Convex Boundary


void convexBound(float *sx, float *syin, int npnts, float fracomit, float pad, float *bx, float *by, int *nvert, float *xcen, float *ycen, int maxverts)
void convexbound(float *sx, float *syin, int *npnts, float *fracomit, float *pad, float *bx, float *by, int *nvert, float *xcen, float *ycen, int *maxverts)

Function for Sorting Points onto Two Surfaces


int setSurfSortParam(int index, float value)
int setsurfsortparam(int *index, float *value)
int surfaceSort(float *xyz, int numPts, int markersInGroup, int *group)
int surfacesort(float *xyz, int *numPts, int *markersInGroup, int *group)
Algorithm and Parameters for surfaceSort

Function for Parsing Lists


int *parselist(const char *line, int *nlist)
int parselistfw(const char *line, int *list, int *nlist, int *limlist, int linelen)

Histogram Functions


void kernelHistogram(float *values, int numVals, float *bins, int numBins, float firstVal, float lastVal, float h, int verbose)
void kernelhistogram(float *values, int *numVals, float *bins, int *numBins, float *firstVal, float *lastVal, float *h, int *verbose)
int scanHistogram(float *bins, int numBins, float firstVal, float lastVal, float scanBot, float scanTop, int findPeaks, float *dip, float *peakBelow, float *peakAbove)
int scanhistogram(float *bins, int *numBins, float *firstVal, float *lastVal, float *scanBot, float *scanTop, int *findPeaks, float *dip, float *peakBelow, float *peakAbove)
int findHistogramDip(float *values, int numVals, int minGuess, float *bins, int numBins, float firstVal, float lastVal, float *histDip, float *peakBelow, float *peakAbove, int verbose)
int findhistogramdip(float *values, int *numVals, int *minGuess, float *bins, int *numBins, float *firstVal, float *lastVal, float *histDip, float *peakBelow, float *peakAbove, int *verbose)

Functions for Analyzing Gold Beads


void makeModelBead(int boxSize, float beadSize, float *array)
void makemodelbead(int *boxSize, float *beadSize, float *array)
double beadIntegral(float *array, int nxdim, int nx, int ny, float rCenter, float rInner, float rOuter, float xcen, float ycen, float *cenmean, float *annmean, float *temp, float annPct, float *median)

Functions for Image Statistics by Sampling


int sampleMeanSD(unsigned char **image, int type, int nx, int ny, float sample, int ixStart, int iyStart, int nxUse, int nyUse, float *mean, float *sd)
int samplemeansd(float *image, int *nx, int *ny, float *sample, int *ixStart, int *iyStart, int *nxUse, int *nyUse, float *mean, float *sd)
int percentileStretch(unsigned char **image, int type, int nx, int ny, float sample, int ixStart, int iyStart, int nxUse, int nyUse, float pctLo, float pctHi, float *scaleLo, float *scaleHi)

Functions for Measuring Volume Mean/SD at Multiple Binnings


int multiBinSetup(int binning[][3], int boxSize[][3], int boxSpacing[][3], int numScales, int startCoord[3], int endCoord[3], int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds)
int multibinsetup(int binning[][3], int boxSize[][3], int boxSpacing[][3], int *numScales, int *startCoord, int *endCoord, int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds)
int multiBinStats(int binning[][3], int boxSize[][3], int boxSpacing[][3], int numScales, int startCoord[3], int endCoord[3], int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds, float *buffer, float *means, float *SDs, int *funcData, int (*getSliceFunc)(int *, int *, float *))
int multibinstats(int binning[][3], int boxSize[][3], int boxSpacing[][3], int *numScales, int *startCoord, int *endCoord, int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds, float *buffer, float *means, float *SDs, int *funcData, int (*getSliceFunc)(int *, int *, float *))

Functions for Working with Montage Piece Lists


int checkPieceList(int *pclist, int stride, int npclist, int redfac, int nframe, int *minpiece, int *npieces, int *noverlap)
void checklist(int *pclist, int *npclist, int *redfac, int *nframe, int *minpiece, int *npieces, int *noverlap)
void adjustPieceOverlap(int *pclist, int stride, int npclist, int nframe, int minpiece, int noverlap, int newOverlap)
void adjustpieceoverlap(int *pclist, int *npclist, int *nframe, int *minpiece, int *noverlap, int *newOverlap)

Functions for Finding Alignments of Montage Pieces


void montXCBasicSizes(int ixy, int nbin, int indentXC, int *nxyPiece, int *nxyOverlap, float aspectMax, float extraWidth, float padFrac, int niceLimit, int *indentUse, int *nxyBox, int *numExtra, int *nxPad, int *nyPad, int *maxLongShift)
void montxcbasicsizes(int *ixy, int *nbin, int *indentXC, int *nxyPiece, int *nxyOverlap, float *aspectMax, float *extraWidth, float *padFrac, int *niceLimit, int *indentUse, int *nxyBox, int *numExtra, int *nxPad, int *nyPad, int *maxLongShift)
void montXCIndsAndCTF(int ixy, int *nxyPiece, int *nxyOverlap, int *nxyBox, int nbin, int indentUse, int *numExtra, int nxPad, int nyPad, int numSmooth, float sigma1, float sigma2, float radius1, float radius2, int evalCCC, int *ind0Lower, int *ind1Lower, int *ind0Upper, int *ind1Upper, int *nxSmooth, int *nySmooth, float *ctf, float *delta)
void montxcindsandctf(int *ixy, int *nxyPiece, int *nxyOverlap, int *nxyBox, int *nbin, int *indentUse, int *numExtra, int *nxPad, int *nyPad, int *numSmooth, float *sigma1, float *sigma2, float *radius1, float *radius2, int *evalCCC, int *ind0Lower, int *ind1Lower, int *ind0Upper, int *ind1Upper, int *nxSmooth, int *nySmooth, float *ctf, float *delta)
int montXCFindBinning(int maxBin, int targetSize, int indentXC, int *nxyPiece, int *nxyOverlap, float aspectMax, float extraWidth, float padFrac, int niceLimit, int *numPaddedPix, int *numBoxedPix)
int montxcfindbinning(int *maxBin, int *targetSize, int *indentXC, int *nxyPiece, int *nxyOverlap, float *aspectMax, float *extraWidth, float *padFrac, int *niceLimit, int *numPaddedPix, int *numBoxedPix)
void montXCorrEdge(float *lowerIn, float *upperIn, int *nxyBox, int *nxyPiece, int *nxyOverlap, int nxSmooth, int nySmooth, int nxPad, int nyPad, float *lowerPad, float *upperPad, float *lowerCopy, int numXcorrPeaks, int legacy, float *ctf, float delta, int *numExtra, int nbin, int ixy, int maxLongShift, int weightCCC, float *xDisplace, float *yDisplace, float *CCC, void (*twoDfft)(float *, int *, int *, int *), void (*dumpEdge)(float *, int *, int *, int *, int *, int *), char *debugStr, int debugLen, int debugLevel)
void montxcorredge(float *lowerIn, float *upperIn, int *nxyBox, int *nxyPiece, int *nxyOverlap, int *nxSmooth, int *nySmooth, int *nxPad, int *nyPad, float *lowerPad, float *upperPad, float *lowerCopy, int *numXcorrPeaks, int *legacy, float *ctf, float *delta, int *numExtra, int *nbin, int *ixy, int *maxLongShift, int *weightCCC, float *xDisplace, float *yDisplace, float *CCC, void (*twoDfft)(float *, int *, int *, int *), void (*dumpEdge)(float *, int *, int *, int *, int *, int *), int *debugLevel)
void montBigSearch(float *array, float *brray, int nx, int ny, int ixBox0, int iyBox0,int ixBox1, int iyBox1, float *dxMin, float *dyMin, float *sdMin, float *ddenMin, int numIter, int limStep)
void montbigsearch(float *array, float *brray, int *nx, int *ny, int *ixBox0, int *iyBox0,int *ixBox1, int *iyBox1, float *dxMin, float *dyMin, float *sdMin, float *ddenMin, int *numIter, int *limStep)
void montSdCalc(float *array, float *brray, int nx, int ny, int ixBox0, int iyBox0, int ixBox1, int iyBox1, float dx, float dy, float *sd, float *dden)
void montsdcalc(float *array, float *brray, int *nx, int *ny, int *ixBox0, int *iyBox0, int *ixBox1, int *iyBox1, float *dx, float *dy, float *sd, float *dden)
int findPieceShifts(int *ivarpc, int nvar, int *indvar, int *ixpclist, int *iypclist, float *dxedge, float *dyedge, int idir, int *pieceLower, int *pieceUpper, int *ifskipEdge, int edgeStep, float *dxyvar, int varStep, int *edgeLower, int *edgeUpper, int pcStep, int *work, int fort, int leaveInd, int skipCrit, float robustCrit, float critMaxMove, float critMoveDiff, int maxIter, int numAvgForTest, int intervalForTest, int *numIter, float *wErrMean, float *wErrMax)

Function for Computing a Scaled, Reduced Power Spectrum


int spectrumScaled(void *image, int type, int nx, int ny, void *spectrum, int padSize, int finalSize, int bkgdGray, float truncDiam, int filtType, void (*twoDfft)(float *, int *, int *, int *))

Color Map Functions


int *cmapStandardRamp()
int *cmapInvertedRamp()
int cmapConvertRamp(int *rampData, unsigned char table[3][256])
int cmapReadConvert(char *filename, unsigned char table[3][256])
Standard color ramps

Function to Test for Point Inside Boundary

int InsideContour(float *ptX, float *ptY, int np, float x, float y)

Returns 1 if the point x, y is inside or on the polygon (contour) whose vertices are in ptX, ptY, where np is the number of points, otherwise returns 0.

int insidecontour(float *ptX, float *ptY, int *np, float *x, float *y)

Fortran wrapper to InsideContour .  But it is more convenient to call
   logical function inside(ptX, ptY, np, x, y)

Function for Finding Convex Boundary

void convexBound(float *sx, float *syin, int npnts, float fracomit, float pad, float *bx, float *by, int *nvert, float *xcen, float *ycen, int maxverts)

Finds the smallest convex boundary around an arbitrary set of points in sx, sy.  npnts is the total number of points. If fracomit is nonzero, it omits that fraction of points from consideration completely, choosing the points that are farthest from the centroid of all the points.  (If fracomit is negative, the Y coordinates of the points will be temporarily scaled to have the same standard deviation as the X range in order to determine the farthest points; it then eliminates -fracomit points.) It returns the centroid of the points under consideration in xcen, ycen.  It returns the set of boundary points in bx, by and the number of points in nvert.  The size of bx and by is specified by maxverts.  If pad is nonzero, it makes the boundary points be that distance away from the outermost sx, sy points. It uses Graham's scan algorithm.  nvert is returned with -1 for an error allocating temporary arrays, or -2 for the boundary arrays not large enough.

void convexbound(float *sx, float *syin, int *npnts, float *fracomit, float *pad, float *bx, float *by, int *nvert, float *xcen, float *ycen, int *maxverts)

Fortran wrapper for convexBound

Function for Sorting Points onto Two Surfaces

int setSurfSortParam(int index, float value)

Sets one parameter for surfaceSort .

int setsurfsortparam(int *index, float *value)

Fortran wrapper for setSurfSortParam

int surfaceSort(float *xyz, int numPts, int markersInGroup, int *group)

Sorts a set of points in 3D onto two surfaces.  The X, Y, and Z coordinates of the numPts points are packed sequentially into the array xyz.  If markersInGroup is nonzero, then the group array should be passed in with 0 for points eligible to be included in the pairs of points that are used to initiate a cluster, or with a -1 for points that should not be used to initiate clusters.  The surface numbers, 1 for the bottom one and 2 for the top one, are returned in group.  The return value is 1 for a memory allocation error.

int surfacesort(float *xyz, int *numPts, int *markersInGroup, int *group)

Fortran wrapper for surfaceSort, where xyz is dimensioned to (3,*)

Algorithm and Parameters for surfaceSort

This routine works as follows:
First it fits a plane to all of the points and rotates the points to make that plane level.  Analysis proceeds with rotated positions.
For rapid access to neighboring points, it sets up a grid of squares (at sGridSpacing) and list of delta values to grid positions in successively wider rings around one grid square.  Then it makes lists of points within each grid square.
For each point (not marked as an outlier by values of -1 in the group array) it then searches for the neighboring point at the steepest angle with respect to it, looking at up to sMaxAngleNeigh neighbors.
Duplicate pairs are allowed to form, and are eliminated after sorting.
The points are sorted by steepest angle, and then it determines a set of pairs to use for building clusters, as controlled by parameters sSteepestRatio, sNumMinForAmax, sAngleMax, sNumMinForArelax, and sAngleRelax.
Given this set of pairs, it then makes a list of the separation in Z between each pair and analyzes for MAD-Median outliers with criterion sOutlierCrit.  These outliers are eliminated from cluster formation, as are the ones identified by the caller. Pairs are always considered in order from the steepest angle downward.  To build a cluster, it starts with the next point that hasn't been added to a cluster yet, putting it and its pair on a list of cluster points to check.  Each point on the cluster list is checked for whether it occurs in a pair at a steep angle, and if so its mate is then added to the opposite group.
For each cluster, a central position and an average delta Z are computed.
Finally, it searches in progressively bigger rings around the clusters for grid squares where there are unassigned points.  When it finds an unassigned point, it collects the nearest assigned points in the neighborhood up to sMaxNumFit or up to a distance of sMaxFitDist.  If this includes points from both surfaces and there are at least sBiplaneMinFit points, it fits a pair of planes to all points; otherwise it fits a single plane to just the points from one surface if there are at least sPlaneMinFit points for this; otherwise it simply gets a mean of Z positions on the top and bottom if possible.  One way or another this gives an estimate of the top and bottom Z position at the unassigned point, and it assigns the point based on which Z it is closer to.  In fact, this process is done in two rounds, and on the first round assignment is deferred if the disparity between the point and estimated surface Z values as a percentage of the Z separation exceeds sMaxRound1Dist.

The numbers in the comments below are index values for setSurfSortParam

/* Spacing of squares for sorting points and accessing rings of points */
static float sGridSpacing = 50.; /* 0 */
/* Maximum neighbors to evaluate for finding neighbor with steepest angle */
static int sMaxAngleNeigh = 50;  /* 1 */
/* Use pairs with angle more than this fraction of very steepest angle */
static float sSteepestRatio = 0.5f; /* 2 */
/* If there are fewer than this number of pairs, take pairs down to sAngleMax */
static int sNumMinForAmax = 10;  /* 3 */
/* Angle to go down to if sSteepestRatio doesn't give enough pairs */
static float sAngleMax = 20.;    /* 4 */
/* If fewer than this number of pairs, take pairs down to sAngleRelax */
static int sNumMinForArelax = 3; /* 5 */
/* Angle to go down to if sAngleMax doesn't give enough pairs */
static float sAngleRelax = 5.;   /* 6 */
/* Criterion for MAD-Median outlier elimination based on delta Z of a pair */
static float sOutlierCrit = 3.;  /* 7 */
/* Maximum distance to search for neighboring points in plane fits */
static float sMaxFitDist = 2048.; /* 8 */
/* Maximum number of points in plane fits */
static int sMaxNumFit = 15;      /* 9 */
/* Minimum # of points for fitting parallel planes */
static int sBiplaneMinFit = 5;   /* 10 */
/* Minimum # of points for fitting one plane */
static int sPlaneMinFit = 4;     /* 11 */
/* Maximum % distance of Z value between nearest and other plane for deferring to 2nd
   round: values > 50 disable deferring points */
static float sMaxRound1Dist = 100.; /* 12 */
/* 1 for minimal output, 2 for exhaustive output */
static int sDebugLevel = 0;      /* 13 */

Function for Parsing Lists

int *parselist (const char *line, int *nlist)

Converts a list entry in the string line into a set of integers, returns a pointer to the list of numbers, and returns the number of values in nlist.  An example of a list is 1-5,7,9,11,15-20. Numbers separated by dashes are replaced by all of the numbers in the range.  Numbers need not be in any order, and backward ranges (10-5) are handled.  Any characters besides digits are valid separators.  A / at the beginning of the string will return a NULL and a value of -1 in nlist; memory allocation errors return a NULL and -2 in nlist.  Parsing stops and the list is returned when a anything besides space, tab, comma, dash, or a digit occurs after a space or tab; otherwise other characters (as well as comma or dash right after the last number) make it return NULL and -3 in nlist.  Negative numbers can be entered provided that the minus sign immediately precedes the number.  E.g.: -3 - -1 or -3--1 will give -3,-2,-1; -3, -1,1 or -3,-1,1 will give -3,-1,1.

int parselistfw(const char *line, int *list, int *nlist, int *limlist, int linelen)

Fortran wrapper to parselist, called by the Fortran parselist2 subroutine. The list is returned into array list and limlist is the size of that array, or 0 if unknown.  The return value is 1 for a memory allocation error, 2 for a bad character in the entry, and -1 for the list too big for the array; in the latter case the portion of the list that fits is returned but parselist2 will exit with an error.  If a slash is entered, the input list is returned unchanged.

Histogram Functions

void kernelHistogram(float *values, int numVals, float *bins, int numBins, float firstVal, float lastVal, float h, int verbose)

Computes a standard or kernel histogram from the numVals values in the array values.  The histogram is placed in the array bins, and occupies numBins between firstVal and lastVal. The kernel half-width is set by h, where a triweight kernel equal to (1-(x/h)
2)
3 / (0.9143 * h) is added at each point.  Use 0 for a standard binned histogram.  Set verbose to 1 for a list of values, or 2 for a list of bin values.

void kernelhistogram(float *values, int *numVals, float *bins, int *numBins, float *firstVal, float *lastVal, float *h, int *verbose)

Fortran wrapper for kernelHistogram

int scanHistogram(float *bins, int numBins, float firstVal, float lastVal, float scanBot, float scanTop, int findPeaks, float *dip, float *peakBelow, float *peakAbove)

Scans a histogram for peaks and a dip.  The histogram is in numBins elements of bins, extending from firstVal to lastVal.  It will be scanned between values scanBot and scanTop.  If findPeaks is nonzero, it will find the two highest peaks and return their values in peakBelow and peakAbove, then find the dip between them and returns its value in dip; otherwise it just scans the range for the lowest point and returns the value in dip.  The return value is 1 if it fails to find two peaks.

int scanhistogram(float *bins, int *numBins, float *firstVal, float *lastVal, float *scanBot, float *scanTop, int *findPeaks, float *dip, float *peakBelow, float *peakAbove)

Fortran wrapper for scanHistogram

int findHistogramDip(float *values, int numVals, int minGuess, float *bins, int numBins, float firstVal, float lastVal, float *histDip, float *peakBelow, float *peakAbove, int verbose)

Finds a histogram dip by starting with a high smoothing and dropping to a lower one.  
values - Array of values to form histogram from
numVals - Number of values in the array
bins - Array to build histogram in
numBins - Number of bins to divide histogram into
firstVal, lastVal - Starting and ending value for histogram range
minGuess - If non-zero, it specifies an estimate of the minimum number of items above the dip
verbose - Verbose output flag, passed to kernelHistogram
histDip - Returned with the location of the dip
peakBelow, peakAbove - Returned with the peak values below and above the dip
The return value is 1 if a dip cannot be found after 4 trial H values.

int findhistogramdip(float *values, int *numVals, int *minGuess, float *bins, int *numBins, float *firstVal, float *lastVal, float *histDip, float *peakBelow, float *peakAbove, int *verbose)

Fortran wrapper for findHistogramDip

Functions for Analyzing Gold Beads

void makeModelBead(int boxSize, float beadSize, float *array)

Makes a model bead of radius beadSize in array with dimensions boxSize by boxSize.

void makemodelbead(int *boxSize, float *beadSize, float *array)

Fortran wrapper for makeModelbead.

double beadIntegral(float *array, int nxdim, int nx, int ny, float rCenter, float rInner, float rOuter, float xcen, float ycen, float *cenmean, float *annmean, float *temp, float annPct, float *median)

Finds integral above background of a bead located at xcen, ycen in array, an image nx by ny with X dimension nxdim.  rCenter, rInner, and rOuter are radii of the center and the inner and outer radii of the background annulus.  Returns the center and annular means in cenmean and annmean, and if annPct is supplied with a fractional percentile, returns the percentile value in median, using temp as a temporary array.  temp is not used if annPct is 0.  The return value is center minus annular mean.

Functions for Image Statistics by Sampling

int sampleMeanSD(unsigned char **image, int type, int nx, int ny, float sample, int ixStart, int iyStart, int nxUse, int nyUse, float *mean, float *sd)

Estimates mean and SD of a sample of an image.  Returns nonzero for errors.
image = array of pointers to each line of data
type = data type, 0 for bytes, 2 for unsigned shorts, 3 for signed shorts, 6 for floats, 7 for 4-byte integers, 8 for RGB triples of bytes or 9 for RGBA bytes; RGB values are weighted with NTSC weighting
nx, ny = X and Y dimensions of image
sample = fraction of pixels to sample
ixStart, iyStart = starting X and Y index to include
nxUse, nyUse = number of pixels in X and Y to include
mean, sd = returned values of mean and sd

int samplemeansd(float *image, int *nx, int *ny, float *sample, int *ixStart, int *iyStart, int *nxUse, int *nyUse, float *mean, float *sd)

Fortran wrapper for sampleMeanSD with a floating point array in image

int percentileStretch(unsigned char **image, int type, int nx, int ny, float sample, int ixStart, int iyStart, int nxUse, int nyUse, float pctLo, float pctHi, float *scaleLo, float *scaleHi)

Estimates percentile limits for contrast strecthing an image by building a histogram of a sample.
image = line pointers to image
mode = data type, a SLICE_MODE_ value; BYTE, USHORT, SHORT, or FLOAT are allowed
nx, ny = X and Y dimensions of image
sample = fraction of pixels to sample
ixStart, iyStart = starting X and Y coordinates to sample
nxUse, nyUse = extent of image to sample in X and Y
pctLo, pctHi = Percentile levels to find on the low and high end of range (1.0 is 1%)
scaleLo, scaleHi = returned values of those percentile levels
Returns 1 for inadequate number of pixels, 2 for invalid mode, or 3 for memory error.

Functions for Measuring Volume Mean/SD at Multiple Binnings

int multiBinSetup(int binning[][3], int boxSize[][3], int boxSpacing[][3], int numScales, int startCoord[3], int endCoord[3], int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds)

Does initial setup for determining the local mean and standard deviation in an array of boxes within an image volume at a series of binnings (scales).  For all of the array variables showing a dimension of 3, the index in the dimension is 0 for the X, 1 for the Y, and 2 for the Z value.  
Inputs:  
binning         - Binning in X, Y, and Z for each scale  
boxSize         - Size of box to analyze in X, Y, and Z for each scale, in binned pixels   
boxSpacing      - Spacing between adjacent boxes in X, Y, and Z for each scale, in binned pixels  
numScales       - Number of binnings to analyze  
startCoord      - Starting unbinned coordinate of volume to analyze in X, Y, and Z  
endCoord        - Ending coordinate of volume to analyze in X, Y, and Z (inclusive)
Outputs:  
boxStart        - Starting coordinate in binned slices of first box in X, Y, and Z for each scale (equals 0 at start of slice, not startCoord)  
numBoxes        - Number of boxes i nX, Y, and Z for each scale  
bufferStartInds - Array to receive starting index of the binned slice for each scale; must be dimensioned to at least numScales + 1.  The size of a floating point array to allocate for the loaded and binned slices is returned in the extra element.  
statStartInds   - Array to receive starting index of the box mean and SD values for each scale; must be dimensioned to at least numScales + 1.  The size of floating point arrays to allocate for the means and SDs is returned in the extra element.  
Returns 1 for inappropriate coordinate limits, 2 for more than the number of scales defined by MAX_MBS_SCALES, 3 for box size bigger than the binned size of the volume at some scale, or 4 for boxSpacing less than 1.

int multibinsetup(int binning[][3], int boxSize[][3], int boxSpacing[][3], int *numScales, int *startCoord, int *endCoord, int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds)

Fortran wrapper for multiBinSetup.  In Fortran, the 2-D arrays would all be dimensioned (3,*), with the X/Y/Z index in the first dimension.

int multiBinStats(int binning[][3], int boxSize[][3], int boxSpacing[][3], int numScales, int startCoord[3], int endCoord[3], int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds, float *buffer, float *means, float *SDs, int *funcData, int (*getSliceFunc)(int *, int *, float *))

Computes the local mean and standard deviation in an array of boxes within an image volume at a series of binnings.  binning, boxSize, boxSpacing, numScales, startCoord, endCoord, boxStart, numBoxes, bufferStartInds, and statStartInds must be the same variables supplied to multiBinSetup.  buffer is an array allocated to the size given by the value after numScales values in bufferStartInds.  Mean and standard deviation values are returned in means and SDs, which must be allocated to the value after numScales values in statStartInds.  The callback function getSliceFunc is called to get the one slice of data within the given coordinate range by getSliceFunc(&iz, funcData, buffer), where iz is the Z coordinate and funcData is filled with whatever the loading function requires (i.e., if it needs the starting and ending X and Y coordinates, they must be passed in funcData).  Returns an error code from getSliceFunc.

int multibinstats(int binning[][3], int boxSize[][3], int boxSpacing[][3], int *numScales, int *startCoord, int *endCoord, int boxStart[][3], int numBoxes[][3], int *bufferStartInds, int *statStartInds, float *buffer, float *means, float *SDs, int *funcData, int (*getSliceFunc)(int *, int *, float *))

Fortran wrapper to getSliceFunc

Functions for Working with Montage Piece Lists

int checkPieceList(int *pclist, int stride, int npclist, int redfac, int nframe, int *minpiece, int *npieces, int *noverlap)

Examines a list of npclist piece coordinates in pclist for one dimension (x or y) at a spacing of stride between successive values, applies a reduction factor redfac, and given the frame size nframe, it finds the minimum coordinate minpiece, the # of pieces npieces, and the overlap noverlap. It also checks that ALL coordinates on the list are multiples of frame size minus overlap.  Upon failure, it returns the line number at which that  check fails, and also return -1 in npieces.

void checklist(int *pclist, int *npclist, int *redfac, int *nframe, int *minpiece, int *npieces, int *noverlap)

Fortran wrapper for checkPieceList which assumes the stride is 1.  Upon error, it prints a message, and still returns -1 in npieces.

void adjustPieceOverlap(int *pclist, int stride, int npclist, int nframe, int minpiece, int noverlap, int newOverlap)

Adjusts a list of npclist piece coordinates in pclist for one dimension (x or y), at a spacing of stride between successive values, changing the overlap from noverlap to newOverlap, given the frame size nframe and the minimum coordinate minpiece (which is retained).

void adjustpieceoverlap(int *pclist, int *npclist, int *nframe, int *minpiece, int *noverlap, int *newOverlap)

Fortran wrapper for adjustPieceOverlap that assumes a stride of 1.

Functions for Finding Alignments of Montage Pieces

void montXCBasicSizes(int ixy, int nbin, int indentXC, int *nxyPiece, int *nxyOverlap, float aspectMax, float extraWidth, float padFrac, int niceLimit, int *indentUse, int *nxyBox, int *numExtra, int *nxPad, int *nyPad, int *maxLongShift)

Sets up most of the sizes for the overlap zone correlations.  Inputs:  
ixy           - 0 for X edge, 1 for Y edge
nbin          - Binning
indentXC      - Border from edge of image  
nxyPiece      - Size of pieces in X and Y
nxyOverlap    - Overlap between pieces in X and Y
aspectMax     - Maximum aspect ratio to correlate
extraWidth    - Fraction of overlap zone width to add for extra extent
padFrac       - Fraction of box size to pad on each edge
niceLimit     - Largest prime factor for FFT sizes
Outputs:  
indentUse     - Actual border being used
nxyBox        - Binned size  in X and Y of boxes to be extracted from images
numExtra      - Number of extra pixels in X and Y
nxPad         - Size of padded image in X  
nyPad         - Size of padded image in Y  
maxLongShift  - Maximum shift along the long axis of the overlap zone

void montxcbasicsizes(int *ixy, int *nbin, int *indentXC, int *nxyPiece, int *nxyOverlap, float *aspectMax, float *extraWidth, float *padFrac, int *niceLimit, int *indentUse, int *nxyBox, int *numExtra, int *nxPad, int *nyPad, int *maxLongShift)

Fortran wrapper for montXCBasicSizes.  ixy should be 1 or 2.

void montXCIndsAndCTF(int ixy, int *nxyPiece, int *nxyOverlap, int *nxyBox, int nbin, int indentUse, int *numExtra, int nxPad, int nyPad, int numSmooth, float sigma1, float sigma2, float radius1, float radius2, int evalCCC, int *ind0Lower, int *ind1Lower, int *ind0Upper, int *ind1Upper, int *nxSmooth, int *nySmooth, float *ctf, float *delta)

Sets up indices for extracting boxes, and the filter function.  Inputs: ixy           - 0 for X edge, 1 for Y edge
nxyPiece      - Size of pieces in X and Y
nxyOverlap    - Overlap between pieces in X and Y
nxyBox        - Binned size  in X and Y of boxes to be extracted from images
nbin          - Binning
indentUse     - Border from edge of image  
numExtra      - Number of extra pixels in X and Y
nxPad         - Size of padded image in X  
nyPad         - Size of padded image in Y  
numSmooth     - Number of pixels of smoothing if possible  
sigma1        - Sigma for inverted gaussian at low frequency
sigma2        - Sigma for high frequency filter
radius1       - Radius out to which low frequency filter is 0
radius2       - Radius for high frequency filter
evalCCC       - 1 to evaluate correlation coefficients
Outputs:  
ind0Lower     - Starting index to extract in X and Y for lower piece
ind1Lower     - Ending index to extract in X and Y for lower piece
ind0Upper     - Starting index to extract in the ixy dimension for uppoer piece
ind1Upper     - Ending index to extract in the ixy dimension for uppoer piece
nxSmooth      - Size of edge-smoothed image in X
nySmooth      - Size of edge-smoothed image in Y
ctf           - Array of filter values; must be dimensioned at least 8193
delta         - Change in radius per index in the ctf array

void montxcindsandctf(int *ixy, int *nxyPiece, int *nxyOverlap, int *nxyBox, int *nbin, int *indentUse, int *numExtra, int *nxPad, int *nyPad, int *numSmooth, float *sigma1, float *sigma2, float *radius1, float *radius2, int *evalCCC, int *ind0Lower, int *ind1Lower, int *ind0Upper, int *ind1Upper, int *nxSmooth, int *nySmooth, float *ctf, float *delta)

Fortran wrapper for montXCIndsAndCTF.  ixy should be 1 or 2.

int montXCFindBinning(int maxBin, int targetSize, int indentXC, int *nxyPiece, int *nxyOverlap, float aspectMax, float extraWidth, float padFrac, int niceLimit, int *numPaddedPix, int *numBoxedPix)

Finds binning needed to keep boxed out area smaller than a target size.  Inputs: maxBin        - Maximum binning to select
targetSize    - Target size to make the number of boxed pixels be less than the square of
indentXC      - Border from edge of image to assume
nxyPiece      - Size of pieces in X and Y
nxyOverlap    - Overlap between pieces in X and Y
aspectMax     - Maximum aspect ratio to correlate
extraWidth    - Fraction of overlap zone width to add for extra extent
padFrac       - Fraction of box size to pad on each edge
niceLimit     - Largest prime factor for FFT sizes
Outputs:  
numPaddedPix  - Total number of pixels for padded images, with extra allowance
numBoxedPix   - Total number of pixels for boxed images, with extra allowance
The return value is the binning.

int montxcfindbinning(int *maxBin, int *targetSize, int *indentXC, int *nxyPiece, int *nxyOverlap, float *aspectMax, float *extraWidth, float *padFrac, int *niceLimit, int *numPaddedPix, int *numBoxedPix)

Fortran wrapper for montXCFindBinning.

void montXCorrEdge(float *lowerIn, float *upperIn, int *nxyBox, int *nxyPiece, int *nxyOverlap, int nxSmooth, int nySmooth, int nxPad, int nyPad, float *lowerPad, float *upperPad, float *lowerCopy, int numXcorrPeaks, int legacy, float *ctf, float delta, int *numExtra, int nbin, int ixy, int maxLongShift, int weightCCC, float *xDisplace, float *yDisplace, float *CCC, void (*twoDfft)(float *, int *, int *, int *), void (*dumpEdge)(float *, int *, int *, int *, int *, int *), char *debugStr, int debugLen, int debugLevel)

Performs Fourier correlations and evaluation of real-space correlations to find displacement on an edge.  Whether the CCC's are evaluated depends on the value of numXcorrPeaks.  If it is 1, they are not computed, but 16 peaks are found and the first peak is taken that does not exceed the lateral shift in maxLongShift.  If it is greater than 1, at least 16 peaks are still found and checked against maxLongShift, then the CC is evaluated at the given number of peaks. Inputs: lowerIn       - Boxed, possibly binned image extracted from lower piece
upperIn       - Boxed, possibly binned image extracted from upper piece
nxyBox        - Size of boxed input images in X and Y
nxyPiece      - Size of full pieces in X and Y
nxyOverlap    - Overlap between pieces in X and Y
nxSmooth      - Size of edge-smoothed image in X
nySmooth      - Size of edge-smoothed image in Y
nxPad         - Size of padded image in X  
nyPad         - Size of padded image in Y  
lowerPad      - Temporary array for lower padded piece
upperPad      - Temporary array for upper padded piece
lowerCopy     - Temporary array for copy of lower padded piece, needed only if evaluating CCC's at multiple peaks
numXcorrPeaks - Number of correlation peaks to examine
legacy        - Non-zero to do legacy correlations
ctf           - Array with filter function
delta         - Change in radius per index in the ctf array
numExtra      - Number of extra pixels in X and Y
nbin          - Binning
ixy           - 0 for X edge, 1 for Y edge
maxLongShift  - Maximum shift along the long axis of the overlap zone
weightCCC     - Non-zero to weight CCC by power function of overlap area
Outputs, functions, debugging variables:
xDisplace     - Displacement in X: amount to shift upper piece to align to lower
yDisplace     - Displacement in Y
CCC           - Value of CCC at the chosen peak, if computed
twoDfft       - 2-D FFT function to receive array, nxPad, nyPad, and 0 for forward or 1 for inverse
dumpEdge      - If non-null, a function for writing images to file.  Its arguments are the image array, X dimension, X size, Y size, ixy + 1,  and 0 for an image or 1 for a correlation  
debugStr      - String to receive debugging text.  To hold all text with full output, it should be numXcorrPeaks + 1 times MONTXC_MAX_DEBUG_LINE
debugLen      - Length of debugLen
debugLevel    - 1 for summary of selected peak and the displacement, 2 for full output about each peak, 3 to include eliminated peaks

void montxcorredge(float *lowerIn, float *upperIn, int *nxyBox, int *nxyPiece, int *nxyOverlap, int *nxSmooth, int *nySmooth, int *nxPad, int *nyPad, float *lowerPad, float *upperPad, float *lowerCopy, int *numXcorrPeaks, int *legacy, float *ctf, float *delta, int *numExtra, int *nbin, int *ixy, int *maxLongShift, int *weightCCC, float *xDisplace, float *yDisplace, float *CCC, void (*twoDfft)(float *, int *, int *, int *), void (*dumpEdge)(float *, int *, int *, int *, int *, int *), int *debugLevel)

Fortran wrapper for montXCorrEdge.  The wrapper will collect and print the debug strings if debugLevel is non-zero.

void montBigSearch(float *array, float *brray, int nx, int ny, int ixBox0, int iyBox0,int ixBox1, int iyBox1, float *dxMin, float *dyMin, float *sdMin, float *ddenMin, int numIter, int limStep)

Compares images in two arrays array and brray (dimensioned nx, ny, with image dimensions nx by ny), and finds the displacements between the two images with the minimum standard deviation of the pixel-by-pixel difference. ixBox0, iyBox0, ixBox1, and iyBox1 (numbered from 0) define the array index limits of the box to be compared in ARRAY.  dxMin and dyMin, upon entering, should have the starting displacement between the box position in array and the position to be compared in brray.  Initially it will search with a step size of 1, up to a distance of limStep from the initial displacements.  It will do numIter iterations, cutting the step size by a factor of 2 on each iteration.  sxMin and dyMin will be returned with the displacement that has the minimum standard deviation sdMin, and ddenMin is the difference in density at that displacement (B-A).

void montbigsearch(float *array, float *brray, int *nx, int *ny, int *ixBox0, int *iyBox0,int *ixBox1, int *iyBox1, float *dxMin, float *dyMin, float *sdMin, float *ddenMin, int *numIter, int *limStep)

Fortran wrapper to montBigSearch, where the array index limits are numbered from 1.

void montSdCalc(float *array, float *brray, int nx, int ny, int ixBox0, int iyBox0, int ixBox1, int iyBox1, float dx, float dy, float *sd, float *dden)

Compares positions in two image arrays array and brray (each dimensioned to nx, ny, the image dimensions as well), and within a defined box, calculates the standard deviation sd of the pixel-by-pixel difference (B minus A).  ixBox0, iyBox0, ixBox1, and iyBox1 (numbered from 0) define the array index limits of the box in array, and dx and dy define the displacement between that position and the desired position in brray. If dx and dy are integers, the computation is done rapidly between corresponding pixels; otherwise, quadratic interpolation is used to find the pixel value in brray for a given pixel in array.

void montsdcalc(float *array, float *brray, int *nx, int *ny, int *ixBox0, int *iyBox0, int *ixBox1, int *iyBox1, float *dx, float *dy, float *sd, float *dden)

Fortran wrapper to montSdCalc, where the array index limits are numbered from 1.

int findPieceShifts(int *ivarpc, int nvar, int *indvar, int *ixpclist, int *iypclist, float *dxedge, float *dyedge, int idir, int *pieceLower, int *pieceUpper, int *ifskipEdge, int edgeStep, float *dxyvar, int varStep, int *edgeLower, int *edgeUpper, int pcStep, int *work, int fort, int leaveInd, int skipCrit, float robustCrit, float critMaxMove, float critMoveDiff, int maxIter, int numAvgForTest, int intervalForTest, int *numIter, float *wErrMean, float *wErrMax)

Solves for positions of overlapping pieces by adjusting them iteratively to minimize the difference between measured displacement between pairs of pieces and the difference in their positions.
ivarpc: index from included variables (pieces) to overall piece number
nvar: Number of variables
indvar: index from piece number to variable number; must be set to less than the minimum variable number (0 for C, 1 for Fortran) for all pieces potentially connected to the ones included as variables.
ixpclist, iypclist: X and Y piece coordinates.
dxedge, dyedge: the measured difference from the nominal displacement between pieces at an edge
idir: 1 if those values are the amount that the upper piece is displaced from being in register with the lower; -1 if they are the amount that the upper piece needs to be shifted to be in register.
pieceLower, pieceUpper: index of pieces below and above an edge
ifskipEdge: A number indicating whether the edge should be skipped
edgeStep: Stride between X and Y edges at the same edge index in the edge arrays dxedge, dyedge, pieceLower, pieceUpper, ifskipEdge.  If the arrays are accessed in Fortran as array(edge,ixy), then this value is the first dimension of the array.  If they are accessed in C as array[edge][ixy] or array[2*edge+ixy] then the value is 1.
dxyvar: array into which the piece shifts in X and Y are returned for each piece in the variable list
varStep: Stride between X and Y shifts for a given piece in dxyvar; as for edgeStep, this is either 1 or the long dimension of the array.
edgeLower, edgeUpper: edge index of edge below or edge above a given piece in X or Y.
pcStep: Stride between X and Y edges at the same piece in those arrays; again, this is either 1 or the long dimension of the array
work: Array for temporary storage, must be at least 20.25 times the number of variables plus 1.
fort: A nonzero value indicates that the indexes in ivarpc, indvar, pieceLower, pieceUpper, edgeLower, edgeUpper are numbered from 1 instead 0.
leaveInd: Absolute edge index in edge arrays of edge to leave out
skipCrit: An edge is skipped if ifskipEdge is >= this value
critMaxMove: Terminate if the maximum move of any piece in X or Y is less than this criterion.
robustCrit: Use robust fitting to eliminate or down-weight outlying edge displacements; multiply outlier criterion by given factor.
critMoveDiff: Terminate if the average move changes by less than this amount between samples
maxIter: Maximum number of iterations
numAvgForTest: Number of iterations over which to average the average  
moves in X and Y for the test based on critMoveDiff  
intervalForTest: Number of iterations between such tests  
numIter is returned with the number of iteratation.
wErrMean and wErrMax are returned with the mean and maximum of weighted errors.
Returns 1 if the strides are not all 1 or all different from 1.
Called from Fortran with identical arguments and name.

Function for Computing a Scaled, Reduced Power Spectrum

int spectrumScaled(void *image, int type, int nx, int ny, void *spectrum, int padSize, int finalSize, int bkgdGray, float truncDiam, int filtType, void (*twoDfft)(float *, int *, int *, int *))

Computes a symmetrical power spectrum of the image in image, with slice type in type and size in nx and ny.  The image will padded to the square size in padSize, which should be a good size for taking a FFT, and tapered down to its mean inside the edges.  The FFT will be taken, the logarithm of the magnitudes will be computed and scaled to occupy a range of 0 to 32000, and this spectrum will be mirrored through the origin.  If finalSize is less than padSize, the spectrum will be reduced with antialiased reduction, using the filter type in filtType. If bkgdGray is less than or equal to 0, the result will be placed in spectrum as an array of shorts.  If bkgdGray is greater than 0, the spectrum will be scaled to range from 0 to 255, with the mean near the edges scaled to a value of bkgdGray, and the mean of pixels in a ring with diameter truncDiam as a fraction finalSize scaled to 255.  This result will be returned in spectrum as an array of bytes.  twoDfft must be supplied with the name of a 2-D FFT function with that calling sequence, (todfft if using libifft in IMOD).  
Returns -1 or -2 for the negative of an error in calling selectZoomFilter, 1 to 5 for an error in calling zoomWithFilter, -3 for a memory error, or -4 for finalSize greater than padSize, bkgdGray more than 192, or truncDiam negative or more than 0.75.

Color Map Functions

int *cmapStandardRamp()

Returns a pointer to the color ramp specification for the standard false color table used in 3dmod.

int *cmapInvertedRamp()

Returns a pointer to a color ramp specification for an false color table inverted from the one used in 3dmod.

int cmapConvertRamp(int *rampData, unsigned char table[3][256])

Converts a color ramp specification in rampData to a full 256-color table in table.  The first value in the rampData array is the number of color specifications; each set of 4 numbers after that is a red, green, and blue value (range 0-255) and a relative position along the color table.

int cmapReadConvert(char *filename, unsigned char table[3][256])

Reads a color map specification from the file specified by filename and converts it to a complete color table in table.  The file starts with the number of lines of color data.  If there are 256 lines they must be red, green, blue triplets (range 0-255); otherwise each line has a red, green, and blue followed by a relative position in the color table. Returns 1 for error opening a file, 2 for errors reading the file, 3 for an invalid number of lines, or 4 for a memory allocation error.

Standard color ramps

The standard color ramp specification used in 3dmod
static int standardRampData[] =
  { 
    14,
    100,   75,  200,  -530,
    60,    96,  255,  -469,
    0,    175,  177,  -400,
    0,    191,  143,  -383,
    0,    207,   78,  -361,
    90,   255,   60,  -305,
    191,  255,    0,  -259,
    239,  255,    0,  -240,
    255,  255,    0,  -229,
    255,  175,    0,  -172,
    255,  105,    0,  -122,
    255,   45,   55,   -60,
    255,    0,   90,   -40,
    255,    0,  192,    10
  };

/* An inverted color ramp specification, used in model view for mesh normals */
static int invertedRampData[] =
  { 
    14,
    255,    0,  192,  -10,
    255,    0,   90,   40,
    255,   45,   55,   60,
    255,  105,    0,  122,
    255,  175,    0,  172,
    255,  255,    0,  229,
    239,  255,    0,  240,
    191,  255,    0,  259,
    90,   255,   60,  305,
    0,    207,   78,  361,
    0,    191,  143,  383,
    0,    175,  177,  400,
    60,    96,  255,  469,
    100,   75,  200,  530
  };