/*
medcorrect.c
Last modified 2 April 2003
For brief usage, type "medcorrect". For more extensive help and code
history, type "medcorrect --help" or see the strings usage_str and
help_str below.
See COMMENTS section at the end of source code for additional notes,
todo lists, etc.
TODO:
1) Pedestal - is it necessary? If so when?
2) Obliquity correction has been disabled, since
a - I haven't tested it
b - It's menu driven
c - I don't anticipate SAXS folks needing it (I'm probably wrong)
Note that although the structure oblique is passed to
intcorr_int, intcorr_int only uses one of oblique's fields,
namely oblique.satn_level, which I've alread made a separate
variable (oblique is defined in objects.h).
3) Versions for Finger Lakes, Gruner 1K detector? Perhaps there
will be a correct_base program that all of these other programs
will link to (use argv[0] to set the defaults?)
4) Incorporate Gil's multi-image dezingering code.
5) Mention DEBUG flag in help string probably can define with gcc
or make flag
*/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include "objects.h"
#define MAX_FILENAME_SIZE 100
/* Apparently have to define constants this (rather than const int
...) way if needed for passing 2D arrays to a function (we're not
doing that, however, since infilenames is now type char ** and uses
the space allocated for argv, rather than allocating its own space
as char [][MAX_FILENAME_SIZE]) */
/* change to 1 or issue -DDEBUG at gcc for lots of extra messages */
#ifndef DEBUG
#define DEBUG 0
#endif
#ifdef __GNUC__
const char usage_str[]="
Usage: medcorrect [i|id|d|di] -s[z] image_file(s) [-o output_filename]
[-t int_correct_file] [-x x_distortion_file] [-b back_file(s)]
[--help]
** At least two arguments are required: '-s[z]' and an image filename.**
** If present, [i|id|d|di] must come first. Other arguments in any order.**\n\n";
#else
const char usage_str[]="Usage: medcorrect [i|id|d|di] -s[z] image_file(s) [-o output_filename]\n [-t int_correct_file] [-x x_distortion_file] [-b back_file(s)] \n [--help]\n See source medcorrect.c for more extensive documentation\n";
#endif
#ifdef __GNUC__
const char help_str[]="
medcorrect.c: This program corrects x-ray image files taken by the
medoptics detector, but can be easily adapted for other
detectors.
FEATURES:
1) The images may be dezingered, background subtracted, distortion
and/or intensity corrected. Corrections are chosen by command
line options.
2) The background images may also be dezingered. This is always
done if more than one background image (eventually more than
two will be allowed, but not now)
3) Command line substitution can be used to specificy multiple
image files. See example (4)
NOTE: Although the program allows complete flexibility, the proper
correction order is determined when correction files are made (see
TV6 command 'cath'). For the medoptics detector intensity
correction file medi12b1.tif, this order is raw->dezing -> back
sub -> intensity correct -> distortion correct. If using the
intensity correction file moi12b1.tif, the order of intensity and
distortion correction are swapped.
Usage: medcorrect [i|id|d|di]
-s[z] image_file(s)
[-o output_filename]
[-t int_correct_file]
[-x x_distortion_file]
[-b back_file(s)]
[--help]
At least two arguments are required: '-s[z]' and an image filename.
If present, [i|id|d|di] must come first. Other arguments in any order.
OPTIONS:
i: intensity correction. uses default intensity correction unless
-t is also present
d: distortion correction. uses default files unless -x is also present
NOTE 1: Correction is done in the order specificied. Make sure that
the this order (id vs. di) matches the type of intensity
correction file being used. (For the Medoptics decector,
see C:\\perm\\README on the Medoptics computer).
-t: explicitly name intensity correction file (see NOTE 2)
-x: explicitly name only the x-distortion correction file. (see
NOTE 2) The y-distortion file must have the same name but with
'y' in place of the final appearance of 'x' in the filename ,
e.g. medxbin1.tif & medybin1.tif
-s[z] image_file(s): if more than one is specified, images are
either combined and dezingered (if z is present) or processed
sequentially. Current limit of 2 images for dezingering
-b background file(s). If more than one is specified, images are combined and
dezingered.
-o output_file_name. If present, this takes precedence. If not
present, then the output filename depends on the action. If
multiple images are being corrected but NOT dezingered, then
the name is generated from the input file by the insertion of
\"_c\" before the extension (e.g. fish.tif -> fish_c.tif, fish
-> fish_c). This is useful for correcting a whole sequence of
images that do not need to be dezingered in a single command
(see example below). If two images are being dezingered, the
default name medcorr.out is used unless an explicit filename
is specified. NOTE that medcorrect NEVER overwrites files. If
the output filename it wants to already exists, medcorrect
aborts with a message to that effect.
--help: Print this help message and exit
NOTE 2: If distortion or intensity correction files are not explicitly named:
1) Look for environment variable CCDCORR_FILES to determine directory in which
to look for the correction files named in (default to /usr/local/lib).
2) Filenames found in the CCDCORR_FILES directory are
SUPERCEDED by the contents of medcorrect.cfg, if found inthe
local directory. This file is written whenever intensity and
distortion correction are both performed
Successfully opened correction filenames are saved to
medcorrect.cfg which, if found in the working directory,
supercedes the default files in /usr/local/lib.
EXAMPLES:
1) Dezinger a pair of images:
>medcorrect -sz image1.tif image2.tif -o image1z2.tif
** or using tcsh command line substitution **
>medcorrect -sz image[12].tif -o image1z2.tif
2) Dezinger two images and backgrounds and subtract the
dezingered background from the dezingered image:
>medcorrect -sz image1.tif image2.tif -o image1z2.tif -b b1.tif b2.tif
** or using tcsh command line substitution **
>medcorrect -sz image[12].tif -b b[12].tif -o image1z2.tif
3) The whole shabang (dezinger image & backs, then subtract
background, then intensity and distortion correct).
>medcorrect id -sz image[12].tif -b b[12].tif
** output stored in medcorr.out **
4) Correct but don't dezinger many files at once.
>medcorrect id -s images*.tif -b backround.tif
Originally written by Sandor L. Barna Princeton University Physics
Department Biophysics Group December 1993
Marian Szebenyi (CHESS) made many modifications, including pedestal,
through 2003.
Converted from text-menu/prompt to UNIX command-line-option style
by Arthur Woll, March 2003
NOTE: The UNIX shell expands wildcard characters like *, ?, and []
before sending the command line to the program. See your sh or csh
manual for details. Turbo C++ and Borland C++ expand wildcard
characters if the file WILDARG.OBJ is linked with your program. See
the manual for details. Almost all UNIX commands use a standard
command-line format. This standard has carried over into other
environments. A standard UNIX command has the form:
command options file1 file1 file3
Please refer to TV6 (Gruner, Eikenberry and Tate) for more information
on the correction procedures.\n";
#else
const char help_str[]="medcorrect help string not available on this (non-GNU) system.\n See the source file or the file USAGE in the distribution directory for the full string (help_str[]).\n";
#endif
const char insert[]="_c"; /* For automatic filename generation */
const char dot = '.'; /* Put insert[] just before dot */
const char default_outname[]="medcorr.out";
const char help_flag[]= "--help";
main(int argc, char *argv[]) {
/* -------------------- DECLARATIONS -------------------- */
extern const char help_flag[];
/* extern const int MAX_FILENAME_SIZE; */
const int TV6_HEADER_BYTES = 4096;
/* Header size for all of the correction files */
int bsub_flag = 0;
int intens_flag = 0;
int dist_flag = 0;
int obliq_flag = 0; /* In case we want to introduce it... */
int dist_first = 0;
int explicit_int = 0;
int explicit_dist = 0;
int explicit_outname = 0;
int dezing_images = 0; /* started life as combine_images */
int dezing_back = 0;
int max_images = 1; /* This is only used if either images or
backgrounds will be dezingered. In that case this is the
larger number of images to be dezingered. Right now we can
only dezinger two -- see COMMENTS */
char *x_addr;
int bfilen = 0;
int infilen = 0;
char dxfilename[MAX_FILENAME_SIZE];
char dyfilename[MAX_FILENAME_SIZE];
char intfilename[MAX_FILENAME_SIZE];
FILE *distx=NULL;
FILE *disty=NULL;
FILE *intens=NULL;
float *dxcorr=NULL;
float *dycorr=NULL;
float *intcorr=NULL;
float *tempbuf=NULL; /* Needed for intensity and distortion correction */
char **bfilenames;
char **infilenames;
FILE *bfile=NULL;
FILE *infile=NULL; /* In main for checking existance & reading header
& byte order. Also local in read_data for
reading actual data */
unsigned short *indat=NULL;
unsigned short *bdat=NULL;
unsigned short *dat_ptr; /* Never to have its own allocation. */
unsigned short **image_stack = NULL; /* A 2-D array to be malloc'd
and to hold multiple data arrays in
the case of dezingering either the
background or input files */
char outname[MAX_FILENAME_SIZE] = "";
unsigned short *outdat;
char *header = NULL, *trailer = NULL;
int header_bytes = 0, trailer_bytes = 0;
int err = 0; /* Initialize to zero - be optimistic! */
int length;
int tiff_format;
int swizdetect_flag;
int *swizdetect;
int swizzle=0; /*--- swizzle = 0: little endian/PC byte order
swizzle = 1: big endian as for SGI ----*/
register i;
const char cffilename[] = "medcorrect.cfg";
FILE *config;
char cfilesdir[MAX_FILENAME_SIZE];
char *corrfiles;
char tempname[MAX_FILENAME_SIZE];
char suffix[8];
char hold_file[3][MAX_FILENAME_SIZE];
char yesno[1];
char response[2];
int pedestal = 0;
int saturated;
struct oblparams oblique;
char *program_name;
void usage();
void help();
void my_abort(char message[]);
void gen_outname(char outfilename[],int outchars,
char**filenames, int nfiles);
int tiff_big_endian(FILE *); /* Determines tiff file byte order */
void read_data(char infilename[],char *head, int headn,
unsigned short *data, int datan, char *trail, int trailn);
int correct_data(int bsub_flag, int intens_flag, int dist_flag, int obliq_flag,
int dist_first, int det_size, int pedestal,
unsigned short *bdat, unsigned short *indat,
float *intcorr, float *dxcorr, float *dycorr, float *tempbuf,
struct oblparams oblique);
void save_corrected(char outname[], char *head, int headn,
unsigned short *data, int datan, char *trail, int trailn);
/* ------------------------ BEGIN -------------------- */
printf("MedOptics CCD file correction program\n");
printf("Sandor Barna (Princeton), Marian Szebenyi (CHESS), ");
printf("Arthur Woll (CHESS)\n");
/* ----------------- DETECTOR PARAMETER INITIALIZATION ----------------- */
/* Check default location for correction files */
corrfiles = getenv ("CCDCORR_FILES");
if (corrfiles != NULL)
strcpy (cfilesdir,corrfiles);
else
strcpy (cfilesdir,"/usr/local/lib/");
/* MEDOPTICS PARAMETERS */
det_size = 1024;
header_bytes = 8;
trailer_bytes = 190;
swizdetect_flag = 1;
tiff_format = 1; /* If 0, swizzle must be chosen here. PC byte
order is the default, and corresponds to
swizzle = 0. */
/* dist_first = 0; if dist_first = 1; then use MOI12B1.TIF or MOI18B1.TIF */
strcpy (suffix,".tif");
oblique.wdpix = 0.00472;
oblique.hgtpix = 0.00472;
oblique.satn_level = 16384;
saturated = 16384;
strcpy (phos_name,"current.phs");
/* Taken directly the Correct package directory. -Adjust? */
strcpy (dxfilename,cfilesdir);
strcat (dxfilename,"medxbin1.tif");
strcpy (dyfilename,cfilesdir);
strcat (dyfilename,"medybin1.tif");
strcpy (intfilename,cfilesdir);
strcat (intfilename,"medi12b1.tif");
/*------------------------- Process Args --------------------*/
if (DEBUG) printf("Process Args\n");
program_name=argv[0];
if (argc<2){
printf("Too few input parameters. Type %s --help for more extensive help.\n",
program_name);
usage();
}
/* Process First Argument */
for (i=0; ((i<2) && (i<strlen(argv[1]))) ;i++){
switch (argv[1][i]) {
case 'i':
dist_first = dist_flag;
intens_flag = 1;
break;
case 'd':
dist_flag = 1;
break;
case '-': /* Dezinger and/or background subtract only */
if (i==0) {
goto process_input_files; /* Skip arg advance. Style Choice? */
/* ++argc; --argv; Counteract arg advance occuring after loop
break; */
}
else {
printf("Invalid first argument\n");
usage();
}
default:
printf("Invalid first argument\n");
usage();
}
}
--argc; ++argv;
/* Process all other arguments */
process_input_files:
while( (argc>1) && (argv[1][0] == '-')) {
switch (argv[1][1]) {
case 't':
explicit_int = 1;
strcpy(intfilename,argv[2]);
--argc; ++argv;
break;
case 'x':
explicit_dist = 1;
strcpy(dxfilename,argv[2]);
strcpy(dyfilename, dxfilename);
x_addr = strrchr(dyfilename,'x'); /* same as strrchr */
/* Address of last occurance of 'x' in the filename */
if (x_addr != NULL)
*x_addr = 'y';
else
my_abort("Invalid distortion filename -- must contain an \'x\'\n");
--argc; ++argv;
break;
case 'o':
explicit_outname = 1;
strcpy(outname,argv[2]); /* outname=argv[2]; */
--argc; ++argv;
break;
case 'b':
/* --- bfilen was initizialized to 0 --- */
bfilenames=&argv[2];
bsub_flag = 1;
while ((argc >2) && (argv[2][0] != '-')) {
if (bfilen>1) {
printf("Error,max 2 back files allowed\n");
usage();
}
/* strcpy(bfilenames[bfilen],argv[2]); */
bfilen++; /* On exit, this is the number of background
files. This should never be modified except
here. */
--argc; ++argv;
}
if (bfilen>1) dezing_back = 1;
break;
case 's':
if (argv[1][2] == 'z') dezing_images = 1;
/* Combine all images by dezingering (max 2 for the moment) */
/* infilen was initialized to 0 */
infilenames = &argv[2]; /* infilenames points to part of argv */
while ((argc >2) && (argv[2][0] != '-')) {
if ((dezing_images) && (infilen>1))
my_abort("Error, max 2 input files if combining & dezingering\n");
else if (infilen>99)
my_abort("Error, max 100 input files when not combining\n");
/* strcpy(infilenames[infilen],argv[2]); */
infile = fopen(infilenames[infilen], "rb");
if (infile==NULL) {
printf("Error opening input file %s\n",infilenames[infilen]);
my_abort("");
}
else if ((infilen==0) && swizdetect_flag)
swizzle=tiff_big_endian(infile);
fclose(infile);
infilen++; /* infilen should never be modified except
here, since it keeps track of the number of
infiles that infilenames points to */
--argc; ++argv;
}
break;
case '-':
if (strcmp(argv[1],help_flag)==0)
help();
default:
printf(" %s: Bad Option, try --help for more options\n", argv[1][0]);
usage();
}
--argc; ++argv;
}
/* -------------------- End Arg Processing -------------------- */
/* ---- If appropriate, read correction file names from file ----- */
config=fopen(cffilename,"r");
if ( (config!=NULL) && (!explicit_int) && (!explicit_dist) ) {
/* Marian's code also looked for the config file in the cffilesdir default
directory. Command line over-rides the contents of this file */
printf("Looking in %s for files\n", cffilename);
for (i=0; i<3; i++) {
fscanf(config,"%s",hold_file[i]);
}
if (!explicit_dist) {
strcpy (dxfilename,hold_file[0]);
strcpy (dyfilename,hold_file[1]);
}
if (!explicit_int)
strcpy (intfilename,hold_file[2]);
fclose(config);
}
/* ----------------------------------------------------------------------
DONE WITH FRONT MATTER. STATUS:
intfilename, dxfilename, dyfilename defined, not tested or opened
(if necessary).
bfilenames (<=2) and infilenames (<=100) are defined and tested, but
bfile and infile are closed. (Will probably move test to later).
infilen and bfilen hold the number of input and background files
dezing_images BOOLEAN whether to loop through input or
dezinger a pair then correct
dezing_back BOOLEAN indicating whether background needs dezingering
*/
/* -------------------- Begin Main Program --------------------*/
/* -------------------- INITIAL MEMORY CHECK -------------------- */
if (DEBUG) printf("Initial Memory Check\n");
length=det_size*det_size;
header=(char *)malloc(header_bytes);
if (trailer_bytes > 0)
trailer = (char *)malloc(trailer_bytes);
if (dist_flag || intens_flag || dezing_images || dezing_back){
tempbuf=(float *)malloc(length*4);
if (tempbuf == NULL)
my_abort("No memory for tempbuf -- Try upgrading from your Commodore 64!\n");
}
if ((dezing_images) || (dezing_back)) {
/* Dynamically allocate a 2D array of shorts. (see love_C.ps --
ANSI C for Programmers, Tim Love, p. 57) Currently, the limit
is 2, but we want to be able to dezinger more than one at a
time */
max_images = (infilen>bfilen?infilen:bfilen);
if (DEBUG) printf("About to allocate image_stack for %d images\n",
max_images);
image_stack=(unsigned short **)malloc(max_images*sizeof(unsigned short *));
if (image_stack == NULL)
my_abort("Out of memory when trying to allocate pointers to image_stack\n");
image_stack[0] = (unsigned short *)malloc(max_images * length * sizeof(short));
if (image_stack[0] == NULL)
my_abort("Out of memory when trying to allocate image_stack\n");
for (i=1; i < max_images; i++)
image_stack[i] = image_stack[0] + i * length;
}
/* -------------- BACKGROUND INIT --------------------- */
if (DEBUG) printf("Background Init\n");
if (bsub_flag) {
if (!(dezing_back && !dezing_images)){ /* Case 4 of TRUTH TABLE below */
bdat=(unsigned short *)malloc(length*2);
if (bdat == NULL) {
my_abort("No memory for bdat -- Try upgrading you Commodore 64!\n");
}
}
for(i=0;i<bfilen;i++){
bfile = fopen(bfilenames[i],"rb");
if (bfile==NULL) {
printf("Error opening background file %s\n",bfilenames[i]);
my_abort("");
}
dat_ptr = dezing_back ? image_stack[i] : bdat;
/* If only one background file (not dezingering),don't use image_stack */
if ((fread(header,1,header_bytes,bfile)!=header_bytes) ||
(fread(dat_ptr,2,length,bfile)!=length))
my_abort("Error reading background file!\n");
else fclose(bfile);
}
if (dezing_back){
dezing(image_stack[0],image_stack[1],det_size,det_size,saturated,tempbuf);
if (dezing_images) { /* In this case, free up image_stack to use again */
for(i=0;i<length;i++)
bdat[i]=image_stack[0][i];
}
else
bdat = image_stack[0]; /* For explanation, see comments & TRUTH TABLE
at DATA CORRECTION CODE */
}
/* shortbg=bdat; PURPOSE UNSURE */
/* Result is put into bdat -- which is also the first image */
/* Gil's dizinger code can take more than two images... */
} /* ..... if bsub_flag */
/* ----- END BACKGROUND INIT. if bsub_flag, then bdat[] is defined ------- */
/* ----- CORRECTION FILE CHECK & INIT --------- */
if (DEBUG) printf("Correction File Init\n");
if (intens_flag){
intens=fopen(intfilename,"rb");
if (intens==NULL) {
printf("Error opening intensity correction file\n\t %s\n",intfilename);
usage();
}
intcorr=(float *)malloc(length*4);
if (intcorr == NULL)
my_abort("No memory for intcorr - quit \n");
if ( (fseek(intens,TV6_HEADER_BYTES,0)) || /* fseek returns 0 on success!? */
(fread(intcorr,4,length,intens)!=length)) {
printf("Error reading Intensity Correction file %s\n", intfilename);
my_abort("");
}
fclose(intens);
/* Uncomment below (and updated makefile) to re-insert Obliquity
correct intensity correction. */
/* if ( (obliq_flag) &&
(obliq(intcorr,oblique,det_size)) )
my_abort("Error performing obliquity correction!\n"); */
}
if (dist_flag){
distx=fopen(dxfilename,"rb");
if (distx==NULL) {
printf("Error opening x distortion file %s\n",dxfilename);
usage();
}
disty=fopen(dyfilename,"rb");
if (disty==NULL) {
printf("Error opening y distortion file %s\n",dyfilename);
usage();
}
dxcorr=(float *)malloc(length*4);
dycorr=(float *)malloc(length*4);
if (!dxcorr || !dycorr)
my_abort("No memory for dxcorr & dycorr - quit \n");
if ( (fseek(distx,TV6_HEADER_BYTES,0)) || /* fseek returns 0 on success!? */
(fread(dxcorr,4,length,distx)!=length)) {
printf("Error reading X-Dist file %s!\n", dxfilename);
my_abort(""); /* -- goto qout -- */
}
if ( (fseek(disty,TV6_HEADER_BYTES,0)) || /* fseek returns 0 on success!? */
(fread(dycorr,4,length,disty)!=length)) {
printf("Error reading Y-Dist file %s!\n", dyfilename);
my_abort("");
}
fclose(distx);
fclose(disty);
}
/* If doing distortion and intensity correction, save the filenames to
configuration file */
if (dist_flag && intens_flag) {
strcpy (hold_file[0],dxfilename);
strcpy (hold_file[1],dyfilename);
strcpy (hold_file[2],intfilename);
printf("Writing correction filenames to %s\n", cffilename);
config=fopen(cffilename,"w+");
if (config!=NULL) {
fprintf(config,"%s\n%s\n%s\n",
hold_file[0],hold_file[1],hold_file[2]);
}
fclose(config);
}
/* ----- END CORRECTION INIT, intcorr, dxcorr, dycorr defined if necessary --- */
/* ----- BEGIN CORRECTION -----
dezing_images? Yes -> load & dezing the images, then
correct in proper order, then save
No -> Begin loop to load images, correct, save
Probably define saparate routines correct_dist and correct_int (ala matlab)
*/
if (DEBUG) printf("Swizzle, oblique check\n");
if (swizzle) {
if (bsub_flag)
swizit((char *)bdat,2,length);
if (dist_flag){
swizit((char *)dxcorr,4,length);
swizit((char *)dycorr,4,length);
}
if (intens_flag)
swizit((char *)intcorr,4,length);
}
/* if (pedestal == 0)
shift(bdat,length); Divide background by two to match data!*/
/* ---------> BEGIN NEW DATA CORRECTION CODE <-------------- */
/* If dezingering, it isn't really necessary to allocate separate
space for indat, since we can re-use image_stack[0] (which is
defined when dezingering either images or background. Here is the
following TRUTH TABLE to minimize memory usage:
------------------------------------------------------------
| Boolean Value | Needs own allocation |
| dezing_back |dezing_images |image_stack | indat | bdat |
| 1 | 0 | Yes | Yes | No |
| 0 | 1 | Yes | No | Yes |
| 0 | 0 | No | Yes | Yes |
| 1 | 1 | Yes | No | Yes |
------------------------------------------------------------
When not given their own allocation, bdat and/or indat just use
image_stack[0]. These optimizations rely on:
1) the background being processed before the image
2) dezing doesn't ask for a separate destination array. This is
not true of the multiple image dezinger code - see COMMENTS below
*/
if (DEBUG) printf("Correction begin\n");
if (dezing_images){
if (DEBUG) printf("Dezingering input images\n");
for(i=0;i<infilen;i++){
dat_ptr = image_stack[i];
read_data(infilenames[i], header, header_bytes, dat_ptr, length,
trailer, trailer_bytes);
if (swizzle) swizit((char *)dat_ptr,2,length);
}
dezing(image_stack[0],image_stack[1],det_size,det_size,saturated,tempbuf);
if (DEBUG) printf("Back grom dezing\n");
indat=image_stack[0];
err = correct_data(bsub_flag, intens_flag, dist_flag, obliq_flag,
dist_first, det_size, pedestal,
bdat, indat,
intcorr, dxcorr, dycorr, tempbuf,
oblique);
if (DEBUG) printf("Back grom correct_data\n");
if (err) my_abort("Error performing distortion correction!\n");
outdat=indat;
if (!explicit_outname)
gen_outname(outname, sizeof(outname),infilenames, infilen);
/* Now convert back from float to short integer format */
if ((dist_flag) || (intens_flag))
shrink(tempbuf,outdat,length,pedestal);
/* If necessary, re-swap byte order to that of the input file */
if (swizzle) swizit((char *)outdat,2,length);
save_corrected(outname, header, header_bytes, outdat, length,
trailer, trailer_bytes);
}
else { /* dezing_images == 0 */
if (DEBUG) printf("Not dezingering input images\n");
indat=(unsigned short *)malloc(length*2);
if (indat == NULL)
my_abort("No memory for indat --quit\n");
for(i=0;i<infilen;i++){
read_data(infilenames[i], header, header_bytes, indat, length,
trailer, trailer_bytes);
if (swizzle) swizit((char *)indat,2,length);
err = correct_data(bsub_flag, intens_flag, dist_flag, obliq_flag,
dist_first, det_size, pedestal,
bdat, indat,
intcorr, dxcorr, dycorr, tempbuf,
oblique);
if (err) goto skipthis;
outdat=indat;
if (!explicit_outname){
if (DEBUG) printf("Generating Output filename\n");
gen_outname(outname, sizeof(outname),&infilenames[i], 1);
/* Need &infilenames[i] == infilenames + i, since
gen_outname expects a pointer to a string, not a
string */
}
if ((dist_flag) || (intens_flag)){
if (DEBUG)
printf("Now convert back from float to short integer format\n");
shrink(tempbuf,outdat,length,pedestal);
}
/* If necessary, re-swap byte order to that of the input file */
if (swizzle) {
if (DEBUG) printf("Unswizzling data\n");
swizit((char *)outdat,2,length);
}
save_corrected(outname, header, header_bytes, outdat, length,
trailer, trailer_bytes);
skipthis:
if (err) {
strcpy(response,"Y");
editin("Continue?",response,2);
if (response[0] != 'Y' && response[0] !='y')
my_abort("");
}
err=0;
if (DEBUG) printf("End (dezing_image == 0) branch, i = %d\n", i);
} /* End for loop within dezing_image == 0 branch */
} /* ... else .... */
if (DEBUG) printf("End main()\n");
/* ---------> END NEW DATA CORRECTION CODE <-------------- */
}
int correct_data(int bsub_flag, int intens_flag, int dist_flag, int obliq_flag,
int dist_first, int det_size, int pedestal,
unsigned short *bdat, unsigned short *indat,
float *intcorr, float *dxcorr, float *dycorr, float *tempbuf,
struct oblparams oblique)
{
void intenscorr(float *buffer,float *intcorr,int size,
struct oblparams *params);
int distcorr(unsigned short *indat,float *tempbuf,float *delx,float *dely,
int size, int pedestal);
void my_abort(char message[]);
if (bsub_flag){
printf("Subtracting background.\n");
bsub(indat,bdat,det_size,pedestal);
}
if (dist_first){
printf("Performing distortion correction.\n");
if (distcorr(indat,tempbuf,dxcorr,dycorr,det_size,pedestal))
return(1);
/* returns 0 on success */
printf("Performing intensity correction.\n");
intenscorr_int(indat,intcorr,det_size,&oblique);
}
else if (intens_flag){
printf("Performing intensity correction.\n");
intenscorr_int(indat,intcorr,det_size,&oblique);
}
if (dist_flag){
printf("Performing distortion correction.\n");
if (distcorr(indat,tempbuf,dxcorr,dycorr,det_size,pedestal))
return(1);
}
return(0);
}
void read_data(char infilename[],char *head, int headn,
unsigned short *data, int datan, char *trail, int trailn)
{
FILE *infile = NULL;
void my_abort(char message[]);
infile = fopen(infilename,"rb");
if (infile==NULL) {
printf("Error opening input data file %s\n",infilename);
my_abort("");
}
if ( (fread(head,1,headn,infile)!=headn) ||
(fread(data,2,datan,infile)!=datan) ) {
printf("Error reading input data file %s\n",infilename);
my_abort("");
}
if (trailn > 0) {
if(fread(trail,1,trailn,infile)!=trailn){
printf("Error reading trailer bytes from %s\n",infilename);
my_abort("");
}
}
fclose(infile);
}
void save_corrected(char outname[], char *head, int headn,
unsigned short *data, int datan, char *trail, int trailn)
{
FILE *outfile = NULL;
void my_abort(char message[]);
outfile = fopen(outname,"r");
if (outfile !=NULL){ /* IF this file exists, close file and abort */
fclose(outfile);
printf("Error: Output file %s already exists\n", outname);
my_abort("");
}
if((outfile=fopen(outname,"wb"))==NULL) {
printf("Error opening %s\n",outname);
my_abort("");
}
printf("Now writing data to output file %s\n",outname);
if ( (fwrite(head,1,headn,outfile)!=headn) ||
(fwrite(data,2,datan,outfile)!=datan)) {
printf("Error writing to %s\n",outname);
my_abort("");
}
/* printf("Now writing trailer bytes to output file %s\n",outname); */
if (trailn > 0) {
if (fwrite(trail,1,trailn,outfile)!=trailn) {
printf("Error writing trailer bytes to to %s\n",outname);
my_abort("");
}
}
fclose(outfile);
}
void gen_outname(char outfilename[], int outchars,char **filenames,
int nfilenames)
/* Takes an array of number filenames (infiles[][]), each no more
than outsize long, and generates another string based on an
algorithm. Uses a different algorithm depending on how many
filenames are given
*/
{
extern const char insert[];
extern const char dot;
extern const char default_outname[];
char *insert_addr;
void str_ins(char [], int, char [], const char [], const char);
if (nfilenames == 1)
str_ins(outfilename,outchars,filenames[0],insert,dot);
else
outfilename=strncpy(outfilename,default_outname,outchars);
return;
}
void str_ins(char out[], int outsize,char in[],
const char ins[], const char tok)
/* str_ins inserts the string ins[] just in front of the right-most
instance of tok in string in[], the result being placed in out[].
if tok is not found in in[], in and ins are concatenated */
{
void my_abort(char message[]);
char *ins_addr; /* Where to insert new sting ins []*/
int ins_ind, last_ind, i, slide;
strcpy(out,in);
if (ins_addr=(strrchr(out,tok)))
ins_ind=(int)(ins_addr-out);
else{ /* tok not found in out[], put ins[] at the end */
out=strcat(out,ins);
return;
}
last_ind = strlen(out); /* out[last_ind] is the terminating null char */
slide = strlen(ins);
if (last_ind+slide > outsize-1)
my_abort("in str_ins, out[] not big enough to insert ins[]\n");
for(i=last_ind;i>=ins_ind;i--)
out[i+slide]=out[i];
for(i=ins_ind ; i < (ins_ind + slide) ; i++)
out[i]=ins[i-ins_ind];
return;
/* NOTE: for the above, might be better to concatenate and THEN
move, rather than doing it as I've done here. */
}
void usage()
{
extern const char usage_str[];
printf(usage_str);
exit(0);
}
void help()
{
extern const char help_str[];
printf(help_str);
exit(0);
}
void my_abort(char message[]){
printf(message);
printf("Aborting\n");
exit(0);
}
int tiff_big_endian(FILE *datfile)
/* datfile is an open binary file. Reads the first 4 bytes and,
assuming that they correspond to the Tiff 'magic number' (I
think that's right), determines byte order.
ASSUMES type int has 4 bytes.
PC byte order (default) = little endian
SGI byte order = big endian
*/
{
int magic;
int big_endian = 0;
rewind(datfile);
if (fread(&magic,sizeof(int),1,datfile)!=1)
my_abort("Error reading data file (in tiff_big_endian)\n");
if (magic == 0x002A4949)
big_endian=0; /* PC byte order - little endian */
else if (magic ==0x49492A00) { /* SGI byte order - big endian */
big_endian=1;
printf("Reversing byte order of input files for the calculations.\n");
printf("They will be output in the original format.\n");
}
else
my_abort("Image file doens't appear to be a tiff file after all\n");
return big_endian;
}
/*
COMMENTS
DEZINGERING:
Currently, use dezinger.c declared as:
void dezing(unsigned short \*indat1,
unsigned short \*indat2,
int slow_max,
int fast_max,
int saturated,
float \*tempbuf)
The current code anticipates a switch to the following, which
dezingers >2 files, by implementing the 2D array **image_stack,
required to make this switch
(From Zingers.c on Gil Toombes\' web site -- Gruner group home page)
int remove_zingers(int num_images,
uint16 \*\*data_ptrs,
double \*dest,
int width,
int height,
long bkg,
double sigma,
double sigma_cut,
double gain,
long int *num_zingers)
BITSHIFTING:
Bitshifting is NOT done. Instead, if the a background pixel has more counts than the data pixel, bsub (in medbsub) just sets the pixel to zero
Bitshifting is necessary for any signal that can exceed 2^15 in order to do
background subtraction. The MedOptics detector\'s ADC is only 14
bits and so constitutes a special case in which the output never
exceeds 2^14=16,384.
PEDESTAL: Not used at present -- may have to re-insert.
SATURATION: Any saturated pixel (detector specific) will
maintain its value throughout the correction
procedure.
OBLIQUITY: Disabled, but flag obliq_flag is defined in case we want to
replace. Current version uses prompts (within module GetObl) to
determine parameters.
BYTE ORDER: For tiff files, program detects whether or not swizzling
of the byte order is needed. If so, it is automatically performed
by module \'swizit\': swizit: Swizzles the byte order (used for
SGI computers) for n-byte variables.
*/