← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/general/src/averagw.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/general/src/averagw.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1989-2009 European Southern Observatory (ESO)
 
3
 
 
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.
 
8
 
 
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.
 
13
 
 
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, 
 
17
  MA 02139, USA.
 
18
 
 
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 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
 
 
29
/* ++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
 
 
31
.IDENTIFICATION
 
32
  program AVERAGW                       version 1.00    090729
 
33
  K. Banse                              ESO - Garching
 
34
 
 
35
.KEYWORDS
 
36
  bulk data frames, average
 
37
 
 
38
.PURPOSE
 
39
  average a no. of weighted image frames, the result will either be a frame
 
40
  with size equal to common area of all frames involved or 
 
41
  size equal to union of all frame areas
 
42
 
 
43
.ALGORITHM
 
44
  extract input frames either from given parameter string or from catalogue,
 
45
  add up all frames + divide in the end
 
46
 
 
47
.INPUT/OUTPUT
 
48
  the following keys are used:
 
49
 
 
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
 
54
 
 
55
.VERSIONS
 
56
 
 
57
 090729         last modif
 
58
-------------------------------------------------- */
 
59
 
 
60
#include <midas_def.h>
 
61
#include <stdio.h>
 
62
#include <stdlib.h>
 
63
 
 
64
#define  MAXIMS  300
 
65
#define  MAXIMSONE  MAXIMS+1
 
66
 
 
67
static double wsum;
 
68
static float usrnul;
 
69
 
 
70
 
 
71
/*
 
72
 
 
73
*/
 
74
 
 
75
void wfill(iaux,ww,a,x,z,apix,cpix,npixa,npixc)
 
76
int    *iaux, apix[3][2], *cpix;
 
77
int    npixa, *npixc;
 
78
float   *x;
 
79
float   ww, *a, *z;
 
80
 
 
81
{
 
82
int    frmcnt, count;
 
83
int    nn, nin, nini, jump, indx, cntdx;
 
84
int    nx, ny, xfirst, yfirst, xend, yend;
 
85
register int nr;
 
86
 
 
87
register float  rr;
 
88
 
 
89
 
 
90
 
 
91
/* ---------------------------
 
92
 
 
93
the intermediate z-space is structured as follows:
 
94
   
 
95
pix(1) pix(2) up to  pix(FRMCNT)      1              )
 
96
pix(1) pix(2) up to  pix(FRMCNT)      2              )
 
97
. . .                                                ) = 1 line
 
98
pix(1) pix(2) up to  pix(FRMCNT)      NPIXC(1)       )
 
99
  
 
100
STRIP lines of above
 
101
--------------------------- */
 
102
 
 
103
 
 
104
/*   if here for the first time, init the count `pixels'  */
 
105
 
 
106
frmcnt = iaux[6];
 
107
count = iaux[7];
 
108
 
 
109
if (count == 0)
 
110
   {
 
111
   nn = npixc[0] * npixc[1];
 
112
   if ((iaux[5] == 0) && (iaux[2] == 0))
 
113
      rr = (float) wsum;                /* for `nomerge' and `all data' */
 
114
   else
 
115
      rr = 0.0;
 
116
   for (nr=0; nr<nn; nr++) x[nr] = rr;
 
117
   }
 
118
 
 
119
if (iaux[0] == 0) return;               /* if iaux[0] = 0, not much to do */
 
120
 
 
121
nini = 0;                          /* input frame begins with 1. valid pixel */
 
122
if (iaux[5] == 0)
 
123
   {
 
124
   /* --------------------------------- */
 
125
   /*  here for the non-merging option  */
 
126
   /* --------------------------------- */
 
127
 
 
128
   indx = count;                /* take all pixels */
 
129
   for (nn=0; nn<npixc[1]; nn++)
 
130
      {
 
131
      nin = nini;
 
132
      for (nr=0; nr<npixc[0]; nr++)
 
133
         {
 
134
         z[indx] = ww * a[nin++];
 
135
         indx += frmcnt;
 
136
         }
 
137
      nini += npixa;                     /* follow with input index */
 
138
      }
 
139
   }
 
140
else
 
141
   {
 
142
   /* ----------------------------- */
 
143
   /*  here for the merging option  */
 
144
   /* ----------------------------- */
 
145
 
 
146
   nx = apix[0][1] - apix[0][0];                /*   get overlapping part  */
 
147
   ny = apix[1][1] - apix[1][0];
 
148
   xfirst = cpix[0];
 
149
   yfirst = cpix[1];
 
150
   xend = xfirst + nx;
 
151
   yend = yfirst + ny;
 
152
   jump = frmcnt * npixc[0];
 
153
 
 
154
   cntdx = 0;
 
155
   indx = count;
 
156
 
 
157
   for (nn=0; nn<npixc[1]; nn++)
 
158
      {
 
159
      z[indx] = usrnul;
 
160
      if ((nn >= yfirst) && (nn <= yend)) 
 
161
         {
 
162
         nin = nini;
 
163
         for (nr=0; nr<npixc[0]; nr++)
 
164
            {
 
165
            if ((nr >= xfirst) && (nr <= xend)) 
 
166
               {
 
167
               z[indx] = ww * a[nin++]; 
 
168
               x[cntdx] += ww;          /* increment weight sum */
 
169
               }
 
170
            indx += frmcnt;
 
171
            cntdx ++;
 
172
            }
 
173
         nini += npixa;                 /* follow with input index */
 
174
         }
 
175
      else
 
176
         {
 
177
         indx += jump;                  /* jump a whole line; */
 
178
         cntdx += npixc[0];
 
179
         }
 
180
      }
 
181
   }
 
182
/*
 
183
for (cntdx=0; cntdx<5; cntdx++)
 
184
   {
 
185
   printf("x[%d] = %f\n",cntdx,x[cntdx]);
 
186
   }
 
187
*/
 
188
 
 
189
}
 
