public software.psfex

[/] [branches/] [rhl/] [src/] [sample_utils.c] - Blame information for rev 191

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 191 rhl
/*
2
*                               utils.c
3
*
4
* Miscellaneous functions.
5
*
6
*
7
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8
*
9
*       This file part of:      PSFEx
10
*
11
*       Copyright:              (C) 1997-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
12
*
13
*       License:                GNU General Public License
14
*
15
*       PSFEx is free software: you can redistribute it and/or modify
16
*       it under the terms of the GNU General Public License as published by
17
*       the Free Software Foundation, either version 3 of the License, or
18
*       (at your option) any later version.
19
*       PSFEx is distributed in the hope that it will be useful,
20
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
21
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22
*       GNU General Public License for more details.
23
*       You should have received a copy of the GNU General Public License
24
*       along with PSFEx.  If not, see <http://www.gnu.org/licenses/>.
25
*
26
*       Last modified:          14/07/2013, Bastille Day
27
*
28
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
29
 
30
#ifdef HAVE_CONFIG_H
31
#include        "config.h"
32
#endif
33
 
34
#include        <assert.h>
35
#include        <ctype.h>
36
#include        <math.h>
37
#include        <stdio.h>
38
#include        <stdlib.h>
39
#include        <string.h>
40
#include        <time.h>
41
 
42
#include        "define.h"
43
#include        "globals.h"
44
#include        "prefs.h"
45
 
46
/*****************************************************************************/
47
/*
48
 * Copied here from sample.c as the rest of that file's all about fits I/O
49
 */
50
#include "sample.h"
51
 
52
/****** malloc_samples *******************************************************
53
PROTO   void malloc_samples(setstruct *set, int nsample)
54
PURPOSE Allocate memory for a set of samples.
55
INPUT   set structure pointer,
56
        desired number of samples.
57
OUTPUT  -.
58
NOTES   -.
59
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
60
VERSION 26/03/2008
61
*/
62
void    malloc_samples(setstruct *set, int nsample)
63
 
64
  {
65
   samplestruct *sample;
66
   int          n;
67
 
68
  QMALLOC(set->sample, samplestruct *, nsample);
69
  set->nsamplemax = nsample;
70
  if (set != set->samples_owner)
71
    return;
72
 
73
  QMALLOC(set->s_sample, samplestruct, nsample);
74
  for (n=0; n!=nsample; ++n)
75
    {
76
    sample = set->sample[n] = &set->s_sample[n];
77
    QMALLOC(sample->vig, float, set->nvig);
78
    QMALLOC(sample->vigresi, float, set->nvig);
79
    QMALLOC(sample->vigweight, float, set->nvig);
80
    QMALLOC(sample->vigchi, float, set->nvig);
81
    if (set->ncontext)
82
      QMALLOC(sample->context, double, set->ncontext);
83
    }
84
 
85
  return;
86
  }
87
 
88
 
89
/****** realloc_samples ******************************************************
90
PROTO   void realloc_samples(setstruct *set, int nsample)
91
PURPOSE Re-allocate memory for a set of samples.
92
INPUT   set structure pointer,
93
        desired number of samples.
94
OUTPUT  -.
95
NOTES   -.
96
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
97
VERSION 26/03/2008
98
*/
99
void    realloc_samples(setstruct *set, int nsample)
100
 
