← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/feros/libsrc/fit_back.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/feros/libsrc/fit_back.c
 
1
/* @(#)fit_back.c       19.1 (ESO-IPG) 02/25/03 14:21:24 */
 
2
/* 
 
3
 
 
4
   Otmar Stahl
 
5
 
 
6
   fit_back.c  
 
7
 
 
8
   interpolate background 
 
9
 
 
10
   1) compute median between orders
 
11
   2) smoothing splines along lines, compute background values at grid-points
 
12
   3) bi-cubic splines over the grid-points
 
13
   4) evaluate splines at all points
 
14
 
 
15
*/
 
16
 
 
17
/* general Midas includes */
 
18
 
 
19
#include <midas_def.h>
 
20
#include <tbldef.h>
 
21
 
 
22
/* FEROS specific includes */
 
23
 
 
24
#include <glsp.h>
 
25
#include <proto_nrutil.h>
 
26
#include <proto_mutil.h>
 
27
#include <echelle.h>
 
28
 
 
29
int fit_back
 
30
#ifdef __STDC__
 
31
(
 
32
 float xval[], float yval[], float **ya, float **y2a, 
 
33
 int npix[], int imno, int m, int n, 
 
34
 int nact, float **centers, int xysize[], int bgdist, int fibmode
 
35
 )
 
36
#else
 
37
     (
 
38
      xval, yval, ya, y2a, npix, imno, m, n, nact, 
 
39
      centers, xysize, bgdist, fibmode
 
40
      )
 
41
     int npix[], imno, nact, m, n, xysize[], bgdist, fibmode; 
 
42
     float **ya, **y2a, xval[], yval[], **centers;
 
43
#endif
 
