← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/general/libsrc/lin.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/libsrc/lin.c
 
1
/*=============================================================================
 
2
*
 
3
*   WCSLIB - an implementation of the FITS WCS proposal.
 
4
*   Copyright (C) 1995, Mark Calabretta
 
5
*
 
6
*   This library is free software; you can redistribute it and/or modify it
 
7
*   under the terms of the GNU Library General Public License as published
 
8
*   by the Free Software Foundation; either version 2 of the License, or (at
 
9
*   your option) any later version.
 
10
*
 
11
*   This library is distributed in the hope that it will be useful, but
 
12
*   WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library
 
14
*   General Public License for more details.
 
15
*
 
16
*   You should have received a copy of the GNU Library General Public License
 
17
*   along with this library; if not, write to the Free Software Foundation,
 
18
*   Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
19
*
 
20
*   Correspondence concerning WCSLIB may be directed to:
 
21
*      Internet email: mcalabre@atnf.csiro.au
 
22
*      Postal address: Dr. Mark Calabretta,
 
23
*                      Australia Telescope National Facility,
 
24
*                      P.O. Box 76,
 
25
*                      Epping, NSW, 2121,
 
26
*                      AUSTRALIA
 
27
*
 
28
*=============================================================================
 
29
*
 
30
*   C routines which implement the FITS World Coordinate System (WCS)
 
31
*   convention.
 
32
*
 
33
*   Summary of routines
 
34
*   -------------------
 
35
*   These utility routines apply the linear transformation defined by the WCS
 
36
*   FITS header cards.  There are separate routines for the image-to-pixel,
 
37
*   linfwd(), and pixel-to-image, linrev(), transformations.
 
38
*
 
39
*   An initialization routine, linset(), computes intermediate values from
 
40
*   the transformation parameters but need not be called explicitly - see the
 
41
*   explanation of lin.flag below.
 
42
*
 
43
*   An auxiliary matrix inversion routine, matinv(), is included.  It uses
 
44
*   LU-triangular factorization with scaled partial pivoting.
 
45
*
 
46
*
 
47
*   Initialization routine; linset()
 
48
*   --------------------------------
 
49
*   Initializes members of a linprm data structure which hold intermediate
 
50
*   values.  Note that this routine need not be called directly; it will be
 
51
*   invoked by linfwd() and linrev() if the "flag" structure member is
 
52
*   anything other than a predefined magic value.
 
53
*
 
54
*   Given and/or returned:
 
55
*      lin      linprm*  Linear transformation parameters (see below).
 
56
*
 
57
*   Function return value:
 
58
*               int      Error status
 
59
*                           0: Success.
 
60
*                           1: Memory allocation error.
 
61
*                           2: PC matrix is singular.
 
62
*
 
63
*   Forward transformation; linfwd()
 
64
*   --------------------------------
 
65
*   Compute pixel coordinates from image coordinates.  Note that where
 
66
*   celestial coordinate systems are concerned the image coordinates
 
67
*   correspond to (x,y) in the plane of projection, not celestial (lng,lat).
 
68
*
 
69
*   Given:
 
70
*      imgcrd   d[]      Image (world) coordinate.
 
71
*
 
72
*   Given and returned:
 
73
*      lin      linprm*  Linear transformation parameters (see below).
 
74
*
 
75
*   Returned:
 
76
*      pixcrd   d[]      Pixel coordinate.
 
77
*
 
78
*   Function return value:
 
79
*               int      Error status
 
80
*                           0: Success.
 
81
*                           1: The transformation is not invertible.
 
82
*
 
83
*   Reverse transformation; linrev()
 
84
*   --------------------------------
 
85
*   Compute image coordinates from pixel coordinates.  Note that where
 
86
*   celestial coordinate systems are concerned the image coordinates
 
87
*   correspond to (x,y) in the plane of projection, not celestial (lng,lat).
 
88
*
 
89
*   Given:
 
90
*      pixcrd   d[]      Pixel coordinate.
 
91
*
 
92
*   Given and/or returned:
 
93
*      lin      linprm*  Linear transformation parameters (see below).
 
94
*
 
95
*   Returned:
 
96
*      imgcrd   d[]      Image (world) coordinate.
 
97
*
 
98
*   Function return value:
 
99
*               int      Error status
 
100
*                           0: Success.
 
101
*                           1: Error.
 
102
*
 
103
*   Linear transformation parameters
 
104
*   --------------------------------
 
105
*   The linprm struct consists of the following:
 
106
*
 
107
*      int flag
 
108
*         This flag must be set to zero whenever any of the following members
 
109
*         are set or modified.  This signals the initialization routine,
 
110
*         linset(), to recompute intermediaries.
 
111
*      int naxis
 
112
*         Number of image axes.
 
113
*      double *crpix
 
114
*         Pointer to the first element of an array of double containing the
 
115
*         coordinate reference pixel, CRPIXn.
 
116
*      double *pc
 
117
*         Pointer to the first element of the PC (pixel coordinate)
 
118
*         transformation matrix.
 
119
*      double *cdelt
 
120
*         Pointer to the first element of an array of double containing the
 
121
*         coordinate increments, CDELTn.
 
122
*
 
123
*   The remaining members of the linprm struct are maintained by the
 
124
*   initialization routine and should not be modified.  Storage will be
 
125
*   allocated for these arrays via malloc().
 
126
*
 
127
*      double *piximg
 
128
*         Pointer to the first element of the matrix containing the product
 
129
*         of the CDELTn diagonal matrix and the PC matrix.
 
130
*      double *imgpix
 
131
*         Pointer to the first element of the inverse of the piximg matrix.
 
132
*
 
133
*   Author: Mark Calabretta, Australia Telescope National Facility
 
134
 
 
135
 
 
136
.VERSION
 
137
 
 
138
 090701         last modif
 
139
*===========================================================================*/
 
