← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/spec/src/splocext.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:
stdred/spec/src/splocext.c
 
1
/* @(#)splocext.c       19.1 (ES0-DMD) 02/25/03 14:29:44 */
 
2
/*===========================================================================
 
3
  Copyright (C) 1995 European Southern Observatory (ESO)
 
4
 
 
5
  This program is free software; you can redistribute it and/or 
 
6
  modify it under the terms of the GNU General Public License as 
 
7
  published by the Free Software Foundation; either version 2 of 
 
8
  the License, or (at your option) any later version.
 
9
 
 
10
  This program is distributed in the hope that it will be useful,
 
11
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
  GNU General Public License for more details.
 
14
 
 
15
  You should have received a copy of the GNU General Public 
 
16
  License along with this program; if not, write to the Free 
 
17
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
  MA 02139, USA.
 
19
 
 
20
  Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
        Internet e-mail: midas@eso.org
 
22
        Postal address: European Southern Observatory
 
23
                        Data Management Division 
 
24
                        Karl-Schwarzschild-Strasse 2
 
25
                        D 85748 Garching bei Muenchen 
 
26
                        GERMANY
 
27
===========================================================================*/
 
28
 
 
29
/* Program : splocext.c
 
30
/* Author  : P. Ballester   -     ESO Garching                    */
 
31
/* Date    : 08.09.94      Version 1.0                            */
 
32
/*                                                                */
 
33
/* Purpose :                                                      */
 
34
/*                                                                */
 
35
/* Row by row detection of local extrema (e.g. for wavelet transforms) */
 
36
/*                                                                */
 
37
/* Input:                                                         */  
 
38
/*             - name of input image            : IN_A            */
 
39
/*             - name of output peak detection  : OUT_A           */
 
40
/*             - name of output window mask     : OUT_B           */
 
41
 
 
42
#include <math.h>
 
43
#include <midas_def.h>
 
44
#define  ipos2D(col,row,npix) row*npix[0]+col
 
45
 
 
46
float    *fvector();
 
47
void     free_fvector();
 
48
 
 
49
main ()
 
