« back to all changes in this revision
Viewing changes to src/back.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Package Import Robot
- Author(s): Gürkan Sengün
- Date: 2013-05-04 09:05:03 UTC
- Revision ID: package-import@ubuntu.com-20130504090503-aemv054m5m20oj9d
Tags: upstream-2.19.1
Import upstream version 2.19.1
- files added:
- autoconf
- doc
- man
- src
- src/fits
- src/wcs
- test
- xsl
added
removed
src/back.c
1
/*
2
back.c
3
4
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5
*
6
* Part of: SWarp
7
*
8
* Author: E.BERTIN (IAP)
9
*
10
* Contents: functions dealing with background computation.
11
*
12
* Last modify: 25/08/2010
13
*
14
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
*/
16
17
#ifdef HAVE_CONFIG_H
18
#include "config.h"
19
#endif
20
21
#ifdef HAVE_MATHIMF_H
22
#include <mathimf.h>
23
#else
24
#include <math.h>
25
#endif
26
#include <stdio.h>
27
#include <stdlib.h>
28
#include <string.h>
29
30
#include "define.h"
31
#include "globals.h"
32
#include "fits/fitscat.h"
33
#include "back.h"
34
#include "coadd.h"
35
#include "field.h"
36
#include "prefs.h"
37
#include "weight.h"
38
39
/******************************** make_back **********************************/
40
/*
41
Background maps are established from the images themselves; thus we need to
42
make at least one first pass through the data.
43
*/
44
void make_back(fieldstruct *field, fieldstruct *wfield, int wscale_flag)
45
46
{
47
backstruct *backmesh,*wbackmesh, *bm,*wbm;
48
tabstruct *tab, *wtab;
49
PIXTYPE *buf,*wbuf, *buft,*wbuft;
50
size_t bufsize, bufsize2,
51
bufshift, size,meshsize,jumpsize;
52
off_t fcurpos,wfcurpos, wfcurpos2,fcurpos2;
53
int i,j,k,m,n, step, nlines,
54
w,bw, bh, nx,ny,nb,
55
lflag, nr;
56
float *ratio,*ratiop, *weight, *sigma,
57
sratio;
58
59
/* If the weight-map is not an external one, no stats are needed for it */
60
if (wfield && (wfield->flags&BACKRMS_FIELD))
61
wfield= NULL;
62
tab = field->tab;
63
if (wfield)
64
wtab = wfield->tab;
65
else
66
wtab = NULL; /* to avoid gcc -Wall warnings */
67
w = field->width;
68
bw = field->backw;
69
bh = field->backh;
70
nx = field->nbackx;
71
ny = field->nbacky;
72
nb = field->nback;
73
74
NFPRINTF(OUTPUT, "Setting up background maps");
75
76
/* Decide if it is worth displaying progress each 16 lines */
77
78
lflag = (field->width*field->backh >= (size_t)65536);
79
80
/* Save current positions in files */
81
82
wfcurpos = wfcurpos2 = 0; /* to avoid gcc -Wall warnings */
83
QFSEEK(tab->cat->file, tab->bodypos, SEEK_SET, field->filename);
84
QFTELL(fcurpos, tab->cat->file, field->filename);
85
if (wfield)
86
{
87
QFSEEK(wtab->cat->file, wtab->bodypos, SEEK_SET, wfield->filename);
88
QFTELL(wfcurpos, wtab->cat->file, wfield->filename);
89
}
90
91
/* Allocate a correct amount of memory to store pixels */
92
93
bufsize = (size_t)w*bh;
94
meshsize = bufsize;
95
nlines = 0;
96
if (bufsize > (size_t)BACK_BUFSIZE)
97
{
98
nlines = BACK_BUFSIZE/w;
99
step = (field->backh-1)/nlines+1;
100
bufsize = (size_t)(nlines = field->backh/step)*w;
101
bufshift = (step/2)*(size_t)w;
102
jumpsize = (step-1)*(size_t)w;
103
}
104
else
105
bufshift = jumpsize = 0; /* to avoid gcc -Wall warnings */
106
107
/* Allocate some memory */
108
QMALLOC(backmesh, backstruct, nx); /* background information */
109
QMALLOC(buf, PIXTYPE, bufsize); /* pixel buffer */
110
free(field->back);
111
QMALLOC(field->back, float, nb); /* background map */
112
free(field->backline);
113
QMALLOC(field->backline, PIXTYPE, w); /* current background line */
114
free(field->sigma);
115
QMALLOC(field->sigma, float, nb); /* sigma map */
116
if (wfield)
117
{
118
QMALLOC(wbackmesh, backstruct, nx); /* background information */
119
QMALLOC(wbuf, PIXTYPE, bufsize); /* pixel buffer */
120
free(wfield->back);
121
QMALLOC(wfield->back, float, nb); /* background map */
122
free(wfield->sigma);
123
QMALLOC(wfield->sigma, float, nb); /* sigma map */
124
wfield->sigfac = 1.0;
125
set_weightconv(wfield); /* Prepare conversion */
126
}
127
else
128
{
129
wbackmesh = NULL;
130
wbuf = NULL;
131
}
132
133
/* Loop over the data packets */
134
for (j=0; j<ny; j++)
135
{
136
if (lflag && j)
137
NPRINTF(OUTPUT,
138
"\33[1M> Setting up background map at line:%7d / %-7d\n\33[1A",
139
j*bh, field->height);
140
if (!nlines)
141
{
142
/*---- The image is small enough so that we can make exhaustive stats */
143
if (j == ny-1 && field->npix%bufsize)
144
bufsize = field->npix%bufsize;
145
read_body(tab, buf, bufsize);
146
if (wfield)
147
{
148
read_body(wtab, wbuf, bufsize);
149
weight_to_var(wbuf, bufsize);
150
}
151
/*---- Build the histograms */
152
backstat(backmesh, wbackmesh, buf, wbuf, bufsize,nx, w, bw,
153
wfield?wfield->weight_thresh:0.0);
154
bm = backmesh;
155
for (m=nx; m--; bm++)
156
if (bm->mean <= -BIG)
157
bm->histo=NULL;
158
else
159
QCALLOC(bm->histo, int, bm->nlevels);
160
if (wfield)
161
{
162
wbm = wbackmesh;
163
for (m=nx; m--; wbm++)
164
if (wbm->mean <= -BIG)
165
wbm->histo=NULL;
166
else
167
QCALLOC(wbm->histo, int, wbm->nlevels);
168
}
169
backhisto(backmesh, wbackmesh, buf, wbuf, bufsize,nx, w, bw,
170
wfield?wfield->weight_thresh:0.0);
171
}
172
else
173
{
174
/*---- Image size too big, we have to skip a few data !*/
175
QFTELL(fcurpos2, tab->cat->file, field->filename);
176
if (wfield)
177
QFTELL(wfcurpos2, wtab->cat->file, wfield->filename);
178
if (j == ny-1 && (n=field->height%field->backh))
179
{
180
meshsize = n*(size_t)w;
181
nlines = BACK_BUFSIZE/w;
182
step = (n-1)/nlines+1;
183
bufsize = (nlines = n/step)*(size_t)w;
184
bufshift = (step/2)*(size_t)w;
185
jumpsize = (step-1)*(size_t)w;
186
free(buf);
187
QMALLOC(buf, PIXTYPE, bufsize); /* pixel buffer */
188
if (wfield)
189
{
190
free(wbuf);
191
QMALLOC(wbuf, PIXTYPE, bufsize); /* pixel buffer */
192
}
193
}
194
195
/*---- Read and skip, read and skip, etc... */
196
QFSEEK(tab->cat->file, bufshift*tab->bytepix, SEEK_CUR, field->filename);
197
buft = buf;
198
for (i=nlines; i--; buft += w)
199
{
200
read_body(tab, buft, w);
201
if (i)
202
QFSEEK(tab->cat->file, jumpsize*tab->bytepix, SEEK_CUR,
203
field->filename);
204
}
205
206
if (wfield)
207
{
208
/*------ Read and skip, read and skip, etc... now on the weight-map */
209
QFSEEK(wtab->cat->file,bufshift*wtab->bytepix, SEEK_CUR,
210
wfield->filename);
211
wbuft = wbuf;
212
for (i=nlines; i--; wbuft += w)
213
{
214
read_body(wtab, wbuft, w);
215
weight_to_var(wbuft, w);
216
if (i)
217
QFSEEK(wtab->cat->file, jumpsize*wtab->bytepix, SEEK_CUR,
218
wfield->filename);
219
}
220
}
221
backstat(backmesh, wbackmesh, buf, wbuf, bufsize, nx, w, bw,
222
wfield?wfield->weight_thresh:0.0);
223
QFSEEK(tab->cat->file, fcurpos2, SEEK_SET, field->filename);
224
bm = backmesh;
225
for (m=nx; m--; bm++)
226
if (bm->mean <= -BIG)
227
bm->histo=NULL;
228
else
229
QCALLOC(bm->histo, int, bm->nlevels);
230
if (wfield)
231
{
232
QFSEEK(wtab->cat->file, wfcurpos2, SEEK_SET, wfield->filename);
233
wbm = wbackmesh;
234
for (m=nx; m--; wbm++)
235
if (wbm->mean <= -BIG)
236
wbm->histo=NULL;
237
else
238
QCALLOC(wbm->histo, int, wbm->nlevels);
239
}
240
/*---- Build (progressively this time) the histograms */
241
for(size=meshsize, bufsize2=bufsize; size>0; size -= bufsize2)
242
{
243
if (bufsize2>size)
244
bufsize2 = size;
245
read_body(tab, buf, bufsize2);
246
if (wfield)
247
{
248
read_body(wtab, wbuf, bufsize2);
249
weight_to_var(wbuf, bufsize2);
250
}
251
backhisto(backmesh, wbackmesh, buf, wbuf, bufsize2, nx, w, bw,
252
wfield?wfield->weight_thresh:0.0);
253
}
254
}
255
256
/*-- Compute background statistics from the histograms */
257
bm = backmesh;
258
for (m=0; m<nx; m++, bm++)
259
{
260
k = m+nx*j;
261
backguess(bm, field->back+k, field->sigma+k);
262
free(bm->histo);
263
}
264
if (wfield)
265
{
266
wbm = wbackmesh;
267
for (m=0; m<nx; m++, wbm++)
268
{
269
k = m+nx*j;
270
backguess(wbm, wfield->back+k, wfield->sigma+k);
271
free(wbm->histo);
272
}
273
}
274
}
275
276
/* Free memory */
277
free(buf);
278
free(backmesh);
279
if (wfield)
280
{
281
free(wbackmesh);
282
free(wbuf);
283
}
284
285
/* Go back to the original position */
286
QFSEEK(field->tab->cat->file, fcurpos, SEEK_SET, field->filename);
287
if (wfield)
288
QFSEEK(wfield->tab->cat->file, wfcurpos, SEEK_SET, wfield->filename);
289
290
/* Median-filter and check suitability of the background map */
291
NFPRINTF(OUTPUT, "Filtering background map(s)");
292
filter_back(field);
293
if (wfield)
294
filter_back(wfield);
295
296
/* Compute normalization for variance- or weight-maps*/
297
if (wfield && wscale_flag && wfield->flags&(VAR_FIELD|WEIGHT_FIELD))
298
{
299
nr = 0;
300
QMALLOC(ratio, float, wfield->nback);
301
ratiop = ratio;
302
weight = wfield->back;
303
sigma = field->sigma;
304
for (i=wfield->nback; i--; sigma++)
305
if ((sratio=*(weight++)) > 0.0
306
&& (sratio = *sigma/sqrt(sratio)) > 0.0)
307
{
308
*(ratiop++) = sratio;
309
nr++;
310
}
311
wfield->sigfac = (double)hmedian(ratio, nr);
312
for (i=0; i<nr && ratio[i]<=0.0; i++);
313
if (i<nr)
314
wfield->sigfac = (double)hmedian(ratio+i, nr-i);
315
else
316
{
317
warning("Null or negative global weighting factor:","defaulted to 1");
318
wfield->sigfac = 1.0;
319
}
320
321
/*-- Update weight threshold */
322
wfield->weight_thresh *= (PIXTYPE)(wfield->sigfac*wfield->sigfac);
323
324
free(ratio);
325
}
326
327
/* Compute 2nd derivatives along the y-direction */
328
NFPRINTF(OUTPUT, "Computing backgound d-map");
329
free(field->dback);
330
field->dback = make_backspline(field, field->back);
331
NFPRINTF(OUTPUT, "Computing backgound-noise d-map");
332
free(field->dsigma);
333
field->dsigma = make_backspline(field, field->sigma);
334
/* If asked for, force the backmean parameter to the supplied value */
335
if (field->back_type == BACK_ABSOLUTE)
336
field->backmean = field->backdefault;
337
338
field->fbackmean = field->backmean*field->fascale;
339
field->fbacksig = field->backsig*field->fascale;
340
341
return;
342
}
343
344
345
/******************************** backstat **********************************/
346
/*
347
Compute robust statistical estimators in a row of meshes.
348
*/
349
void backstat(backstruct *backmesh, backstruct *wbackmesh,
350
PIXTYPE *buf, PIXTYPE *wbuf, size_t bufsize,
351
int n, int w, int bw, PIXTYPE wthresh)
352
353
{
354
backstruct *bm, *wbm;
355
double pix,wpix, sig, mean,wmean, sigma,wsigma, step;
356
PIXTYPE *buft,*wbuft, lcut,wlcut, hcut,whcut;
357
int m,h,x,y, npix,wnpix, offset, lastbite;
358
359
h = bufsize/w;
360
bm = backmesh;
361
wbm = wbackmesh;
362
offset = w - bw;
363
step = sqrt(2/PI)*QUANTIF_NSIGMA/QUANTIF_AMIN;
364
wmean = wsigma = wlcut = whcut = 0.0; /* to avoid gcc -Wall warnings */
365
for (m = n; m--; bm++,buf+=bw)
366
{
367
if (!m && (lastbite=w%bw))
368
{
369
bw = lastbite;
370
offset = w-bw;
371
}
372
mean = sigma = 0.0;
373
buft=buf;
374
/*-- We separate the weighted case at this level to avoid penalty in CPU */
375
npix = 0;
376
if (wbackmesh)
377
{
378
wmean = wsigma = 0.0;
379
wbuft = wbuf;
380
for (y=h; y--; buft+=offset,wbuft+=offset)
381
for (x=bw; x--;)
382
{
383
pix = *(buft++);
384
if ((wpix = *(wbuft++)) < wthresh && pix > -BIG)
385
{
386
wmean += wpix;
387
wsigma += wpix*wpix;
388
mean += pix;
389
sigma += pix*pix;
390
npix++;
391
}
392
}
393
}
394
else
395
{
396
for (y=h; y--; buft+=offset)
397
for (x=bw; x--;)
398
if ((pix = *(buft++)) > -BIG)
399
{
400
mean += pix;
401
sigma += pix*pix;
402
npix++;
403
}
404
}
405
406
/*-- If not enough valid pixels, discard this mesh */
407
if ((float)npix < (float)(bw*h*BACK_MINGOODFRAC))
408
{
409
bm->mean = bm->sigma = -BIG;
410
if (wbackmesh)
411
{
412
wbm->mean = wbm->sigma = -BIG;
413
wbm++;
414
wbuf += bw;
415
}
416
continue;
417
}
418
if (wbackmesh)
419
{
420
wmean /= (double)npix;
421
wsigma = (sig = wsigma/npix - wmean*wmean)>0.0? sqrt(sig):0.0;
422
wlcut = wbm->lcut = (PIXTYPE)(wmean - 2.0*wsigma);
423
whcut = wbm->hcut = (PIXTYPE)(wmean + 2.0*wsigma);
424
}
425
mean /= (double)npix;
426
sigma = (sig = sigma/npix - mean*mean)>0.0? sqrt(sig):0.0;
427
lcut = bm->lcut = (PIXTYPE)(mean - 2.0*sigma);
428
hcut = bm->hcut = (PIXTYPE)(mean + 2.0*sigma);
429
mean = sigma = 0.0;
430
npix = wnpix = 0;
431
buft = buf;
432
if (wbackmesh)
433
{
434
wmean = wsigma = 0.0;
435
wbuft=wbuf;
436
for (y=h; y--; buft+=offset, wbuft+=offset)
437
for (x=bw; x--;)
438
{
439
pix = *(buft++);
440
if ((wpix = *(wbuft++))<wthresh && pix<=hcut && pix>=lcut)
441
{
442
mean += pix;
443
sigma += pix*pix;
444
npix++;
445
if (wpix<=whcut && wpix>=wlcut)
446
{
447
wmean += wpix;
448
wsigma += wpix*wpix;
449
wnpix++;
450
}
451
}
452
}
453
}
454
else
455
for (y=h; y--; buft+=offset)
456
for (x=bw; x--;)
457
{
458
pix = *(buft++);
459
if (pix<=hcut && pix>=lcut)
460
{
461
mean += pix;
462
sigma += pix*pix;
463
npix++;
464
}
465
}
466
467
bm->npix = npix;
468
mean /= (double)npix;
469
sig = sigma/npix - mean*mean;
470
sigma = sig>0.0 ? sqrt(sig):0.0;
471
bm->mean = mean;
472
bm->sigma = sigma;
473
if ((bm->nlevels = (int)(step*npix+1)) > QUANTIF_NMAXLEVELS)
474
bm->nlevels = QUANTIF_NMAXLEVELS;
475
bm->qscale = sigma>0.0? 2*QUANTIF_NSIGMA*sigma/bm->nlevels : 1.0;
476
bm->qzero = mean - QUANTIF_NSIGMA*sigma;
477
if (wbackmesh)
478
{
479
wbm->npix = wnpix;
480
wmean /= (double)wnpix;
481
sig = wsigma/wnpix - wmean*wmean;
482
wsigma = sig>0.0 ? sqrt(sig):0.0;
483
wbm->mean = wmean;
484
wbm->sigma = wsigma;
485
if ((wbm->nlevels = (int)(step*wnpix+1)) > QUANTIF_NMAXLEVELS)
486
wbm->nlevels = QUANTIF_NMAXLEVELS;
487
wbm->qscale = wsigma>0.0? 2*QUANTIF_NSIGMA*wsigma/wbm->nlevels : 1.0;
488
wbm->qzero = wmean - QUANTIF_NSIGMA*wsigma;
489
wbm++;
490
wbuf += bw;
491
}
492
}
493
494
return;
495
}
496
497
498
/******************************** backhisto *********************************/
499
/*
500
Compute robust statistical estimators in a row of meshes.
501
*/
502
void backhisto(backstruct *backmesh, backstruct *wbackmesh,
503
PIXTYPE *buf, PIXTYPE *wbuf, size_t bufsize,
504
int n, int w, int bw, PIXTYPE wthresh)
505
{
506
backstruct *bm,*wbm;
507
PIXTYPE *buft,*wbuft,
508
pix;
509
float qscale,wqscale, cste,wcste, wpix;
510
int *histo,*whisto;
511
int h,m,x,y, nlevels,wnlevels, lastbite, offset, bin;
512
513
h = bufsize/w;
514
bm = backmesh;
515
wbm = wbackmesh;
516
offset = w - bw;
517
for (m=0; m++<n; bm++ , buf+=bw)
518
{
519
if (m==n && (lastbite=w%bw))
520
{
521
bw = lastbite;
522
offset = w-bw;
523
}
524
/*-- Skip bad meshes */
525
if (bm->mean <= -BIG)
526
{
527
if (wbackmesh)
528
{
529
wbm++;
530
wbuf += bw;
531
}
532
continue;
533
}
534
nlevels = bm->nlevels;
535
histo = bm->histo;
536
qscale = bm->qscale;
537
cste = 0.499999 - bm->qzero/qscale;
538
buft = buf;
539
if (wbackmesh)
540
{
541
wnlevels = wbm->nlevels;
542
whisto = wbm->histo;
543
wqscale = wbm->qscale;
544
wcste = 0.499999 - wbm->qzero/wqscale;
545
wbuft = wbuf;
546
for (y=h; y--; buft+=offset, wbuft+=offset)
547
for (x=bw; x--;)
548
{
549
bin = (int)((pix=*(buft++))/qscale + cste);
550
if ((wpix = *(wbuft++))<wthresh && pix>-BIG && bin<nlevels && bin>=0)
551
{
552
(*(histo+bin))++;
553
bin = (int)(wpix/wqscale + wcste);
554
if (bin>=0 && bin<wnlevels)
555
(*(whisto+bin))++;
556
}
557
}
558
wbm++;
559
wbuf += bw;
560
}
561
else
562
for (y=h; y--; buft += offset)
563
for (x=bw; x--;)
564
{
565
bin = (int)((pix=*(buft++))/qscale + cste);
566
if (bin>=0 && bin<nlevels && pix>-BIG)
567
(*(histo+bin))++;
568
}
569
}
570
571
return;
572
}
573
574
/******************************* backguess **********************************/
575
/*
576
Estimate the background from a histogram;
577
*/
578
float backguess(backstruct *bkg, float *mean, float *sigma)
579
580
#define EPS (1e-4) /* a small number */
581
582
{
583
int *histo, *hilow, *hihigh, *histot;
584
unsigned int lowsum, highsum, sum;
585
double ftemp, mea, sig, sig1, med, dpix;
586
int i, n, lcut,hcut, nlevelsm1, pix;
587
588
/* Leave here if the mesh is already classified as `bad' */
589
if (bkg->mean<=-BIG)
590
{
591
*mean = *sigma = -BIG;
592
return -BIG;
593
}
594
595
histo = bkg->histo;
596
hcut = nlevelsm1 = bkg->nlevels-1;
597
lcut = 0;
598
599
sig = 10.0*nlevelsm1;
600
sig1 = 1.0;
601
mea = med = bkg->mean;
602
for (n=100; n-- && (sig>=0.1) && (fabs(sig/sig1-1.0)>EPS);)
603
{
604
sig1 = sig;
605
sum = mea = sig = 0.0;
606
lowsum = highsum = 0;
607
histot = hilow = histo+lcut;
608
hihigh = histo+hcut;
609
for (i=lcut; i<=hcut; i++)
610
{
611
if (lowsum<highsum)
612
lowsum += *(hilow++);
613
else
614
highsum += *(hihigh--);
615
sum += (pix = *(histot++));
616
mea += (dpix = (double)pix*i);
617
sig += dpix*i;
618
}
619
620
med = hihigh>=histo?
621
((hihigh-histo)+0.5+((double)highsum-lowsum)/(2.0*(*hilow>*hihigh?
622
*hilow:*hihigh)))
623
: 0.0;
624
625
if (sum)
626
{
627
mea /= (double)sum;
628
sig = sig/sum - mea*mea;
629
}
630
sig = sig>0.0?sqrt(sig):0.0;
631
lcut = (ftemp=med-3.0*sig)>0.0 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5):0;
632
hcut = (ftemp=med+3.0*sig)<nlevelsm1 ?(int)(ftemp>0.0?ftemp+0.5:ftemp-0.5)
633
: nlevelsm1;
634
}
635
636
*mean = fabs(sig)>0.0? (fabs(bkg->sigma/(sig*bkg->qscale)-1) < 0.0 ?
637
bkg->qzero+mea*bkg->qscale
638
:(fabs((mea-med)/sig)< 0.3 ?
639
bkg->qzero+(2.5*med-1.5*mea)*bkg->qscale
640
:bkg->qzero+med*bkg->qscale))
641
:bkg->qzero+mea*bkg->qscale;
642
643
*sigma = sig*bkg->qscale;
644
645
646
return *mean;
647
}
648
649
650
/******************************* filter_back *********************************/
651
/*
652
Median filtering of the background map to remove the contribution from bright
653
sources.
654
*/
655
void filter_back(fieldstruct *field)
656
657
{
658
float *back,*sigma, *back2,*sigma2, *bmask,*smask, *sigmat,
659
d2,d2min, fthresh, med, val,sval;
660
int i,j,px,py, np, nx,ny, npx,npx2, npy,npy2, dpx,dpy, x,y, nmin;
661
662
fthresh = prefs.back_fthresh;
663
nx = field->nbackx;
664
ny = field->nbacky;
665
np = field->nback;
666
npx = field->nbackfx/2;
667
npy = field->nbackfy/2;
668
npy *= nx;
669
670
QMALLOC(bmask, float, (2*npx+1)*(2*npy+1));
671
QMALLOC(smask, float, (2*npx+1)*(2*npy+1));
672
QMALLOC(back2, float, np);
673
QMALLOC(sigma2, float, np);
674
675
back = field->back;
676
sigma = field->sigma;
677
val = sval = 0.0; /* to avoid gcc -Wall warnings */
678
679
/* Look for `bad' meshes and interpolate them if necessary */
680
for (i=0,py=0; py<ny; py++)
681
for (px=0; px<nx; px++,i++)
682
if ((back2[i]=back[i])<=-BIG)
683
{
684
/*------ Seek the closest valid mesh */
685
d2min = BIG;
686
nmin = 0;
687
for (j=0,y=0; y<ny; y++)
688
for (x=0; x<nx; x++,j++)
689
if (back[j]>-BIG)
690
{
691
d2 = (float)(x-px)*(x-px)+(y-py)*(y-py);
692
if (d2<d2min)
693
{
694
val = back[j];
695
sval = sigma[j];
696
nmin = 1;
697
d2min = d2;
698
}
699
else if (d2==d2min)
700
{
701
val += back[j];
702
sval += sigma[j];
703
nmin++;
704
}
705
}
706
back2[i] = nmin? val/nmin: 0.0;
707
sigma[i] = nmin? sval/nmin: 1.0;
708
}
709
memcpy(back, back2, (size_t)np*sizeof(float));
710
711
/* Do the actual filtering */
712
for (py=0; py<np; py+=nx)
713
{
714
npy2 = np - py - nx;
715
if (npy2>npy)
716
npy2 = npy;
717
if (npy2>py)
718
npy2 = py;
719
for (px=0; px<nx; px++)
720
{
721
npx2 = nx - px - 1;
722
if (npx2>npx)
723
npx2 = npx;
724
if (npx2>px)
725
npx2 = px;
726
i=0;
727
for (dpy = -npy2; dpy<=npy2; dpy+=nx)
728
{
729
y = py+dpy;
730
for (dpx = -npx2; dpx <= npx2; dpx++)
731
{
732
x = px+dpx;
733
bmask[i] = back[x+y];
734
smask[i++] = sigma[x+y];
735
}
736
}
737
if (fabs((med=hmedian(bmask, i))-back[px+py])>=fthresh)
738
{
739
back2[px+py] = med;
740
sigma2[px+py] = hmedian(smask, i);
741
}
742
else
743
{
744
back2[px+py] = back[px+py];
745
sigma2[px+py] = sigma[px+py];
746
}
747
}
748
}
749
750
free(bmask);
751
free(smask);
752
memcpy(back, back2, np*sizeof(float));
753
field->backmean = (double)hmedian(back2, np);
754
free(back2);
755
memcpy(sigma, sigma2, np*sizeof(float));
756
field->backsig = (double)hmedian(sigma2, np);
757
758
if (field->backsig<=0.0)
759
{
760
sigmat = sigma2+np;
761
for (i=np; i-- && *(--sigmat)>0.0;);
762
if (i>=0 && i<(np-1))
763
field->backsig = hmedian(sigmat+1, np-1-i);
764
else
765
{
766
if (field->flags&(DETECT_FIELD|MEASURE_FIELD))
767
warning("Image contains mainly constant data; ",
768
"I'll try to cope with that...");
769
field->backsig = 1.0;
770
}
771
}
772
773
free(sigma2);
774
775
776
return;
777
}
778
779
780
/******************************* make_backspline *****************************/
781
/*
782
Pre-compute 2nd derivatives along the y direction at background nodes.
783
*/
784
float *make_backspline(fieldstruct *field, float *map)
785
786
{
787
int x,y, nbx,nby,nbym1;
788
float *dmap,*dmapt,*mapt, *u, temp;
789
790
nbx = field->nbackx;
791
nby = field->nbacky;
792
nbym1 = nby - 1;
793
QMALLOC(dmap, float, field->nback);
794
for (x=0; x<nbx; x++)
795
{
796
mapt = map+x;
797
dmapt = dmap+x;
798
if (nby>1)
799
{
800
QMALLOC(u, float, nbym1); /* temporary array */
801
*dmapt = *u = 0.0; /* "natural" lower boundary condition */
802
mapt += nbx;
803
for (y=1; y<nbym1; y++, mapt+=nbx)
804
{
805
temp = -1/(*dmapt+4);
806
*(dmapt += nbx) = temp;
807
temp *= *(u++) - 6*(*(mapt+nbx)+*(mapt-nbx)-2**mapt);
808
*u = temp;
809
}
810
*(dmapt+=nbx) = 0.0; /* "natural" upper boundary condition */
811
for (y=nby-2; y--;)
812
{
813
temp = *dmapt;
814
dmapt -= nbx;
815
*dmapt = (*dmapt*temp+*(u--))/6.0;
816
}
817
free(u);
818
}
819
else
820
*dmapt = 0.0;
821
}
822
823
return dmap;
824
}
825
826
827
/******************************** backline **********************************/
828
/*
829
Interpolate background at line y (bicubic spline interpolation between
830
background map vertices).
831
*/
832
void backline(fieldstruct *field, int y, PIXTYPE *line)
833
834
{
835
int i,j,x,yl, nbx,nbxm1,nby, nx,width, ystep, changepoint;
836
float dx,dx0,dy,dy3, cdx,cdy,cdy3, temp, xstep,
837
*node,*nodep,*dnode, *blo,*bhi,*dblo,*dbhi, *u;
838
PIXTYPE bval;
839
840
width = field->width;
841
842
if (field->back_type==BACK_ABSOLUTE)
843
{
844
/*-- In absolute background mode, just subtract a cste */
845
bval = (PIXTYPE)field->backmean;
846
for (i=width; i--;)
847
*(line++) = bval;
848
return;
849
}
850
851
nbx = field->nbackx;
852
nbxm1 = nbx - 1;
853
nby = field->nbacky;
854
if (nby > 1)
855
{
856
dy = (float)y/field->backh - 0.5;
857
dy -= (yl = (int)dy);
858
if (yl<0)
859
{
860
yl = 0;
861
dy -= 1.0;
862
}
863
else if (yl>=nby-1)
864
{
865
yl = nby<2 ? 0 : nby-2;
866
dy += 1.0;
867
}
868
/*-- Interpolation along y for each node */
869
cdy = 1 - dy;
870
dy3 = (dy*dy*dy-dy);
871
cdy3 = (cdy*cdy*cdy-cdy);
872
ystep = nbx*yl;
873
blo = field->back + ystep;
874
bhi = blo + nbx;
875
dblo = field->dback + ystep;
876
dbhi = dblo + nbx;
877
QMALLOC(node, float, nbx); /* Interpolated background */
878
nodep = node;
879
for (x=nbx; x--;)
880
*(nodep++) = cdy**(blo++) + dy**(bhi++) + cdy3**(dblo++) + dy3**(dbhi++);
881
882
/*-- Computation of 2nd derivatives along x */
883
QMALLOC(dnode, float, nbx); /* 2nd derivative along x */
884
if (nbx>1)
885
{
886
QMALLOC(u, float, nbxm1); /* temporary array */
887
*dnode = *u = 0.0; /* "natural" lower boundary condition */
888
nodep = node+1;
889
for (x=nbxm1; --x; nodep++)
890
{
891
temp = -1/(*(dnode++)+4);
892
*dnode = temp;
893
temp *= *(u++) - 6*(*(nodep+1)+*(nodep-1)-2**nodep);
894
*u = temp;
895
}
896
*(++dnode) = 0.0; /* "natural" upper boundary condition */
897
for (x=nbx-2; x--;)
898
{
899
temp = *(dnode--);
900
*dnode = (*dnode*temp+*(u--))/6.0;
901
}
902
free(u);
903
dnode--;
904
}
905
}
906
else
907
{
908
/*-- No interpolation and no new 2nd derivatives needed along y */
909
node = field->back;
910
dnode = field->dback;
911
}
912
913
/*-- Interpolation along x */
914
if (nbx>1)
915
{
916
nx = field->backw;
917
xstep = 1.0/nx;
918
changepoint = nx/2;
919
dx = (xstep - 1)/2; /* dx of the first pixel in the row */
920
dx0 = ((nx+1)%2)*xstep/2; /* dx of the 1st pixel right to a bkgnd node */
921
blo = node;
922
bhi = node + 1;
923
dblo = dnode;
924
dbhi = dnode + 1;
925
for (x=i=0,j=width; j--; i++, dx += xstep)
926
{
927
if (i==changepoint && x>0 && x<nbxm1)
928
{
929
blo++;
930
bhi++;
931
dblo++;
932
dbhi++;
933
dx = dx0;
934
}
935
cdx = 1 - dx;
936
*(line++) = (PIXTYPE)(cdx*(*blo+(cdx*cdx-1)**dblo)
937
+ dx*(*bhi+(dx*dx-1)**dbhi));
938
if (i==nx)
939
{
940
x++;
941
i = 0;
942
}
943
}
944
}
945
else
946
for (j=width; j--;)
947
*(line++) = (PIXTYPE)*node;
948
949
if (nby>1)
950
{
951
free(node);
952
free(dnode);
953
}
954
955
return;
956
}
957
958
959
/******************************* backrmsline ********************************
960
PROTO void backrmsline(fieldstruct *field, int y, PIXTYPE *line)
961
PURPOSE Bicubic-spline interpolation of the background noise along the current
962
scanline (y).
963
INPUT Measurement or detection field pointer,
964
Current line position.
965
Where to put the data.
966
OUTPUT -.
967
NOTES Most of the code is a copy of subbackline(), for optimization reasons.
968
AUTHOR E. Bertin (IAP & Leiden & ESO)
969
VERSION 02/02/98
970
***/
971
void backrmsline(fieldstruct *field, int y, PIXTYPE *line)
972
973
{
974
int i,j,x,yl, nbx,nbxm1,nby, nx,width, ystep, changepoint;
975
float dx,dx0,dy,dy3, cdx,cdy,cdy3, temp, xstep,
976
*node,*nodep,*dnode, *blo,*bhi,*dblo,*dbhi, *u;
977
978
nbx = field->nbackx;
979
nbxm1 = nbx - 1;
980
nby = field->nbacky;
981
if (nby > 1)
982
{
983
dy = (float)y/field->backh - 0.5;
984
dy -= (yl = (int)dy);
985
if (yl<0)
986
{
987
yl = 0;
988
dy -= 1.0;
989
}
990
else if (yl>=nby-1)
991
{
992
yl = nby<2 ? 0 : nby-2;
993
dy += 1.0;
994
}
995
/*-- Interpolation along y for each node */
996
cdy = 1 - dy;
997
dy3 = (dy*dy*dy-dy);
998
cdy3 = (cdy*cdy*cdy-cdy);
999
ystep = nbx*yl;
1000
blo = field->sigma + ystep;
1001
bhi = blo + nbx;
1002
dblo = field->dsigma + ystep;
1003
dbhi = dblo + nbx;
1004
QMALLOC(node, float, nbx); /* Interpolated background */
1005
nodep = node;
1006
for (x=nbx; x--;)
1007
*(nodep++) = cdy**(blo++) + dy**(bhi++) + cdy3**(dblo++) + dy3**(dbhi++);
1008
1009
/*-- Computation of 2nd derivatives along x */
1010
QMALLOC(dnode, float, nbx); /* 2nd derivative along x */
1011
if (nbx>1)
1012
{
1013
QMALLOC(u, float, nbxm1); /* temporary array */
1014
*dnode = *u = 0.0; /* "natural" lower boundary condition */
1015
nodep = node+1;
1016
for (x=nbxm1; --x; nodep++)
1017
{
1018
temp = -1/(*(dnode++)+4);
1019
*dnode = temp;
1020
temp *= *(u++) - 6*(*(nodep+1)+*(nodep-1)-2**nodep);
1021
*u = temp;
1022
}
1023
*(++dnode) = 0.0; /* "natural" upper boundary condition */
1024
for (x=nbx-2; x--;)
1025
{
1026
temp = *(dnode--);
1027
*dnode = (*dnode*temp+*(u--))/6.0;
1028
}
1029
free(u);
1030
dnode--;
1031
}
1032
}
1033
else
1034
{
1035
/*-- No interpolation and no new 2nd derivatives needed along y */
1036
node = field->sigma;
1037
dnode = field->dsigma;
1038
}
1039
1040
/*-- Interpolation along x */
1041
width = field->width;
1042
if (nbx>1)
1043
{
1044
nx = field->backw;
1045
xstep = 1.0/nx;
1046
changepoint = nx/2;
1047
dx = (xstep - 1)/2; /* dx of the first pixel in the row */
1048
dx0 = ((nx+1)%2)*xstep/2; /* dx of the 1st pixel right to a bkgnd node */
1049
blo = node;
1050
bhi = node + 1;
1051
dblo = dnode;
1052
dbhi = dnode + 1;
1053
for (x=i=0,j=width; j--; i++, dx += xstep)
1054
{
1055
if (i==changepoint && x>0 && x<nbxm1)
1056
{
1057
blo++;
1058
bhi++;
1059
dblo++;
1060
dbhi++;
1061
dx = dx0;
1062
}
1063
cdx = 1 - dx;
1064
*(line++) = (PIXTYPE)(cdx*(*blo+(cdx*cdx-1)**dblo)
1065
+ dx*(*bhi+(dx*dx-1)**dbhi));
1066
if (i==nx)
1067
{
1068
x++;
1069
i = 0;
1070
}
1071
}
1072
}
1073
else
1074
for (j=width; j--;)
1075
*(line++) = (PIXTYPE)*node;
1076
1077
if (nby>1)
1078
{
1079
free(node);
1080
free(dnode);
1081
}
1082
1083
return;
1084
}
1085
1086
1087
/********************************* end_back **********************************/
1088
/*
1089
Terminate background procedures (mainly freeing memory).
1090
*/
1091
void end_back(fieldstruct *field)
1092
1093
{
1094
free(field->back);
1095
free(field->dback);
1096
free(field->sigma);
1097
free(field->dsigma);
1098
free(field->backline);
1099
1100
return;
1101
}
1102
Loggerhead is a web-based interface for Breezy
Version: 2.0.1