← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/ccdred/libsrc/mo_links.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/ccdred/libsrc/mo_links.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
.IDENTIFIER  mo_links.c
 
30
.AUTHOR      R.H. Warmels IPG-ESO Garching
 
31
.KEYWORDS    alignment software
 
32
.LANGUAGE    C
 
33
.PURPOSE     Routine to compute the shift for each subframe
 
34
.ENVIRONment MIDAS
 
35
             #include <ccd_def.h>      Symbols used by the ccd package
 
36
.VERSION     1.0     16-May-1995   creation
 
37
 
 
38
 090526         last modif
 
39
------------------------------------------------------------*/
 
40
/*
 
41
 * Define _POSIX_SOURCE to indicate
 
42
 * that this is a POSIX program
 
43
 */
 
44
#define  _POSIX_SOURCE 1
 
45
 
 
46
/*
 
47
 * definition of the used functions in this module
 
48
 */
 
49
 
 
50
#include <midas_def.h>
 
51
#include <stdio.h>
 
52
#include <string.h>
 
53
#include <stdlib.h>
 
54
#include <math.h>
 
55
 
 
56
#include <ccd_def.h>
 
57
 
 
58
int  MO_GET_SHIFTS();
 
59
void MO_MEDR();
 
60
 
 
61
/*
 
62
 MO_LINK -- Routine to compute the shifts directly. Input is given
 
63
            by an input table containing the coordinates of common objects
 
64
 */
 
65
 
 
66
int MO_LINKS(tblco, col, xrshift, yrshift, xcshift, ycshift, nrshift,
 
67
        ncshift, ncols, nrows, nxrsub, nyrsub, nxsub, nysub, nxoverlap,
 
68
        nyoverlap, order, nshifts)
 
69
int     tblco;
 
70
int     *col;                   /* columns containing the shifts */
 
71
float   (*xrshift)[MAXFRM];
 
72
float   (*yrshift)[MAXFRM];
 
73
float   (*xcshift)[MAXFRM];
 
74
float   (*ycshift)[MAXFRM];
 
75
int     (*nrshift)[MAXFRM];     /* number of row shifts */
 
76
int     (*ncshift)[MAXFRM];     /* number of column shifts */
 
77
int     ncols;                  /* number of columns per subraster */
 
78
int     nrows;                  /* number of rows per subraster */
 
79
int     nxrsub;                 /* column index of reference subraster */
 
80
int     nyrsub;                 /* column index of reference subraster */
 
81
int     nxsub;                  /* number of subrasters in X */
 
82
int     nysub;                  /* number of subrasters in y */
 
83
int     nxoverlap;              /* number of columns of overlap */
 
84
int     nyoverlap;              /* number of rows of overlap */
 
85
char    *order;                 /* row or column order */
 
86
int     *nshifts;
 