190
 
 
191
/*
 
192
 
 
193
*/
 
194
 
 
195
void wadd(flag,iaux,x,z,c,cuts,npixc,nc)
 
196
int   *npixc;
 
197
int   flag, *iaux, *nc;
 
198
float   *x, *z, *c, *cuts;
 
199
 
 
200
{
 
201
int    frmcnt, n, mcc, indx, size;
 
202
register int cntdx;
 
203
 
 
204
float  ws, rval;
 
205
  
 
206
double      sum;
 
207
 
 
208
 
 
209
 
 
210
mcc = 0;                                /*  init  */
 
211
indx = 0;
 
212
frmcnt = iaux[6];
 
213
size = npixc[0] * npixc[1];
 
214
  
 
215
 
 
216
   /* -------------------------------*/
 
217
   /*   here for the normal average  */
 
218
   /* -------------------------------*/
 
219
   
 
220
   for (cntdx=0; cntdx<size; cntdx++)           /* loop through count buffer */
 
221
      {
 
222
      ws = x[cntdx];
 
223
  
 
224
      if (ws < 0.000001) 
 
225
         {
 
226
         rval = usrnul;
 
227
         mcc ++;                                /* increment null count  */
 
228
         } 
 
229
      else
 
230
         {
 
231
         sum = z[indx];
 
232
         for (n=indx+1; n<indx+frmcnt; n++) sum += z[n];
 
233
         rval = sum / ws;
 
234
         }
 
235
 
 
236
      c[cntdx] = rval;
 
237
      if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
 
238
      if (cuts[1] < rval) cuts[1] = rval;
 
239
      indx += frmcnt;
 
240
      }
 
241
 
 
242
*nc = mcc;
 
243
}
 
244
 
 
245
/*
 
246
 
 
247
*/
 
248
 
 
249
int main()
 