101
  {
102
   samplestruct *sample;
103
   int          n;
104
 
105
  if (set != set->samples_owner)
106
    {
107
    assert(nsample <= set->nsamplemax);
108
    return;
109
    }
110
 
111
/* If we want to reallocate 0 samples, better free the whole thing! */
112
  if (!nsample)
113
    free_samples(set);
114
 
115
/* Two cases: either more samples are required, or the opposite! */
116
  if (nsample>set->nsamplemax)
117
    {
118
    QREALLOC(set->s_sample, samplestruct, nsample);
119
    free(set->sample);
120
    QMALLOC(set->sample, samplestruct *, nsample);
121
    for (n = set->nsamplemax; n<nsample; n++)
122
      {
123
      sample = &set->s_sample[n];
124
      QMALLOC(sample->vig, float, set->nvig);
125
      QMALLOC(sample->vigresi, float, set->nvig);
126
      QMALLOC(sample->vigchi, float, set->nvig);
127
      QMALLOC(sample->vigweight, float, set->nvig);
128
      if (set->ncontext)
129
        QMALLOC(sample->context, double, set->ncontext);
130
      }
131
    }
132
  else if (nsample<set->nsamplemax)
133
    {
134
    for (n = nsample; n<set->nsamplemax; n++)
135
      {
136
      sample = &set->s_sample[n];
137
      free(sample->vig);
138
      free(sample->vigresi);
139
      free(sample->vigchi);
140
      free(sample->vigweight);
141
      if (set->ncontext)
142
        free(sample->context);
143
      }
144
    QREALLOC(set->s_sample, samplestruct, nsample);
145
    free(set->sample);
146
    QMALLOC(set->sample, samplestruct *, nsample);
147
    }
148
 
149
  set->nsamplemax = nsample;
150
 
151
  for (n = 0; n<nsample; n++)
152
     set->sample[n] = &set->s_sample[n];
153
 
154
  return;
155
  }
156
 
157
 
158
/****** free_samples *********************************************************
159
PROTO   void free_samples(setstruct *set, int nsample)
160
PURPOSE free memory for a set of samples.
161
INPUT   set structure pointer,
162
        desired number of samples.
163
OUTPUT  -.
164
NOTES   -.
165
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
166
VERSION 26/03/2008
167
*/
168
void    free_samples(setstruct *set)
169
 
170
  {
171
   samplestruct *sample;
172
   int          n;
173
 
174
   if (set->samples_owner == set)
175
     {
176
     for (n = 0; n<set->nsamplemax; ++n )
177
       {
178
       sample = set->sample[n];
179
       free(sample->vig);
180
       free(sample->vigresi);
181
       free(sample->vigweight);
182
       free(sample->vigchi);
183
       if (set->ncontext)
184
         free(sample->context);
185
       }
186
     free(set->s_sample);
187
     }
188
 
189
  free(set->sample);
190
  set->sample = NULL;
191
  set->nsample = set->nsamplemax = 0;
192
 
193
  return;
194
  }
195
 
196
 
197
/****** remove_sample ********************************************************
198
PROTO   samplestruct *remove_sample(setstruct *set, int isample)
199
PURPOSE Remove an element from a set of samples.
200
INPUT   set structure pointer,
201
        sample number.
202
OUTPUT  The new pointer for the element that replaced the removed one.
203
NOTES   -.
204
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
205
VERSION 01/03/99
206
*/
207
samplestruct    *remove_sample(setstruct *set, int isample)
208
 
209
  {
210
   static samplestruct  exsample;
211
   samplestruct         *sample;
212
   int                  nsample;
213
 
214
/* If we want to reallocate 0 samples, better free the whole thing! */
215
  nsample = set->nsample-1;
216
  if (nsample>0)
217
    {
218
    sample = set->sample[isample];
219
    exsample = *set->sample[nsample];
220
    *set->sample[nsample] = *sample;
221
    *sample = exsample;
222
    }
223
   else
224
     nsample=0;
225
  realloc_samples(set, nsample);
226
  set->nsample = nsample;
227
 
228
  return set->sample[isample];
229
  }
230
 
231
 
232
/****** init_set ************************************************************
233
PROTO   setstruct *init_set()
234
PURPOSE Allocate and initialize a set structure.
235
INPUT   -.
236
OUTPUT  -.
237
NOTES   See prefs.h.
238
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
239
VERSION 26/03/2008
240
*/
241
setstruct       *init_set(contextstruct *context)
242
 
243
  {
244
   setstruct    *set;
245
   int          i;
246
 
247
  QCALLOC(set, setstruct, 1);
248
  set->nsample = set->nsamplemax = 0;
249
  set->samples_owner = set;             /* we'll own the samples, so we'll need to free them */
250
  set->vigdim = 2;
251
  QMALLOC(set->vigsize, int, set->vigdim);
252
  set->ncontext = context->ncontext;
253
  if (set->ncontext)
254
    {
255
    QMALLOC(set->contextoffset, double, set->ncontext);
256
    QMALLOC(set->contextscale, double, set->ncontext);
257
    QMALLOC(set->contextname, char *, set->ncontext);
258
    for (i=0; i<set->ncontext; i++)
259
      QMALLOC(set->contextname[i], char, 80);
260
    }
261
 
262
  return set;
263
  }