140
 
 
141
#include <stdlib.h>
 
142
#include <math.h>
 
143
#include <lin.h>
 
144
 
 
145
int linset(lin)
 
146
 
 
147
struct linprm *lin;
 
148
 
 
149
{
 
150
   int i, ij, j, mem, n;
 
151
 
 
152
 
 
153
 
 
154
   n = lin->naxis;
 
155
 
 
156
   /* Allocate memory for internal arrays. */
 
157
   mem = n * n * sizeof(double);
 
158
   lin->piximg = (double*)malloc((size_t)mem);
 
159
   if (lin->piximg == (double*)0) return 1;
 
160
 
 
161
   lin->imgpix = (double*)malloc((size_t)mem);
 
162
   if (lin->imgpix == (double*)0) 
 
163
      {
 
164
      (void) free(lin->piximg);
 
165
      return 1;
 
166
      }
 
167
 
 
168
   /* Compute the pixel-to-image transformation matrix. */
 
169
   for (i = 0, ij = 0; i < n; i++) {
 
170
      for (j = 0; j < n; j++, ij++) {
 
171
         lin->piximg[ij] = lin->cdelt[i] * lin->pc[ij];
 
172
      }
 
173
   }
 
174
 
 
175
   /* Compute the image-to-pixel transformation matrix. */
 
176
   if (matinv(n, lin->piximg, lin->imgpix)) 
 
177
      {
 
178
      (void) free(lin->piximg);
 
179
      (void) free(lin->imgpix);
 
180
      return 2;
 
181
      }
 
182
 
 
183
   lin->flag = LINSET;
 
184
 
 
185
   return 0;
 
186
}
 
187
 
 
188
/*--------------------------------------------------------------------------*/
 
189
 
 
190
int linfwd(imgcrd, lin, pixcrd)
 
191
 
 
192
double imgcrd[];
 
193
struct linprm *lin;
 
194
double pixcrd[];
 