250
 
 
251
{
 
252
int    uni, ec, ed, el, ln0, ln1;
 
253
int    imnoc, imnox, imnol[MAXIMS];
 
254
int    felm, sizec, chunk[MAXIMS], chunkc;
 
255
int    naxis[MAXIMS], npix[MAXIMS][3], naxisc, npixc[3];
 
256
int    iseq, iav, nulo, zpix[3];
 
257
int    apix[3][2], cpix[3];
 
258
int    iaux[10];
 
259
int    stat, nolin, begin;
 
260
int    stripe, optio;
 
261
int    frmcnt, nulcnt, kk, m, dim3, floff;
 
262
int    debufl=0, planesize=0;
 
263
register int nr;
 
264
 
 
265
char     *wpntr, *tmpntr;
 
266
char     line[84], action[4];
 
267
char     frame[84], framec[84], catfil[84];
 
268
char     cunit[64], ident[72], output[84];
 
269
static char  error1[] = "operands do not match in stepsize... ";
 
270
static char  error2[] = "operands do not overlap... ";
 
271
static char  error3[] = "stepsizes of different sign... ";
 
272
static char  error4[] = "catalog empty... ";
 
273
static char  mesalloc[] = "could not allocate virtual memory...";
 
274
 
 
275
double   step[MAXIMS][3], start[MAXIMS][3];
 
276
double   stepc[3], ostep[3], endc[3];
 
277
double   startc[3], ostart[3], oldend[3];
 
278
double   aostep;
 
279
 
 
280
float    *pntr[MAXIMS], *pntrc, *xpntr, *zpntr;
 
281
float    dif, eps[3], cuts[4], w[MAXIMS];
 
282
 
 
283
 
 
284
/* set up MIDAS environment + enable automatic error abort */
 
285
 
 
286
SCSPRO("averag");
 
287
 
 
288
pntrc = zpntr = (float *) 0;
 
289
xpntr = (float *) 0;
 
290
 
 
291
for (nr=0; nr<3; nr++)
 
292
   {
 
293
   apix[nr][0] = 0;
 
294
   apix[nr][1] = 0;
 
295
   cpix[nr] = 0;
 
296
   startc[nr] = 0.0;
 
297
   stepc[nr] = 1.0;
 
298
   endc[nr] = 0.0;
 
299
   npixc[nr] = 1;
 
300
   zpix[nr] = 1;
 
301
   ostart[nr] = 0.0;
 
302
   oldend[nr] = 0.0;
 
303
   ostep[nr] = 1.0;
 
304
   npix[0][nr] = 1;
 
305
   start[0][nr] = 0.0;
 
306
   step[0][nr] = 1.0;
 
307
   }
 
308
 
 
309
 
 
310
/* get result frame, input frame list, action flag */
 
311
 
 
312
(void) SCKGETC("OUT_A",1,80,&iav,framec);
 
313
(void) SCKGETC("P3",1,80,&iav,line);
 
314
(void) SCKGETC("ACTION",1,2,&iav,action);
 
315
CGN_UPSTR(action);                      /* convert to upper case */
 
316
nulcnt = 0;
 
317
(void) SCKGETC("MID$SPEC",1,5,&iav,output);
 
318
output[5] = '\0';
 
319
if ( strcmp(output,"DEBUG") == 0) debufl = 1;           /* set debug flag */
 
320
 
 
321
 
 
322
/* ---------------------------
 
323
 
 
324
the auxiliary arrays iaux contains the following:
 
325
 
 
326
iaux[0] = 1, input data is inside result space
 
327
          0, only initialize counter pixels if necessary
 
328
iaux[1,2,3,4] = not used 
 
329
iaux[5] = 0, `nomerge' option
 
330
        = 1, `merge' option
 
331
iaux[6] = frame count
 
332
iaux[7-9] = not used
 
333
 
 
334
--------------------------- */
 
335
 
 
336
for (nr=0; nr<10; nr++) iaux[nr] = 0;
 
337
 
 
338
 
 
339
/* get Null value */
 
340
 
 
341
(void) SCKGETC("P5",1,40,&iav,output);
 
342
if ((output[0] == '+') && (output[1] == '\0'))          /* use key NULL(2) */
 
343
   (void) SCKRDR("NULL",2,1,&iav,&usrnul,&uni,&nulo);
 
344
else
 
345
   {
 
346
   iav = CGN_CNVT(output,2,1,npixc,&usrnul,&aostep);
 
347
   if (iav < 1)
 
348
      SCETER(19,"invalid `null value' ...");
 
349
   }
 
350
 
 
351
optio = 1;
 
352
 
 
353
/*  test, if we have list of frames or catalog */
 
354
 
 
355
kk = CGN_INDEXS(line,".cat");
 
356
if (kk <= 0) kk = CGN_INDEXS(line,".CAT");
 
357
 
 
358
 
 
359
/*   here we handle input from a catalog - get name of first image file 
 
360
     in catalog                                                         */
 
361
 
 
362
if (kk > 0) 
 
363
   {
 
364
   if ((int)strlen(line) > 63)
 
365
      SCETER(3,"catalog name too long...");
 
366
   else
 
367
      (void) strcpy(catfil,line);
 
368
 
 
369
   iseq = 0;
 
370
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)          /* max. MAXIMS frames */
 
371
      {
 
372
      (void) SCCGET(catfil,0,frame,output,&iseq);
 
373
      if (frame[0] == ' ') break;               /*  indicates the end ... */
 
374
    
 
375
      imnol[frmcnt] = 0;
 
376
      (void) SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
 
377
      }
 
378
   sprintf(output,"%d images from catalog to be processed",frmcnt);
 
379
   SCTPUT(output);
 
380
   }
 