264
 
265
 
266
/****** end_set *************************************************************
267
PROTO   void end_set(setstruct *set)
268
PURPOSE free memory allocated by a complete set structure.
269
INPUT   set structure pointer,
270
OUTPUT  -.
271
NOTES   -.
272
AUTHOR  E. Bertin (IAP, Leiden observatory & ESO)
273
VERSION 26/03/2008
274
*/
275
void    end_set(setstruct *set)
276
 
277
  {
278
   int  i;
279
 
280
  free_samples(set);
281
  free(set->vigsize);
282
  if (set->ncontext)
283
    {
284
    for (i=0; i<set->ncontext; i++)
285
      free(set->contextname[i]);
286
    free(set->contextname);
287
    free(set->contextoffset);
288
    free(set->contextscale);
289
    }
290
  free(set->head);
291
  free(set);
292
 
293
  return;
294
  }
295
 
296
 
297
/****** make_weights *********************************************************
298
PROTO   void make_weights(setstruct *set, samplestruct *sample)
299
PURPOSE Produce a weight-map for each sample vignet.
300
INPUT   set structure pointer,
301
        sample structure pointer.
302
OUTPUT  -.
303
NOTES   -.
304
AUTHOR  E. Bertin (IAP,Leiden observatory & ESO)
305
VERSION 13/08/2007
306
*/
307
void make_weights(setstruct *set, samplestruct *sample)
308
 
309
  {
310
   float        *vig, *vigweight,
311
                backnoise2, gain, noise2, profaccu2, pix;
312
   int          i;
313
 
314
/* Produce a weight-map */
315
  profaccu2 = prefs.prof_accuracy*prefs.prof_accuracy;
316
  gain = sample->gain;
317
  backnoise2 = sample->backnoise2;
318
  for (vig=sample->vig, vigweight=sample->vigweight, i=set->nvig; i--;)
319
    {
320
    if (*vig <= -BIG)
321
      *(vig++) = *(vigweight++) = 0.0;
322
    else
323
      {
324
      pix = *(vig++);
325
      noise2 = backnoise2 + profaccu2*pix*pix;
326
      if (pix>0.0 && gain>0.0)
327
        noise2 += pix/gain;
328
      *(vigweight++) = 1.0/noise2;
329
      }
330
    }
331
 
332
  return;
333
  }
334
 
335
 
336
/****** recenter_sample ******************************************************
337
PROTO   void recenter_samples(samplestruct sample,
338
                setstruct *set, float fluxrad)
339
PURPOSE Recenter sample image using windowed barycenter.
340
INPUT   sample structure pointer,
341
        set structure pointer,
342
        flux radius.
343
OUTPUT  -.
344
NOTES   -.
345
AUTHOR  E. Bertin (IAP)
346
VERSION 16/11/2009
347
*/
348
void    recenter_sample(samplestruct *sample, setstruct *set, float fluxrad)
349
 