44
{
 
45
  int i, ii, j, jj, k, iarr1, iy, actvals, center;
 
46
  int splwidth[2], splwin[2];
 
47
  float *arr1, *xpos, *ypos, *imagex1;
 
48
  double xxd, *xn, *fn, *w, *a, *b, *c, *d, ausg[3];
 
49
  double xposs, weight;
 
50
  int status;
 
51
 
 
52
  weight = 1.0e-06; /* smoothing parameter, this is a magic number ? */  
 
53
 
 
54
  splwin[0] = xysize[0];     
 
55
  splwin[1] = xysize[1];
 
56
  splwidth[0] = splwin[0] * 2 + 1;
 
57
  splwidth[1] = splwin[0] * 2 + 1;
 
58
 
 
59
  xn = dvector(0, 2*nact + 2);                  /* alloc mem for vars */
 
60
  fn = dvector(0, 2*nact + 2);
 
61
  w = dvector(0, 2*nact + 2);
 
62
  a = dvector(0, 2*nact + 2);
 
63
  b = dvector(0, 2*nact + 2);
 
64
  c = dvector(0, 2*nact + 2);
 
65
  d = dvector(0, 2*nact + 2);
 
66
 
 
67
  imagex1 = vector(0, npix[0] * splwidth[1]);
 
68
 
 
69
  arr1 = vector (0, splwidth[0] * splwidth[1] + 1);
 
70
  xpos = vector (1, 2*nact + 1);
 
71
  ypos = vector (1, 2*nact + 1);
 
72
 
 
73
  /* init various arrays */
 
74
 
 
75
  for(k = 1; k <= n; k++ )                         /* for each row */
 
76
    {
 
77
      j = (int) yval[k];
 
78
 
 
79
      /* read one line  */
 
80
 
 
81
      SCFGET(imno, (j - splwin[1]) * npix[0] + 1, npix[0] * splwidth[1], 
 
82
             &actvals,  (char *) imagex1); /* data are at *imagex1 */
 
83
 
 
84
      /* loop over the orders */
 
85
 
 
86
      ii = 0;
 
87
      for (i = 1; i <= nact; i++)           /* for all orders  */
 
88
        {
 
89
 
 
90
          /* between two orders if sufficient space */
 
91
 
 
92
          if ( (i > 1) && (centers[i][j] - centers[i-1][j]) > bgdist)
 
93
            {
 
94
              xposs = (centers[i][j] + centers[i-1][j]) / 2;
 
95
              center = (int) xposs;
 
96
              
 
97
              if (center > splwin[0] && center < npix[0] - splwin[0])
 
98
                {
 
99
                  ii++;
 
100
                  iarr1 = 0;
 
101
                  for(jj = -splwin[0]; jj <= splwin[0]; jj++)
 
102
                    {
 
103
                      for (iy = 0; iy < splwidth[1]; iy++)
 
104
                        arr1[++iarr1] = imagex1[center + jj + iy * npix[0]];
 
105
                      
 
106
                    }
 
107
                  xpos[ii] = xposs;
 
108
                  ypos[ii] = select_pos((int) iarr1/2 , iarr1 , arr1);
 
109
                }
 
110
            }
 
111
 
 
112
          /* between the fibers for fiber mode 2 */
 
113
          
 
114
          if (fibmode > 1)
 
115
            {
 
116
              xposs = (centers[i][j]);
 
117
              
 
118
              center = (int) xposs;
 
119
              
 
120
              /* 
 
121
                 populate array for median filtering,
 
122
                 this needs more work, e.g.
 
123
                 median over a number of lines
 
124
              */ 
 
125
              
 
126
              if (center > splwin[0] && center < npix[0] - splwin[0])
 
127
                {
 
128
                  ii++;
 
129
                  iarr1 = 0;
 
130
                  for(jj = -splwin[0]; jj <= splwin[0]; jj++)
 
131
                    {
 
132
                      for (iy = 0; iy < splwidth[1]; iy++)
 
133
                        arr1[++iarr1] = imagex1[center + jj + iy * npix[0]];
 
134
                      
 
135
                    }
 
136
                  xpos[ii] = xposs;
 
137
                  ypos[ii] = select_pos((int) iarr1/2 , iarr1 , arr1);
 
138
                }
 
139
            }
 
140
        }
 
141
      
 
142
      for(i = 0; i < ii; i++)
 
143
        {
 
144
          xn[i] = (double) xpos[i + 1];
 
145
          fn[i] = (double) ypos[i + 1];
 
146
          w[i] = weight;
 
147
        }
 
148
 
 
149
      /* smoothing spline in X-direction and compute new grid:
 
150
         see Engeln-Muellges & Reuther, 1990; BI Wissenschaftsverlag Mannheim;
 
151
         Formelsammlung zur numerischen Mathematik mit C-Programmen, page 236f
 
152
      */
 
153
    
 
154
      status = glspnp(ii - 1, xn, fn, w, 2, 0.0, 0.0, a, b, c, d);
 
155
 
 
156
      for(i = 1; i <= m; i++)
 
157
        {
 
158
          xxd =  spval(ii - 1, (double) xval[i], a, b, c, d, xn, ausg);
 
159
          ya[i][k] = (float) xxd;
 
160
 
 
161
          /* smoothed function values at x[i = 1 to m] are in ya[i][k] for any row 
 
162
             k (= 1 .. n) */
 
163
        }
 
164
    }
 
165
 
 
166
  /* to compute second derivatives, see Recipes, pg. 128 - we need
 
167
     an rectangular grid.*/
 
168
 
 
169
  splie2(xval, yval, ya, m, n, y2a); 
 
170
  
 
171
  free_dvector(xn, 0, 2*nact + 2);
 
172
  free_dvector(fn, 0, 2*nact + 2);
 
173
  free_dvector(w, 0, 2*nact + 2);
 
174
  free_dvector(a, 0, 2*nact + 2);
 
175
  free_dvector(b, 0, 2*nact + 2);
 
176
  free_dvector(c, 0, 2*nact + 2);
 
177
  free_dvector(d, 0, 2*nact + 2);
 
178
 
 
179
  free_vector(imagex1, 0, npix[0] * splwidth[1]);
 
180
  free_vector(arr1, 0, (splwidth[0] * splwidth[1] + 1));
 
181
  free_vector(xpos, 1,2*nact + 1);
 
182
  free_vector(ypos, 1, 2*nact + 1);
 
183
 
 
184
  return 0;
 
185
}