Skip to content

Commit d2ef21f

Browse files
mdlpstscimdlpstsci
committed
Updated routines to accommdate support subarray processing using the new bias level processing which previously only applied to full frame data.
1 parent 68c27eb commit d2ef21f

File tree

6 files changed

+764
-816
lines changed

6 files changed

+764
-816
lines changed

‎pkg/acs/calacs/acsccd/blev_funcs_postsm4.c

Lines changed: 200 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,11 @@
22
#include <stdlib.h>
33
#include <string.h>
44
#include <stdio.h>
5+
#include <stdarg.h> /* Supports the variable number of parameters for bias_shift_corr() */
6+
#include <ctype.h>
7+
#include <assert.h>
58

9+
#include "hstcal_memory.h"
610
#include "hstcal.h"
711
#include "hstio.h"
812

@@ -36,152 +40,243 @@ static int make_amp_array(const int arr_rows, const int arr_cols, SingleGroup *
3640
static int unmake_amp_array(const int arr_rows, const int arr_cols, SingleGroup * im,
3741
int amp, double * array);
3842

43+
/* Remove the signal dependent bias shift from post-SM4 full frame WFC data
44+
* (based on ISR http://www.stsci.edu/hst/acs/documents/isrs/isr1202.pdf),
45+
* as well as from supported subarray data (~post-mid-2016) as identified by
46+
* specific aperture names
47+
* (based on ISR http://www.stsci.edu/hst/acs/documents/isrs/isr1703.pdf).
48+
*
49+
* chip2 consists of amps C & D, chip1 consists of amps A & B.
50+
*
51+
* Because full frame data uses both chips, but supported subarray data
52+
* is only from a single amp (therefore, a portion of a single chip), this
53+
* routine has a variable number of input SingleGroups.
54+
*
55+
* *** This routine should be modified according to issues discussed in PR #312 and noted in IT #334. ***
56+
*
57+
* 16-May-2018 M.D. De La Pena: Generalized bias_shift_corr() to handle subarray data.
58+
*/
59+
int bias_shift_corr(ACSInfo *acs, int nGroups, ...) {
60+
61+
/* This function has built in assumptions regarding the data being passed */
62+
assert (acs != NULL); /* Make sure the ACSInfo structure is defined */
63+
assert ((nGroups == 1) || (nGroups == 2)); /* Must be at least one and no more than two SingleGroups */
64+
assert (((nGroups == 1) && (acs->subarray == YES)) || ((nGroups == 2) && (acs->subarray == NO))); /* Subarray has one and full frame has two SingleGroups */
65+
66+
/* Code to handle variable number of input parameters. The number of SingleGroup is 2 when
67+
processing a full frame (2 chips - all 4 amps) and only 1 when processing a supported subarray
68+
(1 chip from a single amp).
69+
*/
70+
SingleGroup *sg[nGroups];
71+
va_list arguments;
72+
73+
/* Initializing arguments to store all the values being passed to this routine via "...".
74+
* The nGroups value and the number of SingleGroups actually passed must match.
75+
*/
76+
va_start (arguments, nGroups);
77+
{
78+
/* sg[0] is chip2 and sg[1] is chip 1 for a full frame image */
79+
{unsigned int i;
80+
for (i = 0; i < nGroups; i++) {
81+
sg[i] = NULL;
82+
sg[i] = va_arg (arguments, SingleGroup *);
83+
if (!sg[i]) {
84+
return (status = INTERNAL_ERROR);
85+
}
86+
}}
87+
}
88+
va_end (arguments);
3989

40-
/* remove signal dependent bias shift from post-SM4 full frame WFC data.
41-
* based on ISR http://www.stsci.edu/hst/acs/documents/isrs/isr1202.pdf
42-
* chip2 is amps C & D, chip1 is amps A & B. */
43-
int bias_shift_corr(ACSInfo *acs, SingleGroup *chip2, SingleGroup *chip1) {
4490
extern int status;
45-
46-
int i, j, k; /* iteration variables */
47-
48-
/* array for amp data and amp data + gap pixels */
49-
double * ampdata, * ampdata_gap;
50-
91+
5192
const double serial_freq = 1000./22.; /* serial pixel frequency */
5293
const double parallel_shift = 3.212; /* parallel shift time */
5394

5495
/* number of virtual pixels at end of each row */
5596
const int ngap_pix = (int) (serial_freq * parallel_shift + 0.5);
5697

57-
const int arr_rows = chip2->sci.data.ny;
58-
const int arr_cols = chip2->sci.data.nx/2;
98+
/* There is always at least one SingleGroup - use sg[0] to get the image size of a single amp.
99+
* If this is a full frame image, then "nx" represents all the columns across both amps plus
100+
* the physical overscan and must be divided by 2. If this is a subarray, then "nx" represents
101+
* just the columns for a single amp plus physical overscan.
102+
*/
103+
const int arr_rows = sg[0]->sci.data.ny;
104+
const int arr_cols = (acs->subarray == NO) ? sg[0]->sci.data.nx/2 : sg[0]->sci.data.nx;
59105

60106
/* total number of real and gap pixels per quadrant */
61107
const int nquad_pix = (arr_cols + ngap_pix) * arr_rows;
62108

63-
/* array of true DC bias levels */
64-
double * dc_bias_levels;
65-
66109
/* arrays below are in order of amp A, B, C, D */
67110

68-
/* time constant of AC high pass filter in external pre-amp. */
111+
/* time constant of AC high pass filter in external pre-amp in microseconds */
69112
const double time_const[NAMPS] = {37.290251, 36.180001, 37.867770, 42.461249};
70113

71114
/* ratio of DC offset shift and pixel signal */
72115
const double dc_ratio[NAMPS] = {0.3, 0.3, 0.3, 0.3};
73116

74117
/* DSI sensitivity */
75118
const double dsi_sens[NAMPS] = {2.9188919e-3, 10.805754e-3, 12.432145e-3, 3.8596253e-3};
119+
120+
PtrRegister ptrReg;
121+
initPtrRegister(&ptrReg);
122+
123+
/* array of true DC bias levels */
124+
double * dc_bias_levels = NULL;
125+
dc_bias_levels = malloc((nquad_pix + 1) * sizeof(double));
126+
if (!dc_bias_levels) {
127+
freeOnExit (&ptrReg);
128+
return (status = OUT_OF_MEMORY);
129+
}
130+
addPtr(&ptrReg, dc_bias_levels, &free);
76131

77-
/* factor combining time constant and clocking frequency */
78-
double factor;
79-
80-
/* summation variables */
81-
double sum;
82-
int num;
83-
double magic_square_mean;
84-
85-
/* allocate space for data arrays */
86-
ampdata = malloc(arr_rows * arr_cols * sizeof(double));
87-
ampdata_gap = malloc(nquad_pix * sizeof(double));
132+
/* array for amp data and amp data + gap pixels */
133+
double *ampdata = malloc(arr_rows * arr_cols * sizeof(double));
134+
if (!ampdata) {
135+
freeOnExit (&ptrReg);
136+
return (status = OUT_OF_MEMORY);
137+
}
138+
addPtr(&ptrReg, ampdata, &free);
139+
double *ampdata_gap = malloc(nquad_pix * sizeof(double));
140+
if (!ampdata_gap) {
141+
freeOnExit (&ptrReg);
142+
return (status = OUT_OF_MEMORY);
143+
}
144+
addPtr(&ptrReg, ampdata_gap, &free);
88145

89-
/* allocate space for true DC bias levels array */
90-
dc_bias_levels = malloc((nquad_pix + 1) * sizeof(double));
146+
/* For full frame data, there are four amps in use. For supported subarray
147+
* data, there is only one amp in use. The actual amp in use must be determined.
148+
*/
149+
int numAmpsInUse = (acs->subarray == YES) ? 1 : NAMPS;
91150

92-
for (i = 0; i < NAMPS; i++) {
93-
/* put a single amp's data into ampdata */
94-
if (i < 2) {
95-
make_amp_array(arr_rows, arr_cols, chip1, i, ampdata);
96-
} else {
97-
make_amp_array(arr_rows, arr_cols, chip2, i, ampdata);
98-
}
151+
/* Loop over all the amps for full frame data (total of 4 amps, 0 to 3) or just use
152+
* the single amp value corresponding to the supported subarray (1 amp, value defined in acs.h).
153+
*/
154+
int ampInUse = 0;
155+
{unsigned int i;
156+
for (i = 0; i < numAmpsInUse; i++) {
157+
ampInUse = i;
158+
if ((i == 0) && (acs->subarray == YES)) {
159+
if (strcmp (acs->ccdamp, "A") == 0)
160+
ampInUse = AMP_A;
161+
else if (strcmp (acs->ccdamp, "B") == 0)
162+
ampInUse = AMP_B;
163+
else if (strcmp (acs->ccdamp, "C") == 0)
164+
ampInUse = AMP_C;
165+
else if (strcmp (acs->ccdamp, "D") == 0)
166+
ampInUse = AMP_D;
167+
}
168+
169+
/* put a single amp's data into ampdata */
170+
if (acs->subarray == YES) {
171+
make_amp_array(arr_rows, arr_cols, sg[0], ampInUse, ampdata);
172+
} else {
173+
if (i < 2) {
174+
make_amp_array(arr_rows, arr_cols, sg[1], ampInUse, ampdata);
175+
} else {
176+
make_amp_array(arr_rows, arr_cols, sg[0], ampInUse, ampdata);
177+
}
178+
}
99179

100-
/* calculate "magic square mean" */
101-
sum = 0.0;
102-
num = 0;
180+
/* calculate "magic square mean" - this works for any single amp */
181+
double sum = 0.0;
182+
int num = 0;
183+
double magic_square_mean = 0.0;
103184

104-
for (j = 2057; j <= 2066; j++) {
105-
for (k = 13; k <= 22; k++) {
106-
sum += ampdata[arr_cols*j + k];
107-
num++;
108-
}
109-
}
185+
{unsigned int j, k;
186+
for (j = 2057; j <= 2066; j++) {
187+
for (k = 13; k <= 22; k++) {
188+
sum += ampdata[arr_cols*j + k];
189+
num++;
190+
}
191+
}}
110192

111-
magic_square_mean = sum / (double) num;
193+
magic_square_mean = sum / (double) num;
112194

113-
/* make amp + gap array */
114-
for (j = 0; j < arr_rows; j++) {
115-
for (k = 0; k < (arr_cols + ngap_pix); k++) {
116-
if (k < arr_cols) {
117-
ampdata_gap[(arr_cols + ngap_pix)*j + k] = ampdata[arr_cols*j + k];
118-
} else {
119-
ampdata_gap[(arr_cols + ngap_pix)*j + k] = magic_square_mean;
195+
/* make amp + gap array */
196+
{unsigned int j, k;
197+
for (j = 0; j < arr_rows; j++) {
198+
for (k = 0; k < (arr_cols + ngap_pix); k++) {
199+
if (k < arr_cols) {
200+
ampdata_gap[(arr_cols + ngap_pix)*j + k] = ampdata[arr_cols*j + k];
201+
} else {
202+
ampdata_gap[(arr_cols + ngap_pix)*j + k] = magic_square_mean;
203+
}
120204
}
121-
}
122-
}
205+
}}
123206

124-
/* calculate true DC bias levels */
125-
factor = 1.0 - exp(-1.0 / (time_const[i] * serial_freq));
207+
/* factor combining time constant and clocking frequency */
208+
double factor = 1.0 - exp(-1.0 / (time_const[ampInUse] * serial_freq));
126209

127-
dc_bias_levels[0] = magic_square_mean * dc_ratio[i];
210+
/* calculate true DC bias levels */
211+
dc_bias_levels[0] = magic_square_mean * dc_ratio[ampInUse];
128212

129-
for (j = 1; j < nquad_pix + 1; j++) {
130-
dc_bias_levels[j] = ampdata_gap[j-1] * factor * dc_ratio[i] +
213+
{unsigned int j;
214+
for (j = 1; j < nquad_pix + 1; j++) {
215+
dc_bias_levels[j] = ampdata_gap[j-1] * factor * dc_ratio[ampInUse] +
131216
(1.0 - factor) * dc_bias_levels[j-1];
132-
}
217+
}}
133218

134-
/* calculate correction to data */
135-
for (j = 0; j < nquad_pix; j++) {
136-
ampdata_gap[j] = (ampdata_gap[j] - dsi_sens[i] * dc_bias_levels[j+1]) -
219+
/* calculate correction to data */
220+
{unsigned int j;
221+
for (j = 0; j < nquad_pix; j++) {
222+
ampdata_gap[j] = (ampdata_gap[j] - dsi_sens[ampInUse] * dc_bias_levels[j+1]) -
137223
(10./22.) * (dc_bias_levels[j+1] - dc_bias_levels[j]);
138-
}
224+
}}
139225

140-
/* copy corrected data back to ampdata */
141-
for (j = 0; j < arr_rows; j++) {
142-
for (k = 0; k < arr_cols; k++) {
143-
ampdata[arr_cols*j + k] = ampdata_gap[(arr_cols + ngap_pix)*j + k];
144-
}
145-
}
226+
/* copy corrected data back to ampdata */
227+
{unsigned int j, k;
228+
for (j = 0; j < arr_rows; j++) {
229+
for (k = 0; k < arr_cols; k++) {
230+
ampdata[arr_cols*j + k] = ampdata_gap[(arr_cols + ngap_pix)*j + k];
231+
}
232+
}}
146233

147-
/* re-calculate "magic square mean" */
148-
sum = 0.0;
149-
num = 0;
234+
/* re-calculate "magic square mean" */
235+
sum = 0.0;
236+
num = 0;
237+
magic_square_mean = 0.0;
150238

151-
for (j = 2057; j <= 2066; j++) {
152-
for (k = 13; k <= 22; k++) {
153-
sum += ampdata[arr_cols*j + k];
154-
num++;
155-
}
156-
}
239+
{unsigned int j, k;
240+
for (j = 2057; j <= 2066; j++) {
241+
for (k = 13; k <= 22; k++) {
242+
sum += ampdata[arr_cols*j + k];
243+
num++;
244+
}
245+
}}
157246

158-
magic_square_mean = sum / (double) num;
159-
/* Keep track of the total bias correction applied to each AMP region */
160-
acs->blev[i] += magic_square_mean;
161-
/* Report to the user the contribution to the bias level correction
162-
made by this processing step.
163-
*/
164-
sprintf(MsgText, "Bias shift correcting for bias level in Amp %c of %0.4f electrons.",AMPSORDER[i], magic_square_mean);
165-
trlmessage(MsgText);
247+
magic_square_mean = sum / (double) num;
248+
249+
/* Keep track of the total bias correction applied to each AMP region */
250+
acs->blev[ampInUse] += magic_square_mean;
251+
252+
/* Report to the user the contribution to the bias level correction
253+
made by this processing step.
254+
*/
255+
sprintf(MsgText, "Bias shift correcting for bias level in Amp %c of %0.4f electrons.",AMPSORDER[ampInUse], magic_square_mean);
256+
trlmessage(MsgText);
166257

167-
/* subtract "magic square mean" from data*/
168-
for (j = 0; j < arr_rows; j++) {
169-
for (k = 0; k < arr_cols; k++) {
170-
ampdata[arr_cols*j + k] -= magic_square_mean;
171-
}
172-
}
258+
/* subtract "magic square mean" from data*/
259+
{unsigned int j, k;
260+
for (j = 0; j < arr_rows; j++) {
261+
for (k = 0; k < arr_cols; k++) {
262+
ampdata[arr_cols*j + k] -= magic_square_mean;
263+
}
264+
}}
173265

174-
/* copy modified data back to SingleGroup structs */
175-
if (i < 2) {
176-
unmake_amp_array(arr_rows, arr_cols, chip1, i, ampdata);
177-
} else {
178-
unmake_amp_array(arr_rows, arr_cols, chip2, i, ampdata);
179-
}
180-
}
266+
/* copy modified data back to SingleGroup structs */
267+
268+
if (acs->subarray == YES) {
269+
unmake_amp_array(arr_rows, arr_cols, sg[0], ampInUse, ampdata);
270+
} else {
271+
if (i < 2) {
272+
unmake_amp_array(arr_rows, arr_cols, sg[1], ampInUse, ampdata);
273+
} else {
274+
unmake_amp_array(arr_rows, arr_cols, sg[0], ampInUse, ampdata);
275+
}
276+
}
277+
}}
181278

182-
free(ampdata);
183-
free(ampdata_gap);
184-
free(dc_bias_levels);
279+
freeOnExit(&ptrReg);
185280

186281
return status;
187282
}
@@ -894,6 +989,7 @@ static int make_amp_array(const int arr_rows, const int arr_cols, SingleGroup *i
894989
array[i*arr_cols + j] = Pix(im->sci.data, c, r);
895990
}
896991
}
992+
trlmessage (MsgText);
897993

898994
return status;
899995
}

0 commit comments

Comments
 (0)