381
 
 
382
 
 
383
/*  here we handle input from single input line - pull out operands */
 
384
 
 
385
else
 
386
   {
 
387
   begin = 0;
 
388
   m = (int) strlen(line);
 
389
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
 
390
      {
 
391
      kk = CGN_EXTRSS(line,m,',',&begin,output,60);
 
392
      if (kk <= 0) break;
 
393
       
 
394
      CGN_FRAME(output,F_IMA_TYPE,frame,0);     /* convert frames  */
 
395
      imnol[frmcnt] = 0;
 
396
      stat = SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
 
397
      }
 
398
   }
 
399
 
 
400
if (frmcnt <= 0) SCETER(4,error4);       /* there must be at least 1 image  */
 
401
 
 
402
 
 
403
/*  --------------------------------------------*/
 
404
/*   get initial area from 1. frame             */
 
405
/*  --------------------------------------------*/
 
406
 
 
407
(void) SCDRDI(imnol[0],"NAXIS",1,1,&iav,&naxis[0],&uni,&nulo);
 
408
(void) SCDRDI(imnol[0],"NPIX",1,3,&iav,&npix[0][0],&uni,&nulo);
 
409
(void) SCDRDD(imnol[0],"START",1,3,&iav,&start[0][0],&uni,&nulo);
 
410
(void) SCDRDD(imnol[0],"STEP",1,3,&iav,&step[0][0],&uni,&nulo);
 
411
(void) SCDGETC(imnol[0],"CUNIT",1,64,&iav,cunit);
 
412
 
 
413
 
 
414
for (nr=0; nr<3; nr++)
 
415
   {
 
416
   ostart[nr] = start[0][nr];
 
417
   ostep[nr] = step[0][nr];
 
418
   aostep = ostep[nr];
 
419
   if (aostep < 0.0) aostep = - aostep;
 
420
   eps[nr] = 0.0001 * aostep;                   /* take 0.01% of stepsize */
 
421
   oldend[nr] = ostart[nr] + (npix[0][nr]-1)*ostep[nr];
 
422
   }
 
423
cuts[0] = 999999.0;                                     /* cuts[0] > cuts[1] */
 
424
cuts[1] = -999999.0;
 
425
 
 
426
 
 
427
/*   loop + look for descriptor WEIGHT  */
 
428
 
 
429
   stat = SCECNT("GET",&ec,&el,&ed);
 
430
   ln1 = 1; ln0 = 0;
 
431
   stat = SCECNT("PUT",&ln1,&ln0,&ln0);
 
432
 
 
433
   wsum = 0.0;
 
434
   for (nr=0; nr< frmcnt; nr++)
 