195
 
 
196
{
 
197
   int i, ij, j, n;
 
198
 
 
199
   n = lin->naxis;
 
200
 
 
201
   if (lin->flag != LINSET) {
 
202
      if (linset(lin)) return 1;
 
203
   }
 
204
 
 
205
   for (i = 0, ij = 0; i < n; i++) {
 
206
      pixcrd[i] = 0.0;
 
207
      for (j = 0; j < n; j++, ij++) {
 
208
         pixcrd[i] += lin->imgpix[ij] * imgcrd[j];
 
209
      }
 
210
   }
 
211
 
 
212
   for (j = 0; j < n; j++) {
 
213
      pixcrd[j] += lin->crpix[j];
 
214
   }
 
215
 
 
216
   return 0;
 
217
}
 
218
 
 
219
/*--------------------------------------------------------------------------*/
 
220
 
 
221
int linrev(pixcrd, lin, imgcrd)
 
222
 
 
223
double pixcrd[];
 
224
struct linprm *lin;
 
225
double imgcrd[];
 
226
 
 
227
{
 
228
   int i, ij, j, n;
 
229
   double temp;
 
230
 
 
231
   n = lin->naxis;
 
232
 
 
233
   if (lin->flag != LINSET) {
 
234
      if (linset(lin)) return 1;
 
235
   }
 
236
 
 
237
   for (i = 0; i < n; i++) {
 
238
      imgcrd[i] = 0.0;
 
239
   }
 
240
 
 
241
   for (j = 0; j < n; j++) {
 
242
      temp = pixcrd[j] - lin->crpix[j];
 
243
      for (i = 0, ij = j; i < n; i++, ij+=n) {
 
244
         imgcrd[i] += lin->piximg[ij] * temp;
 
245
      }
 
246
   }
 
247
 
 
248
   return 0;
 
249
}
 
250
 
 
251
/*--------------------------------------------------------------------------*/
 
252
 
 
253
int matinv(n, mat, inv)
 
254
 
 
255
int n;
 
256
double mat[], inv[];
 