50
 
 
51
{
 
52
 
 
53
    char   inframe[100], outpeak[100], outzero[100];
 
54
    char   outwind[100], outmask[100];
 
55
    char   cunit[64], ident[72];
 
56
    int    imnin, imnp, imnz, imnw, imnm, imns, imnb, unit, null;
 
57
    int    actvals, naxis, npix[2];
 
58
    float  *pntrin, *pntrp, *pntrw, *pntrz, *pntrm, *pntrs, *pntrb;
 
59
    double start[2], step[2];
 
60
 
 
61
    int    sign_ind1, sign_ind2, slope_ind, i, j, x, row;
 
62
    int    sign_previous, sign_current, rows[2];
 
63
    float  previous, current, next, zero_cross;
 
64
    float  max_value, abs_value, kappa;
 
65
 
 
66
    SCSPRO ("splocext");
 
67
 
 
68
    SCKGETC ("IN_A",   1, 60, &actvals, inframe);
 
69
    SCKGETC ("IN_B",   1, 60, &actvals, outpeak);
 
70
    SCKGETC ("OUT_A",  1, 60, &actvals, outzero);
 
71
    SCKGETC ("OUT_B",  1, 60, &actvals, outwind);
 
72
    SCKGETC ("INPUTC", 1, 20, &actvals, outmask);
 
73
    SCKRDI  ("INPUTI", 1, 2,  &actvals, rows,   &unit, &null);
 
74
    SCKRDR  ("INPUTR", 1, 1,  &actvals, &kappa, &unit, &null);
 
75
 
 
76
    strcpy(ident,""), strcpy(cunit,"");
 
77
 
 
78
    SCIGET(inframe, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 2, 
 
79
           &naxis, npix, start, step, ident, cunit, (char **)&pntrin, &imnin);
 
80
 
 
81
 
 
82
    if (naxis == 1) npix[1]=1;
 
83
 
 
84
    strcpy(ident,"Local extrema detection from: ");
 
85
    strcat(ident,inframe);
 
86
    SCIPUT(outpeak, D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
 
87
           npix, start, step, ident, cunit, (char **)&pntrp, &imnp);
 
88
 
 
89
    strcpy(ident,"Zero Crossings from: ");
 
90
    strcat(ident,inframe);
 
91
    SCIPUT(outzero, D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
 
92
           npix, start, step, ident, cunit, (char **)&pntrz, &imnz);
 
93
 
 
94
    strcpy(ident,"Windows from: ");
 
95
    strcat(ident,inframe);
 
96
    SCIPUT(outwind, D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
 
97
           npix, start, step, ident, cunit, (char **)&pntrw, &imnw);
 
98
 
 
99
    strcpy(ident,"Segmented local extrema from: ");
 
100
    strcat(ident,inframe);
 
101
    SCIPUT("segment.bdf", D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
 
102
           npix, start, step, ident, cunit, (char **)&pntrs, &imns);
 
103
 
 
104
    strcpy(ident,"Window fusion buffer from: ");
 
105
    strcat(ident,inframe);
 
106
    SCIPUT("buffer.bdf", D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, naxis, 
 
107
           npix, start, step, ident, cunit, (char **)&pntrb, &imnb);
 
108
 
 
109
    strcpy(ident,"Features mask from: ");
 
110
    strcat(ident,inframe);
 
111
    SCIPUT(outmask, D_R4_FORMAT, F_IO_MODE, F_IMA_TYPE, 1,
 
112
           npix, start, step, ident, cunit, (char **)&pntrm, &imnm);
 
113
 
 
114
    for (row=0; row<npix[1]; row++) {
 
115
         for (i=0; i<npix[0]; i++)  {
 
116
 
 
117
              if (i>0)         previous = pntrin[ipos2D(i-1,row,npix)];
 
118
                               current  = pntrin[ipos2D(i  ,row,npix)];
 
119
              if (i<npix[0]-1) next     = pntrin[ipos2D(i+1,row,npix)];
 
120
 
 
121
              if (i == 0)          previous = next;
 
122
              if (i == npix[0]-1)  next     = previous;
 
123
 
 
124
              /* Detection of local extrema */
 
125
              slope_ind = compare(current,previous) + compare(current,next);
 
126
              if (slope_ind == 2 || slope_ind == -2)
 
127
                   pntrp[ipos2D(i,row,npix)] = current;
 
128
              else
 
129
                  pntrp[ipos2D(i,row,npix)] = 0.;
 
130
 
 
131
              /* Detection of zero-crossings */
 
132
              zero_cross = 0.;
 
133
              if (i>0 && i<npix[0]-1) {
 
134
                 sign_previous = (previous > 0.) ? 1 : -1;
 
135
                 sign_current  = (current  > 0.) ? 1 : -1;
 
136
                 if ((sign_previous + sign_current) == 0) 
 
137
                     zero_cross = (float) compare(previous,current);
 
138
               }
 
139
              pntrz[ipos2D(i,row,npix)] = zero_cross;
 
140
 
 
141
            }
 
142
 
 
143
         /* Segmentation of each row to generate window mask */
 
144
         for (i=0; i<npix[0]; i++) {
 
145
                pntrm[i] = pntrp[ipos2D(i,row,npix)];
 
146
                pntrw[ipos2D(i,row,npix)] = 0.;
 
147
                pntrb[ipos2D(i,row,npix)] = 0.;
 
148
              }
 
149
 
 
150
         segment(pntrm, npix[0], kappa);
 
151
 
 
152
         for (i=0; i<npix[0]; i++) {
 
153
             pntrs[ipos2D(i,row,npix)] = pntrm[i];
 
154
             pntrw[ipos2D(i,row,npix)] = pntrm[i];
 
155
           }
 
156
 
 
157
          /* Generation of window mask */
 
158
 
 
159
           for (i=0; i<npix[0]; i++) {
 
160
            if (pntrm[i] != 0.) {
 
161
              for(j=0; 
 
162
                  (i+j)<npix[0] && pntrz[ipos2D(i+j,row,npix)]==0.
 
163
                  && (pntrp[ipos2D(i+j,row,npix)]==0. || j==0); 
 
164
                  j++)
 
165
                  pntrw[ipos2D(i+j,row,npix)] = pntrm[i];
 
166
                  pntrw[ipos2D(i+j,row,npix)] = pntrm[i];
 
167
              for(j=0; 
 
168
                  (i-j)>=0  && pntrz[ipos2D(i-j,row,npix)]==0.
 
169
                  && (pntrp[ipos2D(i+j,row,npix)]==0. || j==0); 
 
170
                  j++)
 
171
                  pntrw[ipos2D(i-j,row,npix)] = pntrm[i];
 
172
                  pntrw[ipos2D(i-j,row,npix)] = pntrm[i];
 
173
            }
 
174
         }
 
175
 
 
176
    /* Removing artifact features */
 
177
 
 
178
         for (i=1; i<npix[0]; i++) {
 
179
 
 
180
              sign_ind2 = sign(pntrw[ipos2D((i-1),row,npix)]) - 
 
181
                          sign(pntrw[ipos2D((i),  row,npix)]);
 
182
 
 
183
              sign_ind2 *= sign((float)sign_ind2);
 
184
 
 
185
              if (sign_ind2 == 2) {
 
186
 
 
187
                  printf("Artifact at %d\n",i);
 
188
                  previous = pntrw[ipos2D((i-1),row,npix)];
 
189
                  previous *= sign(previous);
 
190
                  next     = pntrw[ipos2D((i),row,npix)];
 
191
                  next     *= sign(next);
 
192
 
 
193
                  if (previous > next) {
 
194
                      for(j=0; 
 
195
                         (i+j)<npix[0] && 
 
196
                         (pntrw[ipos2D((i+j),row,npix)] != 0. || j==0);
 
197
                         j++)
 
198
                         pntrw[ipos2D((i+j),row,npix)] = 0.;
 
199
                         printf("Erased until %d\n",i+j-1);
 
200
                    }
 
201
                  else {
 
202
                      for(j=1; 
 
203
                         (i-j)>=0 && 
 
204
                         (pntrw[ipos2D((i-j),row,npix)] != 0. || j==1);
 
205
                         j++)
 
206
                         pntrw[ipos2D((i-j),row,npix)] = 0.;
 
207
                         printf("Erased until %d\n",i-j+1);
 
208
                  }}}
 
209
 
 
210
    } /* For row = .. */
 
211
 
 
212
    /* Multi-resolution scale fusion */
 
213
    /* for (row=0; row<npix[1]; row++)
 
214
        fusion_scales_down (pntrw, pntrb, pntrz, npix, row); 
 
215
     for (i=0; i<npix[0]; i++) pntrm[i] = pntrb[ipos2D(i,(npix[1]-1),npix)]; */
 
216
 
 
217
    for (row=npix[1]-1; row>=0; row--)
 
218
        fusion_scales_up (pntrw, pntrb, pntrz, npix, row);
 
219
     for (i=0; i<npix[0]; i++) pntrm[i] = pntrb[ipos2D(i,0,npix)];
 
220
 
 
221
    SCSEPI();
 
222
  }
 