87
 
 
88
{       
 
89
int     i, j, nxsize, nysize, ilimit, olimit;
 
90
int     ns;
 
91
int     *nrowavg, *ncolavg;
 
92
 
 
93
float   *xcolavg, *ycolavg, *xrowavg, *yrowavg;
 
94
float   isign, jsign, xrmed, yrmed, xcmed, ycmed;
 
95
 
 
96
char  *osmmget();
 
97
 
 
98
 
 
99
 
 
100
if (strcmp(order,MO_COLUMN) == 0)
 
101
   {
 
102
   ilimit = nysub;
 
103
   olimit = nxsub;
 
104
   }
 
105
else 
 
106
   {
 
107
   ilimit = nxsub;
 
108
   olimit = nysub;
 
109
   }
 
110
        
 
111
/*
 
112
 Accumulate the shifts.
 
113
 */
 
114
nxsize  = ncols - nxoverlap;
 
115
nysize  = nrows - nyoverlap;
 
116
MO_GET_SHIFTS(tblco, col, xrshift, yrshift, nrshift, xcshift,
 
117
                        ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, 
 
118
                        nxoverlap, nyoverlap, nxsize, nysize, &ns);
 
119
if (ns == 0)
 
120
   return (0);
 
121
 
 
122
xcolavg = (float *) osmmget( olimit * sizeof( float ));   
 
123
ycolavg = (float *) osmmget( olimit * sizeof( float ));   
 
124
xrowavg = (float *) osmmget( olimit * sizeof( float ));   
 
125
yrowavg = (float *) osmmget( olimit * sizeof( float ));   
 
126
nrowavg = (int *) osmmget( olimit * sizeof( int ));   
 
127
ncolavg = (int *) osmmget( olimit * sizeof( int ));   
 
128
 
 
129
for (i = 0; i < olimit; i++) *(xcolavg+i) = 0.0;
 
130
for (i = 0; i < olimit; i++) *(ycolavg+i) = 0.0;
 
131
for (i = 0; i < olimit; i++) *(xrowavg+i) = 0.0;
 
132
for (i = 0; i < olimit; i++) *(yrowavg+i) = 0.0;
 
133
for (i = 0; i < olimit; i++) *(nrowavg+i) = 0.0;
 
134
for (i = 0; i < olimit; i++) *(ncolavg+i) = 0.0;
 
135
 
 
136
/*
 
137
 Compute the row or column sums.
 
138
 */
 
139
if (strcmp(order,MO_COLUMN) == 0)
 
140
   {
 
141
   for (i = 0; i < nxsub; i++) 
 
142
      {
 
143
      for (j = 0; j < nysub; j++)
 
144
          {
 
145
          if (nrshift[i][j] > 0) 
 
146
             {
 
147
             *(xrowavg+i) = *(xrowavg+i) + fabs(xrshift[i][j]);
 
148
             *(yrowavg+i) = *(yrowavg+i) + fabs(yrshift[i][j]);
 
149
             *(nrowavg+i) = *(nrowavg+i) + 1;
 
150
             } 
 
151
                    
 
152
          if (ncshift[i][j] > 0) 
 
153
             {
 
154
             *(xcolavg+i) = *(xcolavg+i) + fabs(xcshift[i][j]);
 
155
             *(ycolavg+i) = *(ycolavg+i) + fabs(ycshift[i][j]);
 
156
             *(ncolavg+i) = *(ncolavg+i) + 1;
 
157
             }
 
158
          }
 
159
      }
 
160
   } 
 
161
 
 
162
else 
 
163
   {
 
164
   for (i = 0; i < nysub; i++) 
 
165
      {
 
166
      for (j = 0; j < nxsub; j++)
 
167
          {
 
168
          if (nrshift[j][i] > 0) 
 
169
             {
 
170
             *(xrowavg+i) = *(xrowavg+i) + fabs(xrshift[j][i]);
 
171
             *(yrowavg+i) = *(yrowavg+i) + fabs(yrshift[j][i]);
 
172
             *(nrowavg+i) = *(nrowavg+i) + 1;
 
173
             } 
 
174
                   
 
175
          if (ncshift[j][i] > 0) 
 
176
             {
 
177
             *(xcolavg+i) = *(xcolavg+i) + fabs(xcshift[j][i]);
 
178
             *(ycolavg+i) = *(ycolavg+i) + fabs(ycshift[j][i]);
 
179
             *(ncolavg+i) = *(ncolavg+i) + 1;
 
180
             }
 
181
          }
 
182
      }
 
183
   }
 
184
        
 
185
/*
 
186
 Compute the averages.
 
187
 */
 
188
for (i = 0; i < olimit; i++)
 
189
   {
 
190
   if (*(nrowavg+i) > 0) 
 
191
      {
 
192
      *(xrowavg+i) = *(xrowavg+i) / *(nrowavg+i);
 
193
      *(yrowavg+i) = *(yrowavg+i) / *(nrowavg+i);
 
194
      }
 
195
   if ( (*ncolavg+i) > 0) 
 
196
      { 
 
197
      *(xcolavg+i) = *(xcolavg+i) / *(ncolavg+i);
 
198
      *(ycolavg+i) = *(ycolavg+i) / *(ncolavg+i);
 
199
      }
 
200
   }
 
201
 
 
202
/*
 
203
 Compute the medians of the row and column averages.
 
204
 */
 
205
MO_MEDR(xrowavg, nrowavg, olimit, &xrmed);
 
206
MO_MEDR(yrowavg, nrowavg, olimit, &yrmed);
 
207
MO_MEDR(xcolavg, ncolavg, olimit, &xcmed);
 
208
MO_MEDR(ycolavg, ncolavg, olimit, &ycmed);
 
209
 
 
210
/*
 
211
 Use the average shifts for subrasters with no information.
 
212
 */
 
213
for (j = 1; j <= nysub; j++)
 
214
   {
 
215
   if (j == nyrsub)
 
216
      jsign = 0.0;
 
217
   else if (j < nyrsub)
 
218
      jsign = 1.0;
 
219
   else
 
220
      jsign = -1.0;
 
221
 
 
222
   for (i = 1; i <= nxsub; i++)
 
223
      {
 
224
      if (i == nxrsub)
 
225
         isign = 0.0;
 
226
      else if (i < nxrsub)
 
227
         isign = 1.0;
 
228
      else
 
229
         isign = -1.0;
 
230
 
 
231
      if (nrshift[i-1][j-1] <= 0) 
 
232
         {
 
233
         if ( *(nrowavg+i-1) <= 0) 
 
234
            {
 
235
            xrshift[i-1][j-1] = isign * xrmed;
 
236
            yrshift[i-1][j-1] = jsign * yrmed;
 
237
            } 
 
238
         else if (strcmp(order,MO_COLUMN) == 0)
 
239
            {
 
240
            xrshift[i-1][j-1] = isign * *(xrowavg+i-1);
 
241
            yrshift[i-1][j-1] = jsign * *(yrowavg+i-1);
 
242
            } 
 
243
         else 
 
244
            {
 
245
            xrshift[i-1][j-1] = isign * *(xrowavg+j-1);
 
246
            yrshift[i-1][j-1] = jsign * *(yrowavg+j-1);
 
247
            }
 
248
         }
 
249
 
 
250
      if (ncshift[i-1][j-1] <= 0) 
 
251
         {
 
252
         if ( *(ncolavg+i-1) <= 0) 
 
253
            {
 
254
            xcshift[i-1][j-1] = isign * xcmed;
 
255
            ycshift[i-1][j-1] = jsign * ycmed;
 
256
            } 
 
257
         else if (strcmp(order,MO_COLUMN) == 0)
 
258
            {
 
259
            xcshift[i-1][j-1] = isign * *(xcolavg+i-1);
 
260
            ycshift[i-1][j-1] = jsign * *(ycolavg+i-1);
 
261
            } 
 
262
         else 
 
263
            {
 
264
            xcshift[i-1][j-1] = isign * *(xcolavg+j-1);
 
265
            ycshift[i-1][j-1] = jsign * *(ycolavg+j-1);
 
266
            }
 
267
         }
 
268
      }
 
269
   }
 
270
 
 
271
*nshifts = ns;
 
272
osmmfree((char *) xcolavg);
 
273
osmmfree((char *) xrowavg);
 
274
osmmfree((char *) ycolavg);
 
275
osmmfree((char *) yrowavg);
 
276
osmmfree((char *) nrowavg);
 
277
osmmfree((char *) ncolavg);
 
278
 
 
279
return 0;
 
280
}
 