435
      {
 
436
      stat = SCDRDR(imnol[nr],"WEIGHT",1,1,&iav,&w[nr],&uni,&nulo);
 
437
      if (stat != ERR_NORMAL) 
 
438
         {
 
439
         w[nr] = 1.0;
 
440
         (void) sprintf(line,"frame no. %d: descr. WEIGHT missing - set to 1.0 ",imnol[nr]);
 
441
         SCTPUT(line);
 
442
         }
 
443
 
 
444
      wsum += w[nr];
 
445
      }
 
446
 
 
447
 
 
448
   stat = SCECNT("PUT",&ec,&el,&ed);
 
449
 
 
450
/*   now loop through the other input frames  */
 
451
  
 
452
for (nr=1; nr< frmcnt; nr++)
 
453
   {
 
454
   for (m=0; m<3; m++)                          /* init values... */
 
455
      {
 
456
      start[nr][m] = 0.0;
 
457
      step[nr][m] = 1.0;
 
458
      npix[nr][m] = 1;
 
459
      }
 
460
  
 
461
   (void) SCDRDI(imnol[nr],"NAXIS",1,1,&iav,&naxis[nr],&uni,&nulo);
 
462
   (void) SCDRDI(imnol[nr],"NPIX",1,3,&iav,&npix[nr][0],&uni,&nulo);
 
463
   (void) SCDRDD(imnol[nr],"START",1,3,&iav,&start[nr][0],&uni,&nulo);
 
464
   (void) SCDRDD(imnol[nr],"STEP",1,3,&iav,&step[nr][0],&uni,&nulo);
 
465
 
 
466
 
 
467
/*  stepsizes should have same sign and not differ too much...  */
 
468
 
 
469
   for (m=0; m<3; m++)
 
470
      {
 
471
      if ((ostep[m]*step[nr][m]) <= 0.) SCETER(1,error3);
 
472
      aostep = step[nr][m] - ostep[m];
 
473
      if (aostep < 0.) aostep = -aostep;
 
474
      
 
475
      if (aostep > eps[m]) SCETER(5,error1);
 
476
      }
 
477
 
 
478
 
 
479
/*   get intersection or union of image areas  */
 
480
 
 
481
   if (action[0] != 'M')        /* intersection */ 
 
482
      {
 
483
      for (m=0; m<3; m++)
 
484
         {
 
485
         dif = start[nr][m] + (npix[nr][m]-1)*step[nr][m];
 
486
         if (ostep[m] > 0.) 
 
487
            {
 
488
            if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m];   /* MAX */
 
489
            if (oldend[m] > dif) oldend[m] = dif;                     /* MIN */
 
490
            }
 
491
         else
 
492
            {
 
493
            if (ostart[m] > start[nr][m]) ostart[m] = start[nr][m];   /* MIN */
 
494
            if (oldend[m] < dif) oldend[m] = dif;                     /* MAX */
 
495
            }
 
496
         }
 
497
      }
 
498
   else                         /* union */
 
499
      {
 
500
      for (m=0; m<3; m++)
 
501
         {
 
502
         dif = start[nr][m] + (npix[nr][m]-1)*step[nr][m];
 
503
         if (ostep[m] < 0.) 
 
504
            {
 
505
            if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m];   /* MAX */
 
506
            if (oldend[m] > dif) oldend[m] = dif;                     /* MIN */
 
507
            }
 
508
         else
 
509
            {
 
510
            if (ostart[m] > start[nr][m]) ostart[m] = start[nr][m];   /* MIN */
 
511
            if (oldend[m] < dif) oldend[m] = dif;                     /* MAX */
 
512
            }
 
513
         }
 
514
      iaux[5] = 1;
 
515
      }
 
516
 
 
517
   }
 
518
 
 
519
 
 
520
/*  test, if something is left...  */
 
521
 
 
522
for (m=0; m<3; m++)
 
523
   {
 
524
   if (ostep[m]*(oldend[m]-ostart[m]) < 0.) SCETER(2,error2);
 
525
   }
 
526
 
 
527
 
 
528
/*   create new result frame with dimension as intersection 
 
529
     or union of input frames                                   */
 
530
 
 
531
naxisc = naxis[0];
 
532
if (action[0] == 'M') 
 