223
 
 
224
/***********************************************************************/
 
225
/* Subroutine: compare(a, b) */
 
226
 
 
227
static int compare(a,b)
 
228
float a,b;
 
229
{
 
230
int ret;
 
231
ret = (a >= b) ? -1 : 1;
 
232
if (a == b)     ret = 0;
 
233
return(ret);
 
234
}
 
235
 
 
236
 
 
237
/***********************************************************************/
 
238
/* Subroutine: sign(a) */
 
239
 
 
240
static int sign(a)
 
241
float a;
 
242
{
 
243
if (a<0.) return(-1);
 
244
if (a>0.) return(1);
 
245
if (a==0.) return(0);
 
246
}
 
247
 
 
248
/*************************************************************************/
 
249
/*  Sorting algorithm:      Heapsort method                      */
 
250
/*  From: Numerical Recipes, Cambridge University Press, p. 226  */
 
251
 
 
252
static int sort(n,ra)
 
253
 
 
254
int n;
 
255
float ra[];
 
256
{
 
257
        int l,j,ir,i;
 
258
        float rra;
 
259
 
 
260
        l=(n >> 1)+1;
 
261
        ir=n;
 
262
        for (;;) {
 
263
                if (l > 1)
 
264
                        rra=ra[--l];
 
265
                else {
 
266
                        rra=ra[ir];
 
267
                        ra[ir]=ra[1];
 
268
                        if (--ir == 1) {
 
269
                                ra[1]=rra;
 
270
                                return;
 
271
                        }
 
272
                }
 
273
                i=l;
 
274
                j=l << 1;
 
275
                while (j <= ir) {
 
276
                        if (j < ir && ra[j] < ra[j+1]) ++j;
 
277
                        if (rra < ra[j]) {
 
278
                                ra[i]=ra[j];
 
279
                                j += (i=j);
 
280
                        }
 
281
                        else j=ir+1;
 
282
                }
 
283
                ra[i]=rra;
 
284
        }
 
285
}
 