281
 
 
282
 
 
283
/*
 
284
 MO_CLINKS -- Routine to compute the shift for each subframe. 
 
285
              Input is given by a fixed shift in x and y
 
286
*/
 
287
 
 
288
int MO_CLINKS(xrshift, yrshift, xcshift, ycshift, nxrsub, nyrsub,
 
289
          nxsub, nysub, xshift, yshift, nshifts)
 
290
 
 
291
float    (*xrshift)[MAXFRM];
 
292
float    (*yrshift)[MAXFRM];
 
293
float    (*xcshift)[MAXFRM];
 
294
float    (*ycshift)[MAXFRM];
 
295
int      nxrsub;
 
296
int      nyrsub;
 
297
int      nxsub;
 
298
int      nysub;
 
299
float    xshift;
 
300
float    yshift;
 
301
int      *nshifts;
 
302
 
 
303
{
 
304
int      i, j;
 
305
int      isign, jsign;
 
306
 
 
307
for (j = 0; j < nysub; j++)
 
308
    {
 
309
    if (j+1 == nyrsub)
 
310
       jsign = 0;
 
311
    else if (j+1 < nyrsub)
 
312
       jsign = 1;
 
313
    else
 
314
       jsign = -1;
 
315
 
 
316
    for (i = 0; i < nxsub; i++)
 
317
        {
 
318
        if (i+1 == nxrsub)
 
319
           isign = 0;
 
320
        else if (i+1 < nxrsub)
 
321
           isign = 1;
 
322
        else
 
323
           isign = -1;
 
324
 
 
325
        xrshift[i][j] = isign * fabs (xshift);
 
326
        yrshift[i][j] = 0.0;
 
327
        xcshift[i][j] = 0.0; 
 
328
        ycshift[i][j] = jsign * fabs (yshift);
 
329
        }
 
330
    }
 
331
*nshifts = 1;
 
332
return 0;
 
333
}
 
