* ImgLib.c * Please note last saved date below:* $Date: 2008/05/16 17:19:08 $ * A set of primitives for use in image processing. * c 1996-2002 Mark S. Cohen * This software is distributed as is with no guarantees. * Please report any errors or bug fixes to mscohen@ucla.edu * * For academic use only. Commercial users should contact * Mark Cohen for licensing information. * * Keep track of the differences between vlen and npts as the description * of vector sizes. While vlen is supposed to be used to represent the * length of the input vector, npts is to represent the number of points * actually processed. Use might not be entirely consistent. * * Most of the routines that ask for an output vector cannot be used as * in-place operations. * * ORIGINAL REVISION DATE: 9/15/96 * MAJOR REVISION: 9/21/96 * BACK CONVERSION TO UNIX 3/24/98 * *********************************************************************** * Keep control of ImgLib revisions * void GetImgLibVer( char *buff ) * * Returns the version number for ImgLib in buff * *********************************************************************** * OpenProcFile * * Open a standard proc file for tracking user inputs and settings * * TYPICAL CALLING SEQUENCE * * static char id[] = "$Revision: 2.96 $$Date: 2008/05/16 17:19:08 $"; * error = OpenProcFile( argc, argv, outFilename, id, &ProcFile ); * * INPUTS: argc -- from main command line * argv -- from main command line * outFilename -- basename of output file (typical) * ProcFile -- the freshly opened ProcFile * * FUTURE: capture command line in users home directory ************************************************************************ * ReportFileType * Return the name of the input file type, returning it as a text string * * char *ReportFileType( int file_type ) * ************************************************************************ * xinterp and yinterp * Bilinear interpolation of a 2 D image along a single dimension * Only integer interpolation ratios are supported. * * CALLING SEQUENCE: yinterp(*invec, xs, ys, Ratio, *outvec, data_type); * xinterp(*invec, xs, ys, Ratio, *outvec, data_type); * * INPUTS: invec -- vector of arbitrary data type * xs -- number of points in x dimension (input) * ys -- number of points in y dimension (input) * Ratio -- amount by which the image should be enlarged * data_type -- standard UC data type * OUTPUTS: None * outvec -- output vector of same data type as input * * Note that xinterp and yinterp are the same, but for a swap of x & y variables. * finterpx and finterpy are floating point versions ***************************************************************** * interp2D * 2-dimensional (x and y) interpolation for images * * CALLING SEQUENCE: * int interp2D( * void *invec // pointer to input data * int xsize, int ysize // image dimensions (input) * int Ratio // interpolation ratio - integers only * void *outvec // pointer to output data. Calling routine * must allocate sufficient space... * int data_type ) // standard UC type T_FLOAT, T_SHORT & T_USHORT supported * ***************************************************************** OSErr finterpX(void *invec, int xs, int ys, int outxs, void *outvec, int data_type) ***************************************************************** OSErr finterpY(void *invec, int xs, int ys, int outys, void *outvec, int data_type) ***************************************************************** OSErr finterp2D(void *invec, int xs, int ys, int outxs, int outys, void *outvec, int data_type) ***************************************************************** * pixreplicate * 2-dimensional (x and y) pixel replication for images * * CALLING SEQUENCE: * int pixreplicate( * void *invec // pointer to input data * int xsize, int ysize // image dimensions * int Ratio // interpolation ratio - integers only * void *outvec // pointer to output data. Calling routine * // must allocate sufficient space... * int data_type ) // standard UC type * ***************************************************************** * pixrepX * 1-dimensional (x) pixel replication for images * * CALLING SEQUENCE: * int pixrepX( * void *invec // pointer to input data * int xsize, int ysize // image dimensions * int Ratio // interpolation ratio - integers only * void *outvec // pointer to output data. Calling routine * // must allocate sufficient space... * int data_type ) // standard UC type * ***************************************************************** * pixrepY * 1-dimensional (y) pixel replication for images * * CALLING SEQUENCE: * int pixrepY( * void *invec // pointer to input data * int xsize, int ysize // image dimensions * int Ratio // interpolation ratio - integers only * void *outvec // pointer to output data. Calling routine * // must allocate sufficient space... * int data_type ) // standard UC type * ***************************************************************** * DownSample * 2-dimensional matrix reduction for images * * CALLING SEQUENCE: * OSErr DownSample( * void *invec // pointer to input data * int xsize // image dimensions for original image * int ysize * int xFactor, yFactor reduction ratio - positive integers only * void *outvec // pointer to output data. Calling routine * // must allocate sufficient space... * int data_type ) // standard UC type * ********************************************************************** * ScaleImageSize * Interpolate an image and place it into newly created memory. Calling routine * should free this memory as soon as possible * * Requires: * float *ScaleImageSize( * float *theImage source image * IMAGE im descriptive header info * int interpFactor interpolation factor: Positive number for increase, negative for decrease * int displayMode how should interpolation be done? * *OSErr error error return value * int data_type ) * * Changes: Allocates memory for interpbuf and returns image in it. ***************************************************************** * ImageFlipY * Flip image top and bottom * OK to use in-place.. * * CALLING SEQUENCE: * OSErr ImageFlipY( * void *invec // pointer to input data * int xsize // image dimensions for original image * int ysize * void *outvec // pointer to output data. Calling routine * // must allocate sufficient space... * int data_type ) // standard UC type * *******************************************************************\ * CropImageSides * * Remove the specified number of points from the left and right sides * of an image, returning the data in the first part of the input buffer * * Calling routine is responsible for disposing of, or resizing memory * * void CropImageSides ( * void *ImageData // the original data * short UnCropdXSize // number of pixels per x line before cropping * short ysize // number of y lines * short LeftSkip // points to skip on image left * short RightSkip // points to skip on image right * short bytesPix ) // number of bytes/pixel * * CHANGES: * Contents of image data ************************************************************************ * xcrop - remove the left and right quarter images from a large image * * CALLING SEQUENCE: * xcrop( * *InImg, Input Image vector * xs, ys, x and y dimensions in ints * *CropImg, Output Image vector * data_type ) Standard UC data type * ******************************************************************************* * Overlay * * merge two float images to make a uchar composite * * CALLING SEQUENCE: * Overlay( float *baseImage, Image to use as base. Must already be * rescaled to range from BMAX to BMIN * float *OverlayImage, Image to overlay onto base * unsigned char *composite, Output * long imSize, * float overlay_max, Maximum value to be used from overlay * float overlay_cutoff, Minimum value to be used from overlay * short overlayMode ) Pos/Neg/Whatever... * ******************************************************************************* * vGaussSmooth * Use convolution to smooth using a Gaussian kernel * Accepts and returns only T_FLOAT at this time * * CALLING SEQUENCE: * vGaussSmooth( void *invec -- only T_FLOAT now * int width -- Width of Gaussian kernel in pixels * int step -- Step size through input vector * long npts -- Number of points to process * int data_type ) -- Only Accepts T_FLOAT * OUTPUTS: * This is an in-place routine. returns the smoothed vector * ************************************************************************ * HanningSmooth * Performs a very rapid in place 2-dimensional convolution smoothing with a Hanning Kernel * * The Hanning filter is a convolution with a 0.25, 0.5, 0.25 kernel. We handle this * as integer math, by using a 1-2-1 kernel and dividing the sum by four. * * CALLING SEQUENCE: * HanningSmooth( short Kernel, // smoothing kernel * unsigned short *image, // input data * unsigned short *smoother, // storage buffer as large as image * int xsize, * int ysize, * int data_type ) * * INPUTS: *image -- pointer to a 2D image array of short ints * xsize, ysize -- x and y dimensions * Kernel == kSmooth9 or kSmooth25 to control degree of smoothing. * Larger smoothing radii should be performed in k-space * * OUTPUTS: Smoothed image is returned in place. * ************************************************************************ * CalcT2 * * Determine the T2 and Rho of a location based on intensity values and * return the Chi2 value for goodness of fit. * * Data dimensions are passed in the IMAGE struct * * CALLING SEQUENCE * OSErr CalcT2( IMAGE *im, // struct with image dimensions. data_type is IGNORED. * float *data, // set of image volumes arranged as time points in order of te * char *Mask, // binary mask indicating which locations are to be fitted * float *TEs, // vector of te values in msec * float *T2, // image vector for T2 outputs * float *Rho, // image vector for Rho outputs * float *chi2 ) // image vector for fit values * * RETURNS: error code * ************************************************************************ * imageDataRange * * Determine the noise level and lower fraction of the nonNoise data * range based on a histogram of pixel intensities * * Assume that the noise is a commonly occuring value in the bottom * NRANGE of values. Use a histogram of HBINS size and threshold at NMULT * times the peak in the lower NRANGE. * * If the data ranges to negative values, simply return the max and min. * * CALLING SEQUENCE * OSErr imageDataRange( void *vec, * long npts, * void *hicutoff, * void *noiseLevel, * long *nonNoise, * int data_type ) * * INPUTS: * vec -- input vector * npts -- vector length (int) * data_type -- output vector of char of same length as input * OUTPUTS: * noiseLevel is updated to contain the noise threshold * hiCutoff contains the level below which fraction of pixel intensities are found * nonNoise is the number of points above threshhold * RETURNS: error code * ************************************************************************ * ThreshMask & vThreshold * Form an image mask having values of 0 for all pixels below threshold and 1 for all * pixels greater than the threshold * * Input must be of type unsigned short. * * CALLING SEQUENCE * long ThreshMask( void *myImage, // form a mask of points above threshold * unsigned char *theMask, * long imageSize, * void *threshold, * int data_type ) * * long vThreshold( void *myImage, // zero out points below threshhold * long vlen, * void *threshold, * int data_type ) * * INPUTS: *myImage -- image vector * imageSize -- size of image * threshold -- value above which points are treated as valid * * OUTPUTS: *theMask -- vector of unsigned char representing the image mask (ThreshMask) * vThreshold zeroes all points in vec greater than threshold * returns number of points over threshold. ************************************************************************ * autoMask & autoThresh * * Determine the noise level and threshold based * on histogram determination * * Assume that the noise is a commonly occuring * value in the bottom NRANGE of values. Use a * histogram of HBINS size and threshold at * NMULT times the peak in the lower NRANGE. * * CALLING SEQUENCE autoMask: // create a mask of pixels over a calculated threshhold * OSErr autoMask( void *invec, * unsigned char *theMask, * long npts, * void *threshold, * long *nonNoise, * int data_type ) * * INPUTS: * invec -- image vector * mask -- output vector of short of same length as input * npts -- vector size * vlen -- vector length (long) * OUTPUTS: * theMask -- pre-allocated vector that will contain 1 for points above thresh and 0 for points below * threshold -- updated to contain the value used * nonNoise -- the number of points above threshold * * CALLING SEQUENCE autoThresh: // determine the noise threshold * OSErr autoThresh( void *invec, * long npts, * void *threshold, * long *nonNoise, * int data_type ) * * INPUTS: * invec -- input vector. T_USHORT and T_FLOAT supported * npts -- vector length (int) * OUTPUTS: * threshold -- contains the calculated value * nonNoise -- the number of points above threshold ************************************************************** * FloatToCImage * * Convert a floating point vector to a vector of char * * void FloatToCImage( * inImage -- pointer to the input image * cimage -- pointer to the output (char) image * len -- the number of image points * *max, *min -- Calculated or input max and min. Same type as inImage * rescale -- kAutoRange, kFileMinMax, kManualRange, kNoScale * data_type )-- data type of input image * ********************************************************************** * * ImageChangeDetect * * Count the number of image points that change from above to below * threshold or vice versa. * * OSErr ImageChangeDetect( * INPUTS * myImage -- pointer to the input image * theMask -- pointer to the previously calculated image mask * threshold -- intensity threshold used to make the mask * npts -- number of points in the image and mask * RETURNS * numlowerPts -- number of points that have dropped below threshold * numGreaterPts -- number of points that have increased above threshold * ********************************************************************** * * MaskChangeDetect * * Count the number of points that differ between two masks. * * long MaskChangeDetect( * INPUTS * oldMask -- pointer to the previously calculated image mask * newMask -- pointer to the current image mask * npts ) -- number of points in the image and mask * RETURNS * number of points that have changed * ************************************************************************ * vRescale * Adjust the values of the input vector to fall between newMin and newMax * * CALLING SEQUENCE: * double vRescale( void *invec, -- input (and output) data vector * long vlen, -- length of input and ouput vectors * void *oldMax, -- Maximum value of input vector to be used * void *oldMin, -- Minimum value of input to be used * void *newMax, -- New Maximum for data range * void *newMin, -- New Minimum for data range * int data_type ); -- standard data type * * RETURNS: * calculated scale factor. * * Values less than or equal to oldMin in the input vector are set to newMin * Those greater than or equal to oldMax are set to newMax * The rest are scaled linearly between these endpoints. ***************************************************************** * getRowsCols - calculate the number of rows and columns in a multidisplay * image of known xs, ys and slice * * CALLING SEQUENCE: * getRowsCols(short slices, -- Number of slices in input volume * short &Rows, * short &Cols ) -- Returns number of rows and cols in multi display ************************************************************************* * mdisplay - take a vertical stack of images and rearrange into a single * multi-slice display. Memory for the output image must be * pre-allocated. Typically, this will mean a call to getRowsCols * * e.g., getRowsCols( slices, &rows, &cols ); * multImg = malloc( rows*cols*dataSize ); * error = mdisplay( multImg, in, xs, ys, slices, rows, cols, data_type ); * * Returns an error if the number of slices does not fit into the Row x Col space * CALLING SEQUENCE: * OSErr mdisplay( * OutImg -- Storage for output image * InImg -- Vertically arranged image stack or volume * xs, ys -- Dimension of slices in input volume * slices -- Number of slices in input volume * data_type -- Standard UC data_type for input and output * ) * NOT TO BE USED AS AN IN PLACE ROUTINE! ****************************************************************** * unmdisplay * Take multi-display data set and convert it into a * 'volume' structure (xs x ys*slices ) * * OSErr unmdisplay( void *InImg, -- input image * void *OutImg, -- pre-allocated memory for output * int xs, int ys, -- image dimensions of output images * int slices, -- number of valid slices in the input set * int data_type ) * * xs and ys are the sizes of the final images (after unmdisplay) ************************************************************************ * waitfopen * Read successive items from a file. Waits for items to exist, if necessary. * * CALLING SEQUENCE: * FILE *waitfopen( char *name, -- the filename to open * char *mode -- the file mode * int WaitTime ) -- time to wait for file existence (seconds) ************************************************************************ * ReadNextFileString * * Sort of a stripped down fscanf for reading strings from Mac files * * OSErr ReadNextFileString( short fileID, char *myString ) * ************************************************************************ * WaitNextItemFromList * Read successive items from a (paragigm) file. Waits for items to exist, if necessary. * * CALLING SEQUENCE: * (unix) * OSErr WaitNextItemFromList( char *Item, -- the returned item (string) * FILE *fp, -- pointer to an open file * char *fname -- the filename for an error message * int WaitTime ) -- Time to wait for item to be valid ( seconds ); * (mac) * OSErr WaitNextItemFromList( char *Item, short fileID, int WaitSec ) * ************************************************************************ * get_para_file_size * get_paradigm * * OSErr get_para_file_size( char *fname ) * OSErr get_paradigm( char *fname, float *pv, int npts ) * * get_para_file_size -- returns the number of points in a * "paradigm file" of the type used to process scanSTAT images. * INPUTS: char *fname * OUTPUTS: return value is the number of elements. * * get_paradigm -- is usually called after get_para_file_size. * It takes the ascii input in the text file and converts it * to a vector of floats. * INPUTS: char *fname -- name of the paradigm file * float *pv -- a pointer to the vector that will hold the * paradigm file values. * int npts -- number of points in the paradigm file. ****************************************************************************** * irpStringFromIMAGE * Create a string containing the valid contents for an APD header file * These files end in .irp and have names that are four characters shorter * than the corresponding image file. * * Note that more than one image file can be supported by a single .irp * * char *irpStringFromIMAGE( IMAGE *im ) * INPUTS: * a pointer to a populated IMAGE struct * OUTPUTS * a freshly minted string. * Calling routine is responsible for disposing this memory *************************************************************************** * SetImDatatype * * Sets the data type for an image to the selected type * * Calling String: * SetImDatatype( short data_type // base name of the output image * IMAGE *im ); * * Indirect accessor to ensure data consistency. *************************************************************************** * WriteAnalyzeFile * * Save image file with the selected name * * Calling String: * OSErr = ( char *basename // base name of the output image * IMAGE *im, * void *imageData ); * * Output image type is determined from file name and data type. * Data are rescaled if needed *************************************************************************** * CreateHeader * * Write a valid header of the type associated with the selected file type * * Calling String: * OSErr = ( int FileType, // BSHORT, BUCHAR, BFLOAT, ANALYZE, EPI (Genesis not yet supported) * IMAGE im, * char *basename ) // base name of the output image * * CreateHeader will clear the cropping flags that would appear in Analyze header, * as we save the images in cropped format. ****************************************************************************** * FindOutputType * Use the file name and extension to determine the output file type * OSErr FindOutputType( char *fname, int *itsType, short *data_type ) ***************************************************************************** * dsrFromIMAGE * Populate all fields of a dsr from the IMAGE struct * * OSErr dsrFromIMAGE( IMAGE *im, dsr *theDSR ) * * INPUTS: * a pointer to a populated IMAGE struct * a pointer to a dsr; * OUTPUTS: * error * CHANGES: * contents of dsr ***************************************************************************** * VarianFromIMAGE * Populate all fields of a varian header from the IMAGE struct * * OSErr VarianFromIMAGE( IMAGE *im, dsr *theDSR ) * * INPUTS: * a pointer to a populated IMAGE struct * a pointer to a Varian Header; * OUTPUTS: * error * CHANGES: * contents of Varian Header ******************************************************************************** * WriteMGHImage * Mostly as a debugging utility, write one image to a named output * file and create the needed headers. * * At least on initial release, performs no datatype conversions and * only supports output to MGH file types of bshort, buchar and bfloat. * * CALLING SEQUENCE: * WriteMGHImage( * void *imData, -- pointer to image data location * char *filename, -- filename string (base name only) * short xs, -- image dimensions * short ys, * short zs, * int data_type ) -- data type of image data. Five types supported. * * OUTPUTS: * Creates an image file of appropriate type to represent data, and * a corresponding .hdr file. Currently files are converted as follows: * data_type file type * T_USHORT ---> bshort * this may cause trouble. * * T_SHORT ---> bshort * T_CHAR ---> buchar * T_UCHAR ---> buchar * T_FLOAT ---> bfloat * * There is no data range checking and signed integer types may not be displayed * or managed properly (This is a limitation of the MGH header system). ****************************************************************** * InitializeIm * * void InitializeIm( IMAGE *im ) * * Empty an IMAGE struct by clearing all fields then adding a few defaults ************************************************************************ * UC_Readheader and UC_Readimage tools * * Original: Jianxin Wang 2/11/96 * Modification date: 9/17/96 * Bug Fix: 10/4/96 - fix errors on T_FLOAT imports * Major re-write and conversion for Macintosh - M. Cohen 4/97 * fix header bug & (nearly a) complete rewrite. wls 10-11-1997 * * UC_Readimage and UC_Readheader allow the calling routine to import images * from any of the five standard data formats used at UCLA: * * XXX.MR -- GE Signa Genesis format * XXX_YYY.img (with XXX.irp header files) -- ANMR APD2 format * XXX.img (with XXX.hdr header files) -- Analyze 3D or 4D * XXX.buchar (and XXX.hdr) -- blocks of unsigned char * XXX.bshort (and XXX.hdr) -- blocks of unsigned shorts * XXX.bfloat (and XXX.hdr) -- blocks of floats * * The latter three, inherited from MGH, use simple ascii headers with the * format: * ysize, xsize, zsize, swap * Where ysize is the number of points/row, xsize is the number of points/column * and zsize is the number of images in the file. The swap byte is set to 0 for * images created locally and 1 for images on a machine that is byteswapped * relative to the local machine. * * UC_Readheader is used to fill out an IMAGE struct with needed dimensional information * UC_Readimage is used to read in an image. The image is returned to the pointer * passed as *image_data with a data type specified by im->data_type. Note that the * memory for the image must be malloc'd by the calling routine. Usually this means * that UC_Readimage is called first, and the needed memory is calculated from the * x, y and timePts fields of the IMAGE struct. * * The crop flag is used on ANMR images when the FOV is oblong. This is especially * relevant to head images on that system. * * If you extend this routine to include other image types, please let us * know at UCLA so that we can keep it up to date. * * #ifdef UNIX * OSErr UC_Readheader( char *fn, IMAGE *im ) * #endif * * #ifdef MAC * OSErr UC_Readheader( FSSpec *imFileFS, IMAGE *im, Boolean RT ) * #endif * ************************************************************* * ParseXXXHeader * inputs: headerData populated with raw header info * outputs: IMAGE struct with valid contents based on headerData * headerData must contain the complete contents of the corresponding * file type (MGH, Signa, etc...) * * OSErr ParseSignaHeader( void *headerData, IMAGE *im ) * OSErr ParseMGHHeader ( void *headerData, IMAGE *im ) * OSErr ParseANMRHeader ( void *headerData, IMAGE *im ) * OSErr ParseDSR ( char *headerData, IMAGE *im, short fileSize ) *