286
 
 
287
/***********************************************************************/
 
288
static float median(List, N)
 
289
 
 
290
float *List;
 
291
int   N;
 
292
 
 
293
{
 
294
float *Histo, med;
 
295
int   pos, i;
 
296
 
 
297
pos = (N+1)/2;
 
298
Histo = fvector(1, N);
 
299
for (i=1; i<=N; i++) Histo[i] = List[i-1];
 
300
sort(N, Histo);
 
301
med = Histo[pos];
 
302
free_fvector(Histo, 1, N);
 
303
return(med);
 
304
 
 
305
}
 
306
 
 
307
/***********************************************************************/
 
308
static float mean(List, N)
 
309
 
 
310
float *List;
 
311
int   N;
 
312
 
 
313
{
 
314
float av;
 
315
int   pos, i;
 
316
 
 
317
av = 0.;
 
318
for (i=0; i<N; i++) av += List[i];
 
319
av /= (float)N;
 
320
return(av);
 
321
 
 
322
}
 
323
 
 
324
/***********************************************************************/
 
325
/* Subroutine update_window(pntrh, npixh, pntrwin, npix, pos) */
 
326
 
 
327
static int update_window(pntrh, npixh, pntrwin, npix, pos, row)
 
328
float *pntrh, *pntrwin;
 
329
int   npix[], npixh[], pos, row;
 
330
{
 
331
int i = pos, sign_ref;
 
332
float   peak;
 
333
 
 
334
peak = pntrh[ipos2D(pos, row, npixh)]/2.;
 
335
sign_ref = compare(peak,pntrh[ipos2D(pos, row, npix)]);
 
336
 
 
337
while(i < npix[0] && sign_ref == compare(peak,pntrh[ipos2D(i, row, npixh)]))
 
338
          pntrwin[ipos2D(i++, row, npix)] = 1.;
 
339
i = pos;
 
340
while(i>0 && sign_ref == compare(peak,pntrh[ipos2D(i, row, npixh)]))
 
341
          pntrwin[ipos2D(i--, row, npix)] = 1.;
 
342
}
 
343
 
 
344
/***********************************************************************/
 
345
/* Subroutine: Segment(List, N, kappa); */
 
346
 
 
347
static int segment(List, N, kappa)
 
348
float *List, kappa;
 
349
int   N;
 
350
{
 
351
float m, s, t[2], *dev, median();
 
352
int   i, cnt=0, cnt0;
 
353
 
 
354
dev = fvector(0, N-1);
 
355
 
 
356
/* Ignore zeroed values and copies absolute deviations */
 
357
for (i=0; i<N; i++) if (List[i] != 0.) {
 
358
   dev[cnt] = List[i]; 
 
359
   cnt++;
 
360
 }
 
361
 
 
362
 
 
363
m   = median(dev, cnt);
 
364
for (i=0; i<cnt; i++) dev[i] -= m, dev[i] *= dev[i];
 
365
s =  median(dev, cnt);
 
366
s = (float) sqrt((double)s);
 
367
 
 
368
t[0] = m - kappa*s;
 
369
t[1] = m + kappa*s;
 
370
 
 
371
cnt0 = cnt;
 
372
cnt = 0;
 
373
for (i=0; i<N; i++) {
 
374
  if (List[i] > t[0] && List[i] < t[1]) List[i] = 0.;
 
375
  else cnt++, List[i] -= m;
 
376
}
 
377
 
 
378
printf("m = %f  s = %f Range [%f, %f] Kept %d out of %d\n",
 
379
        m,s,t[0],t[1],cnt,cnt0);
 
380
 
 
381
free_fvector(dev, 0, N-1);
 
382
}
 
383
 
 
384
 
 
385
 
 
386
/***********************************************************************/
 
387
static int fusion_scales_down (w, b, z, npix, row)
 
388
 
 
389
float *w, *b, *z;
 
390
int   npix[], row;
 
391
 
 
392
{
 
393
 
 
394
    int       i,j, k, small_sign, large_sign, max_value = 0;
 
395
    float     small_value, large_value;
 
396
    float     D_Scale, Scale;
 
397
 
 
398
    if (row == 0) {
 
399
           for (i=0; i<npix[0]; i++) 
 
400
               b[ipos2D(i, row, npix)] = w[ipos2D(i,row,npix)];
 
401
         }
 
402
    else {
 
403
 
 
404
    Scale     = pow(2.,(double)(row));
 
405
    k         = 4.*Scale;
 
406
 
 
407
    for (i=0; i<k; i++) 
 
408
         b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row-1), npix)];
 