350
  {
351
   double       tv, dxpos, dypos;
352
   float        *ima,*imat,*weight,*weightt,
353
                pix, var, locpix, locarea, sig,twosig2, raper,raper2,
354
                offsetx,offsety, mx,my, mx2ph,my2ph,
355
                rintlim,rintlim2,rextlim2, scalex,scaley,scale2,
356
                dx,dy, dx1,dy2, r2;
357
   int          i, x,y,x2,y2, w,h, sx,sy, xmin,xmax,ymin,ymax, pos;
358
 
359
  sig = fluxrad*2.0/2.35; /* From half-FWHM to sigma */
360
  twosig2 = 2.0*sig*sig;
361
 
362
  w = set->vigsize[0];
363
  h = set->vigsize[1];
364
/* Integration radius */
365
  raper = RECENTER_NSIG*sig;
366
  raper2 = raper*raper;
367
/* Internal radius of the oversampled annulus (<r-sqrt(2)/2) */
368
  rintlim = raper - 0.75;
369
  rintlim2 = (rintlim>0.0)? rintlim*rintlim: 0.0;
370
/* External radius of the oversampled annulus (>r+sqrt(2)/2) */
371
  rextlim2 = (raper + 0.75)*(raper + 0.75);
372
  scaley = scalex = 1.0/RECENTER_OVERSAMP;
373
  scale2 = scalex*scaley;
374
  offsetx = 0.5*(scalex-1.0);
375
  offsety = 0.5*(scaley-1.0);
376
/* Use isophotal centroid as a first guess */
377
  mx = sample->dx + (float)(w/2);
378
  my = sample->dy + (float)(h/2);
379
 
380
  for (i=0; i<RECENTER_NITERMAX; i++)
381
    {
382
    xmin = (int)(mx-raper+0.499999);
383
    xmax = (int)(mx+raper+1.499999);
384
    ymin = (int)(my-raper+0.499999);
385
    ymax = (int)(my+raper+1.499999);
386
    mx2ph = mx*2.0 + 0.49999;
387
    my2ph = my*2.0 + 0.49999;
388
 
389
    if (xmin < 0)
390
      xmin = 0;
391
    if (xmax > w)
392
      xmax = w;
393
    if (ymin < 0)
394
      ymin = 0;
395
    if (ymax > h)
396
      ymax = h;
397
    tv = 0.0;
398
    dxpos = dypos = 0.0;
399
    ima = sample->vig;
400
    weight = sample->vigweight;
401
    for (y=ymin; y<ymax; y++)
402
      {
403
      imat= ima + (pos = y*w + xmin);
404
      weightt = weight + pos;
405
      dy = y - my;
406
      for (x=xmin; x<xmax; x++, imat++, weightt++)
407
        {
408
        dx = x - mx;
409
        if ((r2=dx*dx+dy*dy)<rextlim2)
410
          {
411
          if (RECENTER_OVERSAMP>1 && r2> rintlim2)
412
            {
413
            dx += offsetx;
414
            dy += offsety;
415
            locarea = 0.0;
416
            for (sy=RECENTER_OVERSAMP; sy--; dy+=scaley)
417
              {
418
              dx1 = dx;
419
              dy2 = dy*dy;
420
              for (sx=RECENTER_OVERSAMP; sx--; dx1+=scalex)
421
                if (dx1*dx1+dy2<raper2)
422
                  locarea += scale2;
423
              }
424
            }
425
          else
426
            locarea = 1.0;
427
          locarea *= expf(-r2/twosig2);
428
/*-------- Here begin tests for pixel and/or weight overflows. Things are a */
429
/*-------- bit intricated to have it running as fast as possible in the most */
430
/*-------- common cases */
431
          pix = *imat;
432
          if ((var=*weightt)==0.0)
433
            {
434
            if ((x2=(int)(mx2ph-x))>=0 && x2<w
435
                && (y2=(int)(my2ph-y))>=0 && y2<h
436
                && (var=*(weight + (pos = y2*w + x2)))>0.0)
437
              {
438
              var = 1.0/var;
439
              pix = *(ima+pos);
440
              }
441
            else
442
              pix = var = 0.0;
443
            }
444
          dx = x - mx;
445
          dy = y - my;
446
          locpix = locarea*pix;
447
          tv += locpix;
448
          dxpos += locpix*dx;
449
          dypos += locpix*dy;
450
          }
451
        }
452
      }
453
 
454
    if (tv>0.0)
455
      {
456
      mx += (dxpos /= tv)*RECENTER_GRADFAC;
457
      my += (dypos /= tv)*RECENTER_GRADFAC;
458
      }
459
    else
460
      break;
461
 
462
/*-- Stop here if position does not change */
463
    if (dxpos*dxpos+dypos*dypos < RECENTER_STEPMIN*RECENTER_STEPMIN)
464
      break;
465
    }
466
 
467
  sample->dx = mx - (float)(w/2);
468
  sample->dy = my - (float)(h/2);
469
 
470
  return;
471
  }