1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
|
/****************************************************************************
*
* Synopsis: cf_geometric_distort(fitsfile *header, long nevents,
* float *xfarf, float *yfarf,
* unsigned char *locflags)
*
* Description: Corrects the X and Y positions of the event for geometric
* distortion.
*
* Arguments: fitsfile *header Pointer to the location of the FITS
* header of the Intermediate Data File
* long nevents number of photons in the file
* float *xfarf X positions for each photon event
* float *yfarf Y positions for each photon event
* unsigned char *locflags Location flags for each event
*
* Calls : None
*
* Return: 0 on success
*
* History: 08/09/02 1.1 rdr Begin work
* 11/12/02 1.2 peb Added check to move only events in
* active region. Added locflags to
* function cal.
* 03/11/03 1.3 wvd Changed locflags to unsigned char.
* 05/20/03 1.4 rdr Added call to cf_proc_check
* 07/29/03 1.5 wvd If cf_proc_check fails, return errflg.
* 03/26/04 1.6 rdr Modify specbiny
* 04/05/04 1.7 wvd Modify specbiny only if specbiny > 1.
* 04/07/07 1.8 wvd Clean up compiler warnings.
*
****************************************************************************/
#include <string.h>
#include <stdlib.h>
#include "calfuse.h"
int
cf_geometric_distort(fitsfile *header, long nevents, float *xfarf,
float *yfarf, unsigned char *locflags)
{
char CF_PRGM_ID[] = "cf_geometric_distort";
char CF_VER_NUM[] = "1.8";
char geomfile[FLEN_VALUE]={'\0'};
short *geomx_img, *geomy_img;
int errflg=0, status=0, hdutype, anynull, nullval;
int specbiny ;
long fpixel=1, j;
long xdim_gx, ydim_gx, xdim_gy, ydim_gy, npix_gx, npix_gy;
float scale_gx=1., zero_gx=0., scale_gy=1., zero_gy=0.;
float yexp ;
fitsfile *geomfits=NULL;
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
if ((errflg = cf_proc_check(header, CF_PRGM_ID))) return errflg;
/*
* Get the name of the calibration file and open it
*/
FITS_read_key(header, TINT, "SPECBINY", &specbiny, NULL, &status);
cf_verbose(5,"Initial value of specbiny=%d",specbiny) ;
FITS_read_key(header, TSTRING, "GEOM_CAL", geomfile, NULL, &status);
cf_verbose(4, "Geometric correction file = %s", geomfile);
FITS_open_file(&geomfits, cf_cal_file(geomfile), READONLY, &status);
/* Read the y expansion factor from the calibration header */
FITS_movabs_hdu(geomfits, 1, &hdutype, &status);
FITS_read_key(geomfits, TFLOAT, "YEXPAND", &yexp, NULL, &status);
cf_verbose(3,"Y scale expansion factor = %f",yexp) ;
/* If SPECBINY > 1, scale by Y expansion factor */
if (specbiny > 1) {
specbiny = (int) (0.5 + (float) specbiny * yexp) ;
FITS_update_key(header, TINT, "SPECBINY", &specbiny, NULL, &status) ;
cf_verbose(5,"Corrected value of specbiny = %d",specbiny) ;
}
/*
* Read in the X distortion calibration image
*/
FITS_movabs_hdu(geomfits, 2, &hdutype, &status);
FITS_read_key(geomfits, TLONG, "NAXIS1", &xdim_gx, NULL, &status);
FITS_read_key(geomfits, TLONG, "NAXIS2", &ydim_gx, NULL, &status);
npix_gx = xdim_gx * ydim_gx;
FITS_read_key(geomfits, TFLOAT, "BSCALE", &scale_gx, NULL, &status);
FITS_read_key(geomfits, TFLOAT, "BZERO", &zero_gx, NULL, &status);
geomx_img = (short *) cf_malloc(sizeof(short) * npix_gx);
fits_set_bscale(geomfits, 1., 0., &status);
FITS_read_img(geomfits, TSHORT, fpixel, npix_gx, &nullval, geomx_img,
&anynull, &status);
/*
* Read the Y distortion calibration image
*/
FITS_movabs_hdu(geomfits, 3, &hdutype, &status);
FITS_read_key(geomfits, TLONG, "NAXIS1", &xdim_gy, NULL, &status);
FITS_read_key(geomfits, TLONG, "NAXIS2", &ydim_gy, NULL, &status);
npix_gy = xdim_gy * ydim_gy;
FITS_read_key(geomfits, TFLOAT, "BSCALE", &scale_gy, NULL, &status);
FITS_read_key(geomfits, TFLOAT, "BZERO", &zero_gy, NULL, &status);
geomy_img = (short *) cf_malloc(sizeof(short) * npix_gy);
fits_set_bscale(geomfits, 1., 0., &status);
FITS_read_img(geomfits, TSHORT, fpixel, npix_gy, &nullval, geomy_img,
&anynull, &status);
FITS_close_file(geomfits, &status);
/*
* Make sure that the X and Y distortion images are the same size
*/
if (xdim_gx != xdim_gy || ydim_gx != ydim_gy)
cf_if_warning("X and Y distortion images not the same size.");
/*
* Use the data in the geometric distortion calibration images to
* correct XFARF and YFARF
*/
for (j=0; j<nevents; j++) {
/*
* Move only events in active region.
*/
if (!(locflags[j] & LOCATION_SHLD)) {
int xndx = (int) xfarf[j], yndx = (int) yfarf[j], ndx;
/*
* Check that the y indices are within bounds
*/
if (yndx < 0)
yndx = 0;
else if (yndx > ydim_gy-1)
yndx = ydim_gy-1;
ndx = yndx*xdim_gx + xndx;
xfarf[j] += geomx_img[ndx]*scale_gx + zero_gx;
yfarf[j] += geomy_img[ndx]*scale_gy + zero_gy;
}
}
free(geomx_img);
free(geomy_img);
cf_proc_update(header, CF_PRGM_ID, "COMPLETE");
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done Processing");
return 0;
}
|