533
   {
 
534
   for (nr=1; nr<frmcnt; nr++)
 
535
      {
 
536
      if (naxisc < naxis[nr]) naxisc = naxis[nr];               /* MAX */
 
537
      }
 
538
   }
 
539
else
 
540
   {
 
541
   for (nr=1; nr<frmcnt; nr++)
 
542
      {
 
543
      if (naxisc > naxis[nr]) naxisc = naxis[nr];               /* MIN */
 
544
      }
 
545
   }
 
546
 
 
547
 
 
548
/*  check, if input from a cube */
 
549
 
 
550
if (naxisc > 2) 
 
551
   {
 
552
   dim3 = npix[0][2];
 
553
   naxisc = 2;
 
554
   }
 
555
else
 
556
   dim3 = 0;
 
557
 
 
558
 
 
559
/*  set up standard stuff for result frame        */
 
560
 
 
561
sizec = 1;
 
562
for (nr=0; nr<naxisc; nr++)
 
563
   {
 
564
   startc[nr] = ostart[nr];
 
565
   stepc[nr] = ostep[nr];
 
566
   npixc[nr] = CGN_NINT((oldend[nr]-ostart[nr]) / stepc[nr]) + 1;
 
567
   sizec = sizec * npixc[nr];
 
568
   }
 
569
 
 
570
imnoc = 0;
 
571
(void) SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc);
 
572
(void) SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni);
 
573
(void) SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni);
 
574
(void) SCDWRD(imnoc,"START",startc,1,naxisc,&uni);
 
575
(void) SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni);
 
576
(void) strcpy(ident,"average frame ");
 
577
(void) SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni);
 
578
ln1 = (naxisc+1) * 16;
 
579
(void) SCDWRC(imnoc,"CUNIT",1,cunit,1,ln1,&uni);
 
580
 
 
581
  
 
582
/* in case of debugging save the counts in MIDAS image averdumy.dum  */
 
583
 
 
584
if (debufl == 1)
 
585
   {
 
586
   imnox = 0;
 
587
   (void) strcpy(frame,"averdumy.dum");
 
588
   (void) SCFCRE(frame,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnox);
 
589
   (void) SCDWRI(imnox,"NAXIS",&naxisc,1,1,&uni);
 
590
   (void) SCDWRI(imnox,"NPIX",npixc,1,naxisc,&uni);
 
591
   (void) SCDWRD(imnox,"START",startc,1,naxisc,&uni);
 
592
   (void) SCDWRD(imnox,"STEP",stepc,1,naxisc,&uni);
 
593
   (void) strcpy(ident,"valid pixel counts ");
 
594
   (void) SCDWRC(imnox,"IDENT",1,ident,1,72,&uni);
 
595
   (void) SCDWRC(imnox,"CUNIT",1,cunit,1,ln1,&uni);
 
596
   }
 
597
 
 
598
 
 
599
/* see, if keyword MONITPAR(20) holds amount of virtual memory available */
 
600
 
 
601
(void) SCKRDI("MONITPAR",20,1,&iav,&chunkc,&uni,&nulo);
 
602
ln1 = (chunkc * chunkc) / 5;             /* `chunkc' = x-dim of square image
 
603
                                            which can be fully mapped */
 
604
if (ln1 < npixc[0])
 
605
   stripe = 1;                           /* at least a single line */
 
606
else
 
607
   stripe = ln1 / npixc[0];
 
608
 
 
609
 
 
610
/*  get virtual memory buffers */
 
611
 
 
612
if (npixc[1] < stripe) stripe = npixc[1];
 
613
chunkc = stripe * npixc[0];                 /* size of 1 result strip */
 
614
kk = chunkc * sizeof(float);
 
615
wpntr = malloc((size_t)kk);
 
616
if (wpntr == (char *) 0)
 
617
   SCETER(66,mesalloc);
 
618
else
 
619
   pntrc = (float *) wpntr;
 
620
 
 
621
if (dim3 > 1)
 
