1
/*===========================================================================
2
Copyright (C) 1989-2009 European Southern Observatory (ESO)
4
This program is free software; you can redistribute it and/or
5
modify it under the terms of the GNU General Public License as
6
published by the Free Software Foundation; either version 2 of
7
the License, or (at your option) any later version.
9
This program is distributed in the hope that it will be useful,
10
but WITHOUT ANY WARRANTY; without even the implied warranty of
11
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12
GNU General Public License for more details.
14
You should have received a copy of the GNU General Public
15
License along with this program; if not, write to the Free
16
Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
19
Correspondence concerning ESO-MIDAS should be addressed as follows:
20
Internet e-mail: midas@eso.org
21
Postal address: European Southern Observatory
22
Data Management Division
23
Karl-Schwarzschild-Strasse 2
24
D 85748 Garching bei Muenchen
26
===========================================================================*/
29
/* ++++++++++++++++++++++++++++++++++++++++++++++++++
32
program AVERAG version 1.00 890120
33
K. Banse ESO - Garching
36
bulk data frames, average
39
average a no. of image frames, the result will either be a frame with size
40
equal to common area of all frames involved or size equal to union of all
44
extract input frames either from given parameter string or from catalogue,
45
add up all frames + divide in the end
48
the following keys are used:
50
ACTION/C/1/2 (1:1) = M or N, for merging option
51
(2:2) = I or W, for AVERAGE/IMAGE, /WEIGHTS
52
OUT_A/C/1/60 result frame
53
P3/C/1/80 list of frames to be averaged
58
-------------------------------------------------- */
60
#include <midas_def.h>
65
#define MAXIMSONE MAXIMS+1
75
ndim = dimension of array rka for sorting
76
but we pass the arrays with 1 element in front, so that the algorithm
77
sorts from [1] -> [ndim]
78
this is the Heapsort algorithm from "Numerical Recipes", page 231
90
m = (ndim>>1) + 1; /* ndim/2 + 1 */
94
if (m > 1) /* still hiring */
96
else /* in retirement + promotion phase */
98
zka = rka[ndim]; /* clear a space at end of array rka */
99
rka[ndim] = rka[1]; /* retire the top of the heap into it */
100
if (--ndim == 1) /* done with last promotion */
102
rka[1] = zka; /* the least competent guy of all... */
103
return; /* that's it folks */
107
/* in hiring as well as in promotion phase */
108
/* here set up to shift down zka to its right level */
111
j = m << 1; /* in FORTRAN: j = m + m */
118
if (rka[j] < rka[j+1]) j++;
121
if (zka < rka[j]) /* demote zka */
125
j = j << 1; /* in FORTRAN: j = j + j */
131
rka[i] = zka; /* put zka into its slot */
139
void fill(iaux,faux,a,x,z,apix,cpix,npixa,npixc)
140
int *iaux, apix[3][2], *cpix;
147
int nn, nin, nini, nout, jump, indx, cntdx;
148
int nx, ny, xfirst, yfirst, xend, yend;
155
/* ---------------------------
157
the intermediate z-space is structured as follows:
159
pix(1) pix(2) up to pix(FRMCNT) 1 )
160
pix(1) pix(2) up to pix(FRMCNT) 2 )
162
pix(1) pix(2) up to pix(FRMCNT) NPIXC(1) )
165
--------------------------- */
168
/* if here for the first time, init the count `pixels' */
175
nn = npixc[0] * npixc[1];
176
if ((iaux[5] == 0) && (iaux[2] == 0))
177
nini = frmcnt; /* for `nomerge' and `all data' */
181
for (nr=0; nr<nn; nr++)
185
if (iaux[0] == 0) return; /* if iaux[0] = 0, not much to do */
188
nini = 0; /* input frame begins with 1. valid pixel */
192
/* --------------------------------- */
193
/* here for the non-merging option */
194
/* --------------------------------- */
196
if (iaux[2] == 0) /* take all pixels */
199
for (nn=0; nn<npixc[1]; nn++)
202
for (nr=0; nr<npixc[0]; nr++)
207
nini += npixa; /* follow with input index */
210
else /* take only pixels in valid interval */
214
for (nn=0; nn<npixc[1]; nn++)
217
for (nr=0; nr<npixc[0]; nr++)
220
if ((rr >= faux[2]) && (rr <= faux[3]))
222
nout = indx + x[cntdx];
224
x[cntdx] ++; /* increment count */
229
nini += npixa; /* follow with input index */
237
/* ----------------------------- */
238
/* here for the merging option */
239
/* ----------------------------- */
241
nx = apix[0][1] - apix[0][0]; /* get overlapping part */
242
ny = apix[1][1] - apix[1][0];
247
jump = frmcnt * npixc[0];
252
if (iaux[2] == 0) /* take all pixels */
254
for (nn=0; nn<npixc[1]; nn++)
256
if ((nn >= yfirst) && (nn <= yend))
259
for (nr=0; nr<npixc[0]; nr++)
261
if ((nr >= xfirst) && (nr <= xend))
263
nout = indx + x[cntdx];
265
x[cntdx] ++; /* increment count */
271
nini += npixa; /* follow with input index */
275
indx += jump; /* jump a whole line; */
280
else /* take only pixels in valid interval */
282
for (nn=0; nn<npixc[1]; nn++)
284
if ((nn >= yfirst) && (nn <= yend))
287
for (nr=0; nr<npixc[0]; nr++)
289
if ((nr >= xfirst) && (nr <= xend))
292
if ((rr >= faux[2]) && (rr <= faux[3]))
294
nout = indx + x[cntdx];
296
x[cntdx] ++; /* increment count */
303
nini += npixa; /* follow with input index */
307
indx += jump; /* jump a whole line; */
319
void add(flag,iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
321
int flag, *iaux, *nc;
323
float *faux, *z, *c, usrnul, *cuts;
326
int frmcnt, n, nn, nh, k, mcc, ma, mb, indx, size;
329
register float rval, va, vb;
330
float rbuf[MAXIMSONE]; /* buffer for sorting */
331
static float old = 0.0;
340
size = npixc[0] * npixc[1];
343
/* branch according to FLAG */
345
if (flag == 4) goto median; /* that is the `worst' part */
350
/* -------------------------------*/
351
/* here for the normal average */
352
/* -------------------------------*/
355
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
362
rval = old; /* use last value */
365
mcc ++; /* increment null count */
372
for (n=indx+1; n<indx+nn; n++)
382
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
383
if (cuts[1] < rval) cuts[1] = rval;
390
/* -----------------------*/
391
/* here for the minimum */
392
/* -----------------------*/
395
if (iaux[3] > 0) /* take average of MIN and the next */
396
{ /* `iaux[3]' higher values */
397
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
404
rval = old; /* use last value */
407
mcc ++; /* increment null count */
413
k = 1; /* start at index 1 ... */
414
for (n=indx; n<indx+nn; n++)
416
sortit(rbuf,nn); /* sort array */
422
sum = rbuf[1]; /* avrage over max. iaux[3] values */
423
for (k=2; k<=nh; k++)
433
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
434
if (cuts[1] < rval) cuts[1] = rval;
441
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
448
rval = old; /* use last value */
451
mcc ++; /* increment null count */
458
for (n=indx+1; n<indx+nn; n++)
459
if (rval > z[n]) rval = z[n];
465
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
466
if (cuts[1] < rval) cuts[1] = rval;
474
/* -----------------------*/
475
/* here for the maximum */
476
/* -----------------------*/
479
if (iaux[3] > 0) /* take average of MAX and the next */
480
{ /* `iaux[3]' lower values */
481
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
488
rval = old; /* use last value */
491
mcc ++; /* increment null count */
497
k = 1; /* start at index 1 ... */
498
for (n=indx; n<indx+nn; n++)
500
sortit(rbuf,nn); /* sort array */
506
sum = rbuf[nn]; /* avrage over max. iaux[3] values */
507
for (k=nn-1; k>nn-nh; k--)
517
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
518
if (cuts[1] < rval) cuts[1] = rval;
525
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
532
rval = old; /* use last value */
535
mcc ++; /* increment null count */
542
for (n=indx+1; n<indx+nn; n++)
543
if (rval < z[n]) rval = z[n];
549
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
550
if (cuts[1] < rval) cuts[1] = rval;
563
............................................................
566
here comes the code for handling the median
567
if iaux[1] = 0, use index interval [ iaux[3],iaux[4] ]
568
= 1, use data interval [ faux[0],faux[1] ]
570
............................................................
575
k = iaux[3] + iaux[4];
576
if (k > 0) /* take average over index interval */
577
{ /* [ MEDIAN-iaux[3] , MEDIAN+iaux[4] ] */
578
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
585
rval = old; /* use last value */
588
mcc ++; /* increment null count */
595
rval = (z[indx] + z[indx+1]) / 2;
596
else /* here we have: nn > 2 */
598
k = 1; /* start at index 1 ... */
599
for (n=indx; n<indx+nn; n++)
601
sortit(rbuf,nn); /* sort array */
602
nh = (nn+1) / 2; /* index of median (starting at 1 ...) */
606
if (mb > nn) mb = nn;
609
for (k=ma+1; k<=mb; k++)
617
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
618
if (cuts[1] < rval) cuts[1] = rval;
625
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
632
rval = old; /* use last value */
635
mcc ++; /* increment null count */
645
if (rval > rbuf[0]) rval = rbuf[0]; /* take min of them */
647
else /* here we have: nn > 2 */
649
k = 1; /* start at index 1 ... */
650
for (n=indx; n<indx+nn; n++)
652
sortit(rbuf,nn); /* sort array */
653
nh = (nn+1) / 2; /* index of median (starting at 1 ...) */
660
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
661
if (cuts[1] < rval) cuts[1] = rval;
669
rval = faux[0] + faux[1];
670
if (rval > 0.0) /* take average over data interval */
671
{ /* [ MEDIAN-faux[0] , MEDIAN+faux[4] ] */
672
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
679
rval = old; /* use last value */
682
mcc ++; /* increment null count */
689
rval = (z[indx] + z[indx+1]) / 2;
690
else /* here we have: nn > 2 */
692
k = 1; /* start at index 1 ... */
693
for (n=indx; n<indx+nn; n++)
695
sortit(rbuf,nn); /* sort array */
696
nh = (nn+1) / 2; /* index of median (starting at 1 ...) */
697
va = rbuf[nh] - faux[0];
698
vb = rbuf[nh] + faux[1];
701
for (k=1; k<=nn; k++)
704
if (rval > vb) break; /* array is sorted ... */
717
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
718
if (cuts[1] < rval) cuts[1] = rval;
723
else /* just use median */
725
for (cntdx=0; cntdx<size; cntdx++) /* loop through count buffer */
732
rval = old; /* use last value */
735
mcc ++; /* increment null count */
742
rval = (z[indx] + z[indx+1]) / 2;
743
else /* here we have: nn > 2 */
745
k = 1; /* start at index 1 ... */
746
for (n=indx; n<indx+nn; n++)
748
sortit(rbuf,nn); /* sort array */
749
nh = (nn+1) / 2; /* index of median (starting at 1 ...) */
756
if (cuts[0] > rval) cuts[0] = rval; /* update min + max */
757
if (cuts[1] < rval) cuts[1] = rval;
774
int imnoc, imnox, imnol[MAXIMS];
775
int felm, sizec, chunk[MAXIMS], chunkc;
776
int naxis[MAXIMS], npix[MAXIMS][3], naxisc, npixc[3];
777
int iseq, iav, nulo, zpix[3];
778
int apix[3][2], cpix[3];
780
int stat, nolin, begin;
782
int frmcnt, nulcnt, kk, m, dim3, floff;
783
int debufl=0, planesize=0;
787
char *wpntr, *tmpntr;
788
char line[84], action[4];
789
char frame[84], framec[84], catfil[84];
790
char cunit[64], ident[72], output[84];
791
static char error1[] = "operands do not match in stepsize... ";
792
static char error2[] = "operands do not overlap... ";
793
static char error3[] = "stepsizes of different sign... ";
794
static char error4[] = "catalog empty... ";
795
static char mesalloc[] = "could not allocate virtual memory...";
797
double step[MAXIMS][3], start[MAXIMS][3];
798
double stepc[3], ostep[3], endc[3];
799
double startc[3], ostart[3], oldend[3];
802
float *pntr[MAXIMS], *pntrc, *zpntr;
803
float dif, eps[3], cuts[4];
808
/* set up MIDAS environment + enable automatic error abort */
812
pntrc = zpntr = (float *) 0;
813
xpntr = (short int *) 0;
815
for (nr=0; nr<3; nr++)
834
/* get result frame, input frame list, action flag */
836
(void) SCKGETC("OUT_A",1,80,&iav,framec);
837
(void) SCKGETC("P3",1,80,&iav,line);
838
(void) SCKGETC("ACTION",1,2,&iav,action);
839
CGN_UPSTR(action); /* convert to upper case */
841
(void) SCKGETC("MID$SPEC",1,5,&iav,output);
843
if ( strcmp(output,"DEBUG") == 0) debufl = 1; /* set debug flag */
846
/* ---------------------------
848
the auxiliary arrays iaux + faux contain the following:
850
iaux[0] = 1, input data is inside result space
851
0, only initialize counter pixels if necessary
852
iaux[1] = INDEX/DATA flag - not used here
853
iaux[2] = 0, take all data
854
= 1, take only data in interval [ faux[2],faux[3] ]
855
iaux[3,4] = average limits - not used here
856
iaux[5] = 0, `nomerge' option
858
iaux[6] = frame count
859
iaux[7] = no. of current frame
860
iaux[8] = 0, use NULL value
861
= 1, use last pixel as NULL value
863
faux[0,1] = average limits if iaux[1] == 1 - not used here
864
faux[2,3] = valid_data_interval if iaux[2] = 1
866
--------------------------- */
868
for (nr=0; nr<10; nr++) iaux[nr] = 0;
871
/* get Null value, average option+limits, valid_data_interval */
873
(void) SCKGETC("P5",1,40,&iav,output);
874
if ((output[0] == '+') && (output[1] == '\0'))
875
iaux[8] = 1; /* use `last' value as Null */
878
iav = CGN_CNVT(output,2,1,npixc,&usrnul,&aostep);
880
SCETER(19,"invalid `null value' ...");
884
(void) SCKGETC("P6",1,60,&iav,output);
885
CGN_UPCOPY(frame,output,2); /* upper case -> frame */
891
else if (frame[1] == 'A')
893
else if (frame[1] == 'E')
899
begin = CGN_INDEXC(output,',')+1;
900
if ((begin > 1) && (output[begin] != '\0'))
902
(void) strcpy(output,&output[begin]);
903
if (optio < 4) /* must be MIN,n or MAX,n */
905
kk = CGN_CNVT(output,1,1,npixc,&dif,&aostep);
907
SCETER(8,"invalid syntax in option_string ...");
912
iav = CGN_CNVT(output,2,2,npixc,faux,&aostep);
914
SCETER(8,"invalid syntax in `option_string' ...");
915
if (iav == 1) faux[1] = faux[0];
918
/* skip max. two numbers */
920
m = (int) strlen(output); /* or: MAX,n */
921
kk = CGN_EXTRSS(output,m,',',&begin,cunit,40);
923
kk = CGN_EXTRSS(output,m,',',&begin,cunit,40);
925
if ((output[begin] == 'D') || (output[begin] == 'd'))
926
iaux[1] = 1; /* show that it's DATA method */
928
if ((iaux[3] < 0) || (iaux[4] < 0))
929
SCETER(9,"negative `average limits' ...");
934
/* now let's look if we have a valid data interval */
936
(void) SCKGETC("P7",1,60,&iav,output);
937
if (output[0] != '+')
939
kk = CGN_CNVT(output,2,2,npixc,&faux[2],&aostep);
941
SCETER(10,"invalid syntax in `valid_data_interval' ...");
942
if (faux[2] > faux[3])
943
SCETER(11,"invalid `valid_data_interval' ...");
944
iaux[2] = 1; /* indicate we have a valid_data_interval */
948
/* test, if we have list of frames or catalog */
950
kk = CGN_INDEXS(line,".cat");
951
if (kk <= 0) kk = CGN_INDEXS(line,".CAT");
954
/* here we handle input from a catalog - get name of first image file
959
if ((int)strlen(line) > 63)
960
SCETER(3,"catalog name too long...");
962
(void) strcpy(catfil,line);
965
for (frmcnt=0; frmcnt<MAXIMS; frmcnt++) /* max. MAXIMS frames */
967
(void) SCCGET(catfil,0,frame,output,&iseq);
968
if (frame[0] == ' ') break; /* indicates the end ... */
971
(void) SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
973
sprintf(output,"%d images from catalog to be processed",frmcnt);
978
/* here we handle input from single input line - pull out operands */
983
m = (int) strlen(line);
984
for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
986
kk = CGN_EXTRSS(line,m,',',&begin,output,60);
989
CGN_FRAME(output,F_IMA_TYPE,frame,0); /* convert frames */
991
stat = SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
995
if (frmcnt <= 0) SCETER(4,error4); /* there must be at least 1 image */
998
/* --------------------------------------------*/
999
/* get initial area from 1. frame */
1000
/* --------------------------------------------*/
1002
(void) SCDRDI(imnol[0],"NAXIS",1,1,&iav,&naxis[0],&uni,&nulo);
1003
(void) SCDRDI(imnol[0],"NPIX",1,3,&iav,&npix[0][0],&uni,&nulo);
1004
(void) SCDRDD(imnol[0],"START",1,3,&iav,&start[0][0],&uni,&nulo);
1005
(void) SCDRDD(imnol[0],"STEP",1,3,&iav,&step[0][0],&uni,&nulo);
1006
(void) SCDGETC(imnol[0],"CUNIT",1,64,&iav,cunit);
1009
for (nr=0; nr<3; nr++)
1011
ostart[nr] = start[0][nr];
1012
ostep[nr] = step[0][nr];
1014
if (aostep < 0.0) aostep = - aostep;
1015
eps[nr] = 0.0001 * aostep; /* take 0.01% of stepsize */
1016
oldend[nr] = ostart[nr] + (npix[0][nr]-1)*ostep[nr];
1018
cuts[0] = 999999.0; /* cuts[0] > cuts[1] */
1019
cuts[1] = -999999.0;
1022
/* now loop through the other input frames */
1024
for (nr=1; nr< frmcnt; nr++)
1026
for (m=0; m<3; m++) /* init values... */
1033
(void) SCDRDI(imnol[nr],"NAXIS",1,1,&iav,&naxis[nr],&uni,&nulo);
1034
(void) SCDRDI(imnol[nr],"NPIX",1,3,&iav,&npix[nr][0],&uni,&nulo);
1035
(void) SCDRDD(imnol[nr],"START",1,3,&iav,&start[nr][0],&uni,&nulo);
1036
(void) SCDRDD(imnol[nr],"STEP",1,3,&iav,&step[nr][0],&uni,&nulo);
1039
/* stepsizes should have same sign and not differ too much... */
1043
if ((ostep[m]*step[nr][m]) <= 0.) SCETER(1,error3);
1044
aostep = step[nr][m] - ostep[m];
1045
if (aostep < 0.) aostep = -aostep;
1047
if (aostep > eps[m]) SCETER(5,error1);
1051
/* get intersection or union of image areas */
1053
if (action[0] != 'M')
1057
dif = start[nr][m] + (npix[nr][m]-1)*step[nr][m];
1060
if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m]; /* MAX */
1061
if (oldend[m] > dif) oldend[m] = dif; /* MIN */
1065
if (ostart[m] > start[nr][m]) ostart[m] = start[nr][m]; /* MIN */
1066
if (oldend[m] < dif) oldend[m] = dif; /* MAX */
1074
dif = start[nr][m] + (npix[nr][m]-1)*step[nr][m];
1077
if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m]; /* MAX */
1078
if (oldend[m] > dif) oldend[m] = dif; /* MIN */
1082
if (ostart[m] > start[nr][m]) ostart[m] = start[nr][m]; /* MIN */
1083
if (oldend[m] < dif) oldend[m] = dif; /* MAX */
1092
/* test, if something is left... */
1096
if (ostep[m]*(oldend[m]-ostart[m]) < 0.) SCETER(2,error2);
1100
/* create new result frame with dimension as intersection
1101
or union of input frames */
1104
if (action[0] == 'M')
1106
for (nr=1; nr<frmcnt; nr++)
1108
if (naxisc < naxis[nr]) naxisc = naxis[nr]; /* MAX */
1113
for (nr=1; nr<frmcnt; nr++)
1115
if (naxisc > naxis[nr]) naxisc = naxis[nr]; /* MIN */
1120
/* check, if input from a cube */
1131
/* set up standard stuff for result frame */
1134
for (nr=0; nr<naxisc; nr++)
1136
startc[nr] = ostart[nr];
1137
stepc[nr] = ostep[nr];
1138
npixc[nr] = CGN_NINT((oldend[nr]-ostart[nr]) / stepc[nr]) + 1;
1139
sizec = sizec * npixc[nr];
1143
(void) SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc);
1144
(void) SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni);
1145
(void) SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni);
1146
(void) SCDWRD(imnoc,"START",startc,1,naxisc,&uni);
1147
(void) SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni);
1148
(void) strcpy(ident,"average frame ");
1149
(void) SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni);
1150
mm = (naxisc+1) * 16;
1151
(void) SCDWRC(imnoc,"CUNIT",1,cunit,1,mm,&uni);
1154
/* in case of debugging save the counts in MIDAS image averdumy.dum */
1159
(void) strcpy(frame,"averdumy.dum");
1160
(void) SCFCRE(frame,D_I2_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnox);
1161
(void) SCDWRI(imnox,"NAXIS",&naxisc,1,1,&uni);
1162
(void) SCDWRI(imnox,"NPIX",npixc,1,naxisc,&uni);
1163
(void) SCDWRD(imnox,"START",startc,1,naxisc,&uni);
1164
(void) SCDWRD(imnox,"STEP",stepc,1,naxisc,&uni);
1165
(void) strcpy(ident,"valid pixel counts ");
1166
(void) SCDWRC(imnox,"IDENT",1,ident,1,72,&uni);
1167
(void) SCDWRC(imnox,"CUNIT",1,cunit,1,mm,&uni);
1171
/* see, if keyword MONITPAR(20) holds amount of virtual memory available */
1173
(void) SCKRDI("MONITPAR",20,1,&iav,&chunkc,&uni,&nulo);
1174
mm = (chunkc * chunkc) / 5; /* `chunkc' = x-dim of square image
1175
which can be fully mapped */
1177
stripe = 1; /* at least a single line */
1179
stripe = mm / npixc[0];
1182
/* get virtual memory buffers */
1184
if (npixc[1] < stripe) stripe = npixc[1];
1185
chunkc = stripe * npixc[0]; /* size of 1 result strip */
1186
kk = chunkc * sizeof(float);
1187
wpntr = malloc((size_t)kk);
1188
if (wpntr == (char *) 0)
1189
SCETER(66,mesalloc);
1191
pntrc = (float *) wpntr;
1195
if (npix[0][1] < stripe)
1196
iav = npix[0][1]; /* iav = y-dim */
1200
chunk[0] = iav * npix[0][0]; /* always stripe lines */
1201
for (nr=0; nr<dim3; nr++)
1202
{ /* all chunks like chunk[0] */
1203
chunk[nr] = chunk[0];
1204
kk = chunk[nr] * sizeof(float);
1205
wpntr = malloc((size_t)kk);
1206
if (wpntr == (char *) 0)
1207
SCETER(66,mesalloc);
1209
pntr[nr] = (float *) wpntr;
1210
imnol[nr] = imnol[0];
1211
for (m=0; m<naxisc; m++)
1213
npix[nr][m] = npix[0][m];
1214
start[nr][m] = start[0][m];
1215
step[nr][m] = step[0][m];
1221
for (nr=0; nr<frmcnt; nr++)
1223
if (npix[nr][1] < stripe)
1224
iav = npix[nr][1]; /* iav = y-dim */
1228
chunk[nr] = iav * npix[nr][0]; /* always stripe lines */
1229
kk = chunk[nr] * sizeof(float);
1230
wpntr = malloc((size_t)kk);
1231
if (wpntr == (char *) 0)
1232
SCETER(66,mesalloc);
1234
pntr[nr] = (float *) wpntr;
1239
/* now map chunk for z-direction and count buffer */
1243
planesize = npixc[0] * npixc[1];
1248
m = frmcnt * chunkc;
1249
kk = m * sizeof(float);
1250
wpntr = malloc((size_t)kk);
1251
if (wpntr == (char *) 0)
1252
SCETER(66,mesalloc);
1254
zpntr = (float *) wpntr;
1256
kk = chunkc * sizeof(short int);
1257
wpntr = malloc((size_t)kk);
1258
if (wpntr == (char *) 0)
1259
SCETER(66,mesalloc);
1261
xpntr = (short int *) wpntr;
1264
/* here the main loops over all chunks
1265
first fill the cube chunk, work on it + store result */
1267
endc[0] = startc[0] + (npixc[0]-1)*stepc[0];
1268
endc[1] = startc[1] + (stripe-1)*stepc[1];
1272
for (nolin=0; nolin<npixc[1]; nolin+=stripe)
1274
if ((nolin+stripe) > npixc[1]) /* adjust chunk size */
1276
stripe = npixc[1] - nolin;
1277
chunkc = stripe * npixc[0];
1278
for (nr=0; nr<frmcnt; nr++)
1279
chunk[nr] = stripe * npix[nr][0];
1283
floff = 0; /* for data cube input */
1284
for (nr=0; nr<frmcnt; nr++)
1287
/* convert start + end of overlap region into pixel no.'s */
1289
for (m=0; m<naxisc; m++)
1291
dif = (start[nr][m]-startc[m]) / stepc[m];
1296
cpix[m] = CGN_NINT(dif); /* offset in output */
1297
if (cpix[m] >= zpix[m])
1304
dif = (startc[m]-start[nr][m]) / step[nr][m];
1305
if (dif < 0.1) /* offset in input */
1309
apix[m][0] = CGN_NINT(dif);
1310
if (apix[m][0] >= npix[nr][m])
1317
dif = (endc[m]-start[nr][m]) / step[nr][m];
1318
apix[m][1] = CGN_NINT(dif);
1319
if (apix[m][1] >= npix[nr][m]) apix[m][1] = npix[nr][m] - 1;
1323
/* remember, SCFGET/SCFPUT begins with 1! */
1324
felm = floff + (apix[1][0] * npix[nr][0]) + apix[0][0] + 1;
1325
tmpntr = (char *) pntr[nr];
1326
stat = SCFGET(imnol[nr],felm,chunk[nr],&iav,tmpntr);
1329
sect_5360: /* fill z-buffer */
1331
fill(iaux,faux,pntr[nr],xpntr,zpntr,apix,cpix,npix[nr][0],zpix);
1336
/* now do the calculus */
1338
add(optio,iaux,faux,xpntr,zpntr,pntrc,usrnul,cuts,zpix,&kk);
1341
/* and write results to disk */
1343
felm = nolin*npixc[0] + 1;
1344
tmpntr = (char *) pntrc;
1345
(void) SCFPUT(imnoc,felm,chunkc,tmpntr);
1346
nulcnt += kk; /* update null count */
1349
tmpntr = (char *) xpntr;
1350
(void) SCFPUT(imnox,felm,chunkc,tmpntr); /* if debug, update on disk */
1354
/* follow chunks with start + end value */
1356
aostep = stripe*stepc[1];
1357
startc[1] += aostep;
1362
/* update descriptors + cuts of result frame */
1366
CGN_DSCUPD(imnol[0],imnoc,frame);
1370
(void) SCDWRR(imnoc,"LHCUTS",cuts,1,4,&uni);
1372
cuts[0] = nulcnt; /* update key NULL */
1380
(void) SCKWRR("NULL",cuts,1,mm,&uni);
1384
(void) sprintf(output,
1385
"%d undefined pixels, set to `null value' (= %12.6f)",nulcnt,usrnul);
1387
(void) sprintf(output,
1388
"%d undefined pixels, set to `previous pixel'",nulcnt);
added
removed
prim/general/src/averag.c