334
 
 
335
/*
 
336
 MO_FLINKS -- Routine to fetch the shifts directly from the table
 
337
 */
 
338
 
 
339
int MO_FLINKS(tblco,col,deltax, deltay, deltai, max_nshifts, nshifts)
 
340
int     tblco;                 /* table containing the shifts */
 
341
int     *col;                  /* columns containing the shifts */
 
342
float   *deltax;               /* x shifts */
 
343
float   *deltay;               /* y shifts */
 
344
float   *deltai;               /* intensity shifts */
 
345
int     max_nshifts;           /* maximum number of shifts */
 
346
int     *nshifts;
 
347
 
 
348
{
 
349
int     i;
 
350
int     ns;
 
351
int     ncolco;
 
352
int     nrowco;
 
353
int     nsortco;
 
354
int     allcolco;
 
355
int     allrowco;
 
356
float   data[3];
 
357
int     null;
 
358
 
 
359
(void) TCIGET(tblco, &ncolco, &nrowco, &nsortco, &allcolco, &allrowco);
 
360
ns = 0; 
 
361
 
 
362
for (i = 1; i <= nrowco; i++)
 
363
   {
 
364
   if (ns >= max_nshifts) 
 
365
      break;
 
366
   if (col[0] > 0 && col[1] > 0)
 
367
      {
 
368
      TCERDR(tblco,i,col[0],&data[0],&null);
 
369
      TCERDR(tblco,i,col[1],&data[1],&null);
 
370
      deltax[ns] = data[0];
 
371
      deltay[ns] = data[1];
 
372
      }
 
373
   else
 
374
      continue;
 
375
 
 
376
   if (col[2] == -1)
 
377
      {
 
378
      deltai[ns] = 0.0;
 
379
      }
 
380
   else
 
381
      TCERDR(tblco,i,col[2],&data[2],&null);
 
382
      deltai[ns] = data[2];
 
383
   ns++;
 
384
   }
 
385
*nshifts = ns;
 
386
return 0;
 
387
}
 
388
 
 
389
 
 
390
/*
 
391
 MO_GET_SHIFTS -- Procedure to accumulate shifts for each subraster.
 
392
 */
 
393
 
 
394
int MO_GET_SHIFTS(tblco, col, xrshift, yrshift, nrshift, xcshift,
 
395
        ycshift, ncshift, nxsub, nysub, nxrsub, nyrsub, nxoverlap,
 
396
        nyoverlap, nxsize, nysize, nshift)
 
397
 
 
398
int     tblco;
 
399
int     col[3];
 
400
float   (*xrshift)[MAXFRM];
 
401
float   (*yrshift)[MAXFRM];
 
402
int     (*nrshift)[MAXFRM];
 
403
float   (*xcshift)[MAXFRM];
 
404
float   (*ycshift)[MAXFRM];
 
405
int     (*ncshift)[MAXFRM];     /* number of column shifts */
 
406
int     nxsub;                  /* number of subrasters in x */
 
407
int     nysub;                  /* number of subrasters in y */
 
408
int     nxrsub;                 /* column index of reference subraster */
 
409
int     nyrsub;                 /* column index of reference subraster */
 
410
int     nxoverlap;              /* number of columns of overlap */
 
411
int     nyoverlap;              /* number of rows of overlap */
 