622
   {
 
623
   if (npix[0][1] < stripe) 
 
624
      iav = npix[0][1];                         /* iav = y-dim  */
 
625
   else
 
626
      iav = stripe;
 
627
 
 
628
   chunk[0] = iav * npix[0][0];         /* always stripe lines  */
 
629
   for (nr=0; nr<dim3; nr++)
 
630
      {                                 /* all chunks like chunk[0] */
 
631
      chunk[nr] = chunk[0];
 
632
      kk = chunk[nr] * sizeof(float);
 
633
      wpntr = malloc((size_t)kk);
 
634
      if (wpntr == (char *) 0)
 
635
         SCETER(66,mesalloc);
 
636
      else
 
637
         pntr[nr] = (float *) wpntr;
 
638
      w[nr] = w[0];
 
639
      imnol[nr] = imnol[0];
 
640
      for (m=0; m<naxisc; m++)
 
641
         {
 
642
         npix[nr][m] = npix[0][m];
 
643
         start[nr][m] = start[0][m];
 
644
         step[nr][m] = step[0][m];
 
645
         }
 
646
      }
 
647
   }
 
648
else
 
649
   {
 
650
   for (nr=0; nr<frmcnt; nr++)
 
651
      {
 
652
      if (npix[nr][1] < stripe) 
 
653
         iav = npix[nr][1];                             /* iav = y-dim  */
 
654
      else
 
655
         iav = stripe;
 
656
 
 
657
      chunk[nr] = iav * npix[nr][0];            /* always stripe lines */
 
658
      kk = chunk[nr] * sizeof(float);
 
659
      wpntr = malloc((size_t)kk);
 
660
      if (wpntr == (char *) 0)
 
661
         SCETER(66,mesalloc);
 
662
      else
 
663
         pntr[nr] = (float *) wpntr;
 
664
      }
 
665
   }
 
666
 
 
667
 
 
668
/*  now map chunk for z-direction and count buffer  */
 
669
 
 
670
if (dim3 > 1)
 
671
   {
 
672
   planesize = npixc[0] * npixc[1];
 
673
   frmcnt = dim3;
 
674
   }
 
675
 
 
676
iaux[6] = frmcnt;
 
677
m = frmcnt * chunkc;
 
678
kk = m * sizeof(float);
 
679
wpntr = malloc((size_t)kk);
 
680
if (wpntr == (char *) 0)
 
681
   SCETER(66,mesalloc);
 
682
else
 
683
   zpntr = (float *) wpntr;
 
684
 
 
685
kk = chunkc * sizeof(float);
 
686
wpntr = malloc((size_t)kk);
 
687
if (wpntr == (char *) 0)
 
688
   SCETER(66,mesalloc);
 
689
else
 
690
   xpntr = (float *) wpntr;
 
691
  
 
692
 
 
693
/*   here the main loops over all chunks
 
694
     first fill the cube chunk,  work on it + store result           */
 
695
 
 
696
endc[0] = startc[0] + (npixc[0]-1)*stepc[0];
 
697
endc[1] = startc[1] + (stripe-1)*stepc[1];
 
698
zpix[0] = npixc[0];
 
699
  
 
700
 
 
701
for (nolin=0; nolin<npixc[1]; nolin+=stripe)
 