257
 
 
258
{
 
259
   register int i, ij, ik, j, k, kj, pj;
 
260
   int    status, itemp, mem, *mxl, *lxm, pivot;
 
261
 
 
262
   double colmax, *lu, *rowmax, dtemp;
 
263
 
 
264
 
 
265
 
 
266
   lu = (double *) 0;
 
267
   rowmax = (double *) 0;
 
268
 
 
269
   /* Allocate memory for internal arrays. */
 
270
   status = 1;
 
271
   mem = n * sizeof(int);
 
272
   if ((mxl = (int*)malloc((size_t)mem)) == (int*)0) return status;
 
273
   if ((lxm = (int*)malloc((size_t)mem)) == (int*)0) goto end_of_it;
 
274
 
 
275
   mem = n * sizeof(double);
 
276
   if ((rowmax = (double*)malloc((size_t)mem)) == (double*)0) goto end_of_it;
 
277
 
 
278
   mem *= n;
 
279
   if ((lu = (double*)malloc((size_t)mem)) == (double*)0) goto end_of_it;
 
280
 
 
281
 
 
282
   /* Initialize arrays. */
 
283
   for (i = 0, ij = 0; i < n; i++) {
 
284
      /* Vector which records row interchanges. */
 
285
      mxl[i] = i;
 
286
 
 
287
      rowmax[i] = 0.0;
 
288
 
 
289
      for (j = 0; j < n; j++, ij++) {
 
290
         dtemp = fabs(mat[ij]);
 
291
         if (dtemp > rowmax[i]) rowmax[i] = dtemp;
 
292
 
 
293
         lu[ij] = mat[ij];
 
294
      }
 
295
 
 
296
      /* A row of zeroes indicates a singular matrix. */
 
297
      if (rowmax[i] == 0.0) 
 
298
         {
 
299
         status = 2;
 
300
         goto end_of_it;
 
301
         }
 
302
   }
 
303
 
 
304
 
 
305
   /* Form the LU triangular factorization using scaled partial pivoting. */
 
306
   for (k = 0; k < n; k++) {
 
307
      /* Decide whether to pivot. */
 
308
      colmax = fabs(lu[k*n+k]) / rowmax[k];
 
309
      pivot = k;
 
310
 
 
311
      for (i = k+1; i < n; i++) {
 
312
         ik = i*n + k;
 
313
         dtemp = fabs(lu[ik]) / rowmax[i];
 
314
         if (dtemp > colmax) {
 
315
            colmax = dtemp;
 
316
            pivot = i;
 
317
         }
 
318
      }
 
319
 
 
320
      if (pivot > k) {
 
321
         /* We must pivot, interchange the rows of the design matrix. */
 
322
         for (j = 0, pj = pivot*n, kj = k*n; j < n; j++, pj++, kj++) {
 
323
            dtemp = lu[pj];
 
324
            lu[pj] = lu[kj];
 
325
            lu[kj] = dtemp;
 
326
         }
 
327
 
 
328
         /* Amend the vector of row maxima. */
 
329
         dtemp = rowmax[pivot];
 
330
         rowmax[pivot] = rowmax[k];
 
331
         rowmax[k] = dtemp;
 
332
 
 
333
         /* Record the interchange for later use. */
 
334
         itemp = mxl[pivot];
 
335
         mxl[pivot] = mxl[k];
 
336
         mxl[k] = itemp;
 
337
      }
 
338
 
 
339
      /* Gaussian elimination. */
 
340
      for (i = k+1; i < n; i++) {
 
341
         ik = i*n + k;
 
342
 
 
343
         /* Nothing to do if lu[ik] is zero. */
 
344
         if (lu[ik] != 0.0) {
 
345
            /* Save the scaling factor. */
 
346
            lu[ik] /= lu[k*n+k];
 
347
 
 
348
            /* Subtract rows. */
 
349
            for (j = k+1; j < n; j++) {
 
350
               lu[i*n+j] -= lu[ik]*lu[k*n+j];
 
351
            }
 
352
         }
 
353
      }
 
354
   }
 
355
 
 
356
 
 
357
   /* mxl[i] records which row of mat corresponds to row i of lu.  */
 
358
   /* lxm[i] records which row of lu  corresponds to row i of mat. */
 
359
   for (i = 0; i < n; i++) {
 
360
      lxm[mxl[i]] = i;
 
361
   }
 
362
 
 
363
 
 
364
   /* Determine the inverse matrix. */
 
365
   for (i = 0, ij = 0; i < n; i++) {
 
366
      for (j = 0; j < n; j++, ij++) {
 
367
         inv[ij] = 0.0;
 
368
      }
 
369
   }
 
370
 
 
371
   for (k = 0; k < n; k++) {
 
372
      inv[lxm[k]*n+k] = 1.0;
 
373
 
 
374
      /* Forward substitution. */
 
375
      for (i = lxm[k]+1; i < n; i++) {
 
376
         for (j = lxm[k]; j < i; j++) {
 
377
            inv[i*n+k] -= lu[i*n+j]*inv[j*n+k];
 
378
         }
 
379
      }
 
380
 
 
381
      /* Backward substitution. */
 
382
      for (i = n-1; i >= 0; i--) {
 
383
         for (j = i+1; j < n; j++) {
 
384
            inv[i*n+k] -= lu[i*n+j]*inv[j*n+k];
 
385
         }
 
386
         inv[i*n+k] /= lu[i*n+i];
 
387
      }
 
388
   }
 
389
   status = 0;                          /* all o.k. */
 
390
 
 
391
end_of_it:
 
392
   (void)free(mxl);
 
393
   (void)free(lxm);
 
394
   (void)free(rowmax);
 
395
   (void)free(lu);
 
396
 
 
397
   return status;
 
398
}