409
 
 
410
    for (i=npix[0]-k; i<npix[0]; i++) 
 
411
         b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row-1), npix)];
 
412
 
 
413
    for (i=k; i<npix[0]-k; i++) {
 
414
          if (w[ipos2D(i, row, npix)] == 0.) {
 
415
              max_value=0;
 
416
              b[ipos2D(i, row, npix)] = b[ipos2D(i, (row-1), npix)];
 
417
            }
 
418
          else {
 
419
             if (max_value == 0) {
 
420
                  max_value = 1;
 
421
                  for (j=0; 
 
422
                       (((i+j)<npix[0]) && 
 
423
                       (z[ipos2D((i+j),(row),npix)])==0.);
 
424
                       j++) {
 
425
                     small_value = b[ipos2D((i+j)  , (row-1), npix)];
 
426
                     small_value *= (float) sign(small_value);
 
427
                     large_value = w[ipos2D((i+j)  , row, npix)];
 
428
                     large_value *= (float) sign(large_value);
 
429
                     if (small_value > large_value) max_value = -1;
 
430
                   }}
 
431
 
 
432
              b[ipos2D(i, row, npix)] = 0.;
 
433
              if (max_value == -1) 
 
434
                 b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row-1), npix)];
 
435
              if (max_value == 1) 
 
436
                  b[ipos2D(i, row, npix)] =  w[ipos2D(i, row, npix)];
 
437
 
 
438
           }}}
 
439
  }
 
440
/***********************************************************************/
 
441
static int fusion_scales_up   (w, b, z, npix, row)
 
442
 
 
443
float *w, *b, *z;
 
444
int   npix[], row;
 
445
 
 
446
{
 
447
 
 
448
    int       i,j, k, small_sign, large_sign, max_value = 0;
 
449
    float     small_value, large_value;
 
450
    float     D_Scale, Scale;
 
451
 
 
452
    if (row == npix[1]-1) {
 
453
           for (i=0; i<npix[0]; i++) 
 
454
               b[ipos2D(i, row, npix)] = w[ipos2D(i,row,npix)];
 
455
         }
 
456
    else {
 
457
 
 
458
    Scale     = pow(2.,(double)(row));
 
459
    k         = 4.*Scale;
 
460
 
 
461
    for (i=0; i<k; i++) 
 
462
         b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row+1), npix)];
 
463
 
 
464
    for (i=npix[0]-k; i<npix[0]; i++) 
 
465
         b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row+1), npix)];
 
466
 
 
467
    for (i=k; i<npix[0]-k; i++) {
 
468
          if (w[ipos2D(i, (row  ), npix)] == 0.  && 
 
469
              b[ipos2D(i, (row+1), npix)] == 0.) {
 
470
              max_value=0;
 
471
              b[ipos2D(i, row, npix)] = b[ipos2D(i, (row+1), npix)]; 
 
472
            }
 
473
          else {
 
474
             if (max_value == 0) {
 
475
                  max_value = 1;
 
476
                  for (j=0;
 
477
                       (((i+j)<npix[0]) && 
 
478
                       (b[ipos2D((i+j),(row+1),npix)] !=0. || 
 
479
                        w[ipos2D((i+j),(row  ),npix)] !=0.));
 
480
                       j++) {
 
481
                     small_value = w[ipos2D((i+j)  , (row), npix)];
 
482
                     small_sign  = sign(small_value);
 
483
                     large_value = b[ipos2D((i+j)  , (row+1), npix)];
 
484
                     large_sign  = sign(large_value);
 
485
                     if (small_value*small_sign > large_value*large_sign &&
 
486
                         (small_sign + large_sign) != 0) 
 
487
                            max_value = -1;
 
488
                   }}
 
489
 
 
490
              b[ipos2D(i, row, npix)] = 0.;
 
491
              if (max_value == -1) 
 
492
                 b[ipos2D(i, row, npix)] =  w[ipos2D(i, (row), npix)];
 
493
              if (max_value == 1) 
 
494
                  b[ipos2D(i, row, npix)] =  b[ipos2D(i, (row+1), npix)];
 
495
 
 
496
           }}}
 
497
  }