702
   {
 
703
   if ((nolin+stripe) > npixc[1])               /* adjust chunk size */
 
704
      {
 
705
      stripe = npixc[1] - nolin;
 
706
      chunkc = stripe * npixc[0];
 
707
      for (nr=0; nr<frmcnt; nr++)
 
708
         chunk[nr] = stripe * npix[nr][0];
 
709
      }
 
710
   zpix[1] = stripe;
 
711
 
 
712
   floff = 0;                           /* for data cube input */
 
713
   for (nr=0; nr<frmcnt; nr++)
 
714
      {
 
715
 
 
716
      /*  convert start + end of overlap region into pixel no.'s   */
 
717
 
 
718
      for (m=0; m<naxisc; m++)
 
719
         {
 
720
         dif = (start[nr][m]-startc[m]) / stepc[m];
 
721
         if (dif < 0.1) 
 
722
            cpix[m] = 0;
 
723
         else
 
724
            {
 
725
            cpix[m] = CGN_NINT(dif);                    /* offset in output */
 
726
            if (cpix[m] >= zpix[m]) 
 
727
               {
 
728
               iaux[0] = 0;
 
729
               goto sect_5360;
 
730
               }
 
731
            }
 
732
 
 
733
         dif = (startc[m]-start[nr][m]) / step[nr][m];
 
734
         if (dif < 0.1)                                 /* offset in input */
 
735
            apix[m][0] = 0;
 
736
         else
 
737
            {
 
738
            apix[m][0] = CGN_NINT(dif);
 
739
            if (apix[m][0] >= npix[nr][m]) 
 
740
               {
 
741
               iaux[0] = 0;
 
742
               goto sect_5360;
 
743
               }
 
744
            }
 
745
         
 
746
         dif = (endc[m]-start[nr][m]) / step[nr][m];
 
747
         apix[m][1] = CGN_NINT(dif);
 
748
         if (apix[m][1] >= npix[nr][m]) apix[m][1] = npix[nr][m] - 1;
 
749
         }
 
750
  
 
751
      iaux[0] = 1;
 
752
                                  /* remember, SCFGET/SCFPUT begins with 1!  */
 
753
      felm = floff + (apix[1][0] * npix[nr][0]) + apix[0][0] + 1;
 
754
      tmpntr = (char *) pntr[nr];
 
755
      stat = SCFGET(imnol[nr],felm,chunk[nr],&iav,tmpntr);
 
756
 
 
757
 
 
758
     sect_5360:                                         /* fill z-buffer */
 
759
      iaux[7] = nr;
 
760
      wfill(iaux,w[nr],pntr[nr],xpntr,zpntr,apix,cpix,npix[nr][0],zpix);
 
761
      floff += planesize;
 
762
      }
 
763
  
 
764
 
 
765
   /*   now do the calculus  */
 
766
 
 
767
   wadd(optio,iaux,xpntr,zpntr,pntrc,cuts,zpix,&kk);
 
768
 
 
769
  
 
770
   /*  and write results to disk */
 
771
 
 
772
   felm = nolin*npixc[0] + 1;
 
773
   tmpntr = (char *) pntrc;
 
774
   (void) SCFPUT(imnoc,felm,chunkc,tmpntr);
 
775
   nulcnt += kk;                                /* update null count */
 
776
   if (debufl == 1)
 
777
      {
 
778
      tmpntr = (char *) xpntr;
 
779
      (void) SCFPUT(imnox,felm,chunkc,tmpntr);  /* if debug, update on disk */
 
780
      }
 
781
  
 
782
 
 
783
   /*  follow chunks with start + end value  */
 
784
 
 
785
   aostep = stripe*stepc[1];
 
786
   startc[1] += aostep;
 
787
   endc[1] += aostep;
 
788
   }
 
789
 
 
790
  
 
791
/*  update descriptors + cuts of result frame  */
 
792
 
 
793
frame[0] = ' ';
 
794
frame[1] = '\0';
 
795
CGN_DSCUPD(imnol[0],imnoc,frame);
 
796
 
 
797
cuts[2] = cuts[0];
 
798
cuts[3] = cuts[1];
 
799
(void) SCDWRR(imnoc,"LHCUTS",cuts,1,4,&uni);
 
800
 
 
801
cuts[0] = nulcnt;                       /* update key NULL */
 
802
if (iaux[8] == 0)
 
803
   {
 
804
   cuts[1] = usrnul;
 
805
   ln1 = 2;
 
806
   }
 
807
else
 
808
   ln1 = 1;
 
809
(void) SCKWRR("NULL",cuts,1,ln1,&uni);
 
810
if (nulcnt > 0) 
 
811
   {
 
812
   if (iaux[8] == 0)
 
813
      (void) sprintf(output,
 
814
      "%d undefined pixels, set to `null value' (= %12.6f)",nulcnt,usrnul);
 
815
   else
 
816
      (void) sprintf(output,
 
817
      "%d undefined pixels, set to `previous pixel'",nulcnt);
 
818
   SCTPUT(output);
 
819
   }
 
820
  
 
821
SCSEPI();
 
822
return(0);
 
823
}