412
int     nxsize;                 /* number of columns per subraster */
 
413
int     nysize;                 /* number of rows per subraster */
 
414
int     *nshift;
 
415
 
 
416
{
 
417
int     allcolco, allrowco;
 
418
int     i, irow;
 
419
int     j;
 
420
int     ncolco, nrowco, nsortco;
 
421
int     ns;
 
422
int     nx1, ny1, nx11, ny11;
 
423
int     nx2, ny2, nx21, ny21;
 
424
int     r21, r22;
 
425
int     stat1, stat2;
 
426
float   x1, y1, x2, y2;
 
427
float   xdif, xdifm, ydif, ydifm;
 
428
int     null;
 
429
 
 
430
ns = 0;
 
431
xdifm = ydifm = 0.0;
 
432
nx2 = ny2 = 0;
 
433
 
 
434
(void) TCIGET(tblco, &ncolco, &nrowco, &nsortco, &allcolco, &allrowco);
 
435
irow = 1;
 
436
 
 
437
do {
 
438
   if (col[0] == -1 || col[1] == -1)
 
439
      continue;
 
440
   stat1 = TCERDR(tblco,irow,col[0],&x1,&null);
 
441
   stat2 = TCERDR(tblco,irow,col[1],&y1,&null);
 
442
   if (stat1 !=0 || stat2 != 0) 
 
443
      continue;
 
444
   
 
445
/*
 
446
 Compute which subraster 1 belongs to.
 
447
 */
 
448
   if (fmod( (int) x1, nxsize) == 0)
 
449
      nx1 = (int) x1 / nxsize;
 
450
   else
 
451
      nx1 = (int) x1 / nxsize + 1;
 
452
 
 
453
   if (fmod( (int) y1, nysize) == 0)
 
454
      ny1 = (int) y1 / nysize;
 
455
   else
 
456
      ny1 = (int) y1 / nysize + 1;
 
457
 
 
458
/*
 
459
 Get the second coordinate pair.
 
460
 */
 
461
   irow++;
 
462
   do {
 
463
      if (irow > nrowco) break;
 
464
      stat1 = TCERDR(tblco,irow,col[0],&x2,&null);
 
465
      stat2 = TCERDR(tblco,irow,col[1],&y2,&null);
 
466
 
 
467
/*
 
468
 Compute which subraster 2 belongs to.
 
469
 */
 
470
      if ((stat1 == 0) & (stat2 == 0))
 
471
         {
 
472
         if (fmod( (int) x2, nxsize) == 0)
 
473
            nx2 = (int) x2 / nxsize;
 
474
         else
 
475
            nx2 = (int) x2 / nxsize + 1;
 
476
 
 
477
         if (fmod( (int) y2, nysize) == 0)
 
478
            ny2 = (int) y2 / nysize;
 
479
         else
 
480
            ny2 = (int) y2 / nysize + 1;
 
481
         }
 
482
      }
 
483
   while (stat1 != 0 || stat2 != 0);
 
484
 
 
485
 
 
486
   if ((irow > nrowco) || (stat1 != 0 || stat2 != 0))
 
487
      break;
 
488
 
 
489
   r21 = (nx1 - nxrsub) * (nx1 - nxrsub) + (ny1 - nyrsub) * (ny1 - nyrsub);
 
490
   r22 = (nx2 - nxrsub) * (nx2 - nxrsub) + (ny2 - nyrsub) * (ny2 - nyrsub);
 
491
 
 
492
/*
 
493
 Compute the shift for the first subraster.
 
494
 */
 
495
   if (r21 == r22) 
 
496
      continue;
 
497
 
 
498
   else if (r21 > r22) 
 
499
      {
 
500
      xdif = x2 - x1;
 
501
      if (nxoverlap < 0) 
 
502
         {
 
503
         if (xdif < 0.0)
 
504
            xdifm = xdif - nxoverlap;
 
505
         else if (xdif > 0.0)
 
506
            xdifm = xdif + nxoverlap;
 
507
         } 
 
508
      else
 
509
         xdifm = xdif;
 
510
 
 
511
      ydif = y2 - y1;
 
512
      if (nyoverlap < 0) 
 
513
         {
 
514
         if (ydif < 0.0)
 
515
            ydifm = ydif - nyoverlap;
 
516
         else if (ydif > 0.0)
 
517
            ydifm = ydif + nyoverlap;
 
518
         } 
 
519
      else
 
520
         ydifm = ydif;
 
521
 
 
522
      nx11 = nx1 - 1;
 
523
      ny11 = ny1 - 1;
 
524
      if (nx1 == nx2) 
 
525
         {
 
526
         xcshift[nx11][ny11] = xcshift[nx11][ny11] + xdif;
 
527
         ycshift[nx11][ny11] = ycshift[nx11][ny11] + ydifm;
 
528
         ncshift[nx11][ny11] = ncshift[nx11][ny11] + 1;
 
529
         } 
 
530
      else if (ny1 == ny2) 
 
531
         {
 
532
         xrshift[nx11][ny11] = xrshift[nx1][ny11] + xdifm;
 
533
         yrshift[nx11][ny11] = yrshift[nx1][ny11] + ydif;
 
534
         nrshift[nx11][ny11] = nrshift[nx1][ny11] + 1;
 
535
         } 
 
536
      else
 
537
         continue;
 
538
 
 
539
      ns++;
 
540
      }
 
541
 
 
542
/*
 
543
 Compute the shift for the second subraster.
 
544
 */
 
545
   else 
 
546
      {
 
547
      xdif = x1 - x2;
 
548
      if (nxoverlap < 0) 
 
549
         {
 
550
         if (xdif < 0.0)
 
551
            xdifm = xdif - nxoverlap;
 
552
         else if (xdif > 0.0)
 
553
            xdifm = xdif + nxoverlap;
 
554
         } 
 
555
      else
 
556
         xdifm = xdif;
 
557
 
 
558
      ydif = y1 - y2;
 
559
      if (nyoverlap < 0) 
 
560
         {
 
561
         if (ydif < 0.0)
 
562
            ydifm = ydif - nyoverlap;
 
563
         else if (ydif > 0.0)
 
564
            ydifm = ydif + nyoverlap;
 
565
         } 
 
566
      else
 
567
         ydifm = ydif;
 
568
 
 
569
 
 
570
      nx21 = nx2 - 1;
 
571
      ny21 = ny2 - 1;
 
572
      if (nx1 == nx2) 
 
573
         {
 
574
         xcshift[nx21][ny21] = xcshift[nx21][ny21] + xdif;
 
575
         ycshift[nx21][ny21] = ycshift[nx21][ny21] + ydifm;
 
576
         ncshift[nx21][ny21] = ncshift[nx21][ny21] + 1;
 
577
         } 
 
578
    
 
579
      else if (ny1 == ny2) 
 
580
         {
 
581
         xrshift[nx21][ny21] = xrshift[nx21][ny21] + xdifm;
 
582
         yrshift[nx21][ny21] = yrshift[nx21][ny21] + ydif;
 
583
         nrshift[nx21][ny21] = nrshift[nx21][ny21] + 1;
 
584
         } 
 
585
      else
 
586
         continue;
 
587
      }
 
588
      irow++;
 
589
      ns++;
 
590
   }
 
591
while (irow <=  nrowco);
 
592
 
 
593
/*
 
594
 Compute the final shifts.
 
595
 */
 
596
for (j = 0; j < nysub; j++)
 
597
   {
 
598
   for (i = 0; i < nxsub; i++)
 
599
      {
 
600
      if (nrshift[i][j] > 0) 
 
601
         {
 
602
         xrshift[i][j] = xrshift[i][j] / nrshift[i][j];
 
603
         yrshift[i][j] = yrshift[i][j] / nrshift[i][j];
 
604
         }
 
605
 
 
606
      if (ncshift[i][j] > 0) 
 
607
         {
 
608
         xcshift[i][j] = xcshift[i][j] / ncshift[i][j];
 
609
         ycshift[i][j] = ycshift[i][j] / ncshift[i][j];
 
610
         }
 
611
      }
 
612
   }
 
613
*nshift = ns;
 
614
 
 
615
return 0;
 
616
}
 
617
 
 
618
 
 
619
 
 
620
 
 
621