← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/display/src/smosub3.for

  • 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/display/src/smosub3.for
 
1
C @(#)smosub3.for       13.1.1.1 (ESO-DMD) 06/02/98 18:11:31
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
      SUBROUTINE NORMIT(A,NDIM)
 
30
C
 
31
      IMPLICIT NONE
 
32
C
 
33
      INTEGER  NDIM,N
 
34
      REAL     A(NDIM),R
 
35
      DOUBLE PRECISION   SUM,DR
 
36
C
 
37
      SUM = 0.D0
 
38
      DO 1000, N=1,NDIM
 
39
         SUM = SUM + A(N)
 
40
1000  CONTINUE
 
41
C 
 
42
      IF (SUM.GT.10.D-20) THEN
 
43
         DR = 1.D0 / SUM
 
44
         R = SNGL(DR)
 
45
         DO 2000, N=1,NDIM
 
46
            A(N) = A(N) * R
 
47
2000     CONTINUE
 
48
      ENDIF
 
49
C 
 
50
      RETURN
 
51
      END
 
52
 
 
53
      SUBROUTINE TEMPLB(FLAG,A,B,NPIX,WEITS,RADIUS,CUTS,STAT)
 
54
C
 
55
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
56
C       
 
57
C.IDENTIFICATION
 
58
C  subroutine TEMPLB                    version 1.50    860610
 
59
C  K. Banse                             ESO - Garching
 
60
C
 
61
C.KEYWORDS
 
62
C  template matching, neighbourhood
 
63
C
 
64
C.PURPOSE
 
65
C replace pixels by weighted sum of nw*nw matrix (nw = 2*radius + 1)
 
66
C
 
67
C.ALGORITHM
 
68
C  filter image by taking the weighted sum W(x,y) of all points in neighbourhood
 
69
C  of x,y with given radius, and replace original value I(x,y) by W(x,y)
 
70
C       
 
71
C.INPUT/OUTPUT
 
72
C  call as TEMPLATE(A,B,NPIX,WEITS,RADIUS,STAT)
 
73
C
 
74
C input par:
 
75
C  FLAG:        I*4             function flag:
 
76
C                               = 1 for Gauss filter
 
77
C                               = 2 for digital filter
 
78
C                               = 3 for max filter
 
79
C                               = 4 for min filter
 
80
C  A:           R*4 array       input array
 
81
C  B:           R*4 array       output array
 
82
C  NPIX:        I*4 array       no. of pixels per line, no. of lines
 
83
C  WEITS:       R*4 array       weights for NWX*NWY matrix, ordered as
 
84
C                               1 2 3 ... NWX
 
85
C                               NWX+1 ... 2NWX
 
86
C                               ............
 
87
C                               (NWY-1)*NWX+1 ... NWX*NWY
 
88
C  RADIUS:      I*4 array       radius of template matrix in x,y
 
89
C       
 
90
C  output par:
 
91
C  CUTS:        R*4 array       cut values of result frame
 
92
C  STAT:        I*4             return status, should be = 0, else problems...
 
93
C       
 
94
C --------------------------------------------------------------------
 
95
C
 
96
      IMPLICIT NONE
 
97
C
 
98
      INTEGER     FLAG,NPIX(*),RADIUS(*)
 
99
      INTEGER     NX,NY,DIAMX,DIAMY,RADYNX
 
100
      INTEGER     N,NP,NL,BINDX,STAT
 
101
      INTEGER     WINDX,OFFSET,OFF
 
102
      INTEGER     NRDXP1,RADXP1,NXRADX,LEFT
 
103
C
 
104
      REAL        A(*),B(*),WEITS(*),CUTS(*)
 
105
      REAL        TEST
 
106
C 
 
107
      DOUBLE PRECISION  DTEST
 
108
C
 
109
C  init variables + test size of template
 
110
      BINDX = 1                                !index for frame B
 
111
      NX = NPIX(1)
 
112
      NY = NPIX(2)
 
113
      DIAMX = RADIUS(1) * 2
 
114
      DIAMY = RADIUS(2) * 2
 
115
      IF ((DIAMX.GE.NX).OR.(DIAMY.GE.NY)) THEN        !data matrix at least NW
 
116
         STAT = 1
 
117
         RETURN
 
118
      ELSE
 
119
         STAT = 0
 
120
      ENDIF
 
121
C
 
122
      RADYNX = RADIUS(2) * NX                  !all following constants for
 
123
      RADXP1 = RADIUS(1) + 1                   !speed reasons...
 
124
      NXRADX = NX - RADIUS(1)
 
125
      NRDXP1 = NXRADX + 1
 
126
      LEFT = RADYNX + RADIUS(1)            !lower left corner of template matrix
 
127
      OFFSET = RADYNX
 
128
C
 
129
C  branch according to filter type
 
130
      IF (FLAG.GE.3) GOTO 5000
 
131
C
 
132
C*******
 
133
C
 
134
C  main loop over inner lines
 
135
C
 
136
C*******
 
137
C
 
138
      DO 2000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
 
139
C
 
140
C  step through each pixel of data frame
 
141
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
 
142
C 
 
143
         DO 1500, NP=OFFSET+RADXP1,OFFSET+NXRADX 
 
144
            DTEST = 0.D0                              !clear sum
 
145
            WINDX = 1                                !index for weight matrix
 
146
C
 
147
C  sum + weigh inner part
 
148
            DO 1000, OFF=NP-LEFT,NP+LEFT,NX        !step through lines of window
 
149
               DO 800, N=OFF,OFF+DIAMX
 
150
                  DTEST = DTEST + ( A(N) * WEITS(WINDX) )
 
151
                  WINDX = WINDX + 1
 
152
800            CONTINUE
 
153
1000        CONTINUE
 
154
C
 
155
            B(BINDX) = DTEST                     !fill output frame
 
156
            BINDX = BINDX + 1
 
157
1500     CONTINUE
 
158
C
 
159
C  update main offset
 
160
         OFFSET = OFFSET + NX
 
161
C
 
162
2000  CONTINUE
 
163
C
 
164
C finally update the cuts
 
165
      CALL CHCUTS(B,BINDX,CUTS)
 
166
      RETURN
 
167
C
 
168
C  here for min/max filter
 
169
C 
 
170
5000  IF (FLAG.EQ.3) THEN
 
171
C 
 
172
C  max filter
 
173
         DO 7000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
 
174
C
 
175
C  step through each pixel of data frame
 
176
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
 
177
C
 
178
            DO 6000, NP=OFFSET+RADXP1,OFFSET+NXRADX
 
179
               TEST = A(NP-LEFT)                        !take lower left pixel
 
180
C
 
181
C  find max in inner part
 
182
               DO 5500, OFF=NP-LEFT,NP+LEFT,NX    
 
183
                  DO 5400, N=OFF,OFF+DIAMX
 
184
                     IF (A(N).GT.TEST) TEST = A(N)
 
185
5400              CONTINUE
 
186
5500           CONTINUE
 
187
C
 
188
               B(BINDX) = TEST                     !fill output frame
 
189
               BINDX = BINDX + 1
 
190
6000        CONTINUE
 
191
C
 
192
            OFFSET = OFFSET + NX                   !update main offset
 
193
7000     CONTINUE
 
194
C 
 
195
      ELSE
 
196
C 
 
197
C  min filter
 
198
         DO 9000, NL=RADIUS(2)+1,NY-RADIUS(2)            !step through each line
 
199
C
 
200
C  step through each pixel of data frame
 
201
C  i.e., OFFSET+RADIUS(1)+1  >>>  OFF SET+NX-RADIUS(1)
 
202
C
 
203
            DO 8800, NP=OFFSET+RADXP1,OFFSET+NXRADX
 
204
               TEST = A(NP-LEFT)                        !take lower left pixel
 
205
C
 
206
C  find min in inner part
 
207
               DO 8600, OFF=NP-LEFT,NP+LEFT,NX   
 
208
                  DO 8400, N=OFF,OFF+DIAMX
 
209
                     IF (A(N).LT.TEST) TEST = A(N)
 
210
8400              CONTINUE
 
211
8600           CONTINUE
 
212
C
 
213
               B(BINDX) = TEST                     !fill output frame
 
214
               BINDX = BINDX + 1
 
215
8800        CONTINUE
 
216
C
 
217
            OFFSET = OFFSET + NX                   !update main offset
 
218
9000     CONTINUE
 
219
      ENDIF
 
220
C
 
221
C finally update the cuts
 
222
      CALL CHCUTS(B,BINDX,CUTS)
 
223
      RETURN
 
224
C
 
225
      END
 
226
 
 
227
      SUBROUTINE CHCUTS(A,ASZ1,CUTS)
 
228
C 
 
229
      INTEGER     ASZ1
 
230
      INTEGER     N,NDIM
 
231
C 
 
232
      REAL        A(*),CUTS(*)
 
233
      REAL        TEST
 
234
C 
 
235
      NDIM = ASZ1 - 1
 
236
      IF (NDIM.LT.3) RETURN
 
237
C 
 
238
      IF (A(1).GT.A(2)) THEN
 
239
         CUTS(1) = A(2)
 
240
         CUTS(2) = A(1)
 
241
      ELSE
 
242
         CUTS(1) = A(1)
 
243
         CUTS(2) = A(2)
 
244
      ENDIF
 
245
C 
 
246
      DO 1000, N=3,NDIM
 
247
         TEST = A(N)
 
248
         IF (TEST.LT.CUTS(1)) THEN
 
249
            CUTS(1) = TEST
 
250
         ELSE IF (TEST.GT.CUTS(2)) THEN
 
251
            CUTS(2) = TEST
 
252
         ENDIF
 
253
1000  CONTINUE
 
254
C 
 
255
      RETURN
 
256
      END
 
257
 
 
258
      SUBROUTINE EXPAX(IMNOA,OFFPIX,OUTFR,INBUF,RSTAT)
 
259
C
 
260
C++++++++++++++++++++++++++++++++++++++++++++++++++
 
261
C
 
262
C.IDENTIFICATION
 
263
C  subroutine EXPX
 
264
C  K. Banse                             ESO - Garching
 
265
C  1.00    910517
 
266
C
 
267
C.KEYWORDS
 
268
C  border regions
 
269
C
 
270
C.PURPOSE
 
271
C  expand image by duplicating border lines + columns
 
272
C
 
273
C.ALGORITHM
 
274
C  get the no. of columns to expand at left + at right
 
275
C  get the no. of lines to expand in the beginning + end
 
276
C  flip columns in x, lines in y and add to image
 
277
C
 
278
C.INPUT/OUTPUT
 
279
C  call via
 
280
C   EXPAX(IMNOA,OFFPIX,OUTFR,INBUF,RSTAT)
 
281
C
 
282
C  input par:
 
283
C  IMNOA        integer         file no. of input frame
 
284
C  OFFPIX       integer         offset for 1. pixel to map
 
285
C  OUTFR        char. array     output frame
 
286
C  INBUF        integer array   xl,xr,yb,ye: no. of cols (left/right)
 
287
C                                            no. of lines (beginning/end)
 
288
C  output par:
 
289
C  RSTAT        integer         return status, 0 = o.k.
 
290
C
 
291
C.VERSIONS
 
292
C  1.00   creation
 
293
C       
 
294
C--------------------------------------------------
 
295
C
 
296
      IMPLICIT NONE
 
297
C
 
298
      INTEGER      INBUF(4),RSTAT,OFFPIX
 
299
      INTEGER      IAV,N,NN,LOOPLM,RMAIND,CHUNKI,CHUNKO
 
300
      INTEGER      NAXISA,MOFF,LOFF,NOLIN
 
301
      INTEGER*8    PNTRX,PNTRY
 
302
      INTEGER      IMNOA,IMNOC,IMNOX,IMNOY
 
303
      INTEGER      TRUSIZ,SIZE,STAT
 
304
      INTEGER      NPIXA(2),NPIXC(2)
 
305
      INTEGER      UNI(1),NULO,MADRID(1)
 
306
C
 
307
      CHARACTER*(*) OUTFR
 
308
      CHARACTER    CUNITA*64,IDENTA*72
 
309
C 
 
310
      REAL         CUTVAL(4)
 
311
C
 
312
      DOUBLE PRECISION    STEPA(2),STARTA(2)
 
313
      DOUBLE PRECISION    STARTC(2)
 
314
C
 
315
      COMMON  /VMR/  MADRID
 
316
C
 
317
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
318
C
 
319
      DATA    NPIXA  /1,1/,   NPIXC /1,1/
 
320
      DATA    STARTA /0.0,0.0/,  STEPA  /1.0,1.0/
 
321
      DATA    CUNITA /' '/,   IDENTA /' '/
 
322
C 
 
323
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
324
C
 
325
C 
 
326
      RSTAT = 0
 
327
      CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT)
 
328
      IF (NAXISA.GT.2) NAXISA = 2                                !max 2-dim 
 
329
      CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT)
 
330
      CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT)
 
331
      CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT)
 
332
      CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENTA,UNI,NULO,STAT)
 
333
      N = (NAXISA+1) * 16
 
334
      CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNITA,UNI,NULO,STAT)
 
335
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTVAL,UNI,NULO,STAT)
 
336
C
 
337
C  get size of input + output frame
 
338
      SIZE = NPIXA(1) * NPIXA(2)
 
339
      IF ( (INBUF(1).GT.NPIXA(1)) .OR.
 
340
     +     (INBUF(2).GT.NPIXA(1)) .OR.
 
341
     +     (INBUF(1).LT.0) .OR.
 
342
     +     (INBUF(2).LT.0) ) THEN
 
343
         RSTAT = 1
 
344
         RETURN
 
345
      ENDIF
 
346
C 
 
347
      NPIXC(1) = NPIXA(1) + INBUF(1) + INBUF(2)
 
348
      IF (NAXISA.GT.1) THEN
 
349
         IF ( (INBUF(3).GT.NPIXA(2)) .OR.
 
350
     +        (INBUF(4).GT.NPIXA(2)) .OR.
 
351
     +        (INBUF(3).LT.0) .OR.
 
352
     +        (INBUF(4).LT.0) ) THEN
 
353
            RSTAT = 2
 
354
            RETURN
 
355
         ENDIF
 
356
         NPIXC(2) = NPIXA(2) + INBUF(3) + INBUF(4)
 
357
      ENDIF
 
358
C
 
359
C create result frame
 
360
      IF (OFFPIX.LE.0) THEN
 
361
         TRUSIZ = NPIXC(1) * NPIXC(2)
 
362
         CALL STFCRE(OUTFR,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
363
     +               TRUSIZ,IMNOC,STAT)
 
364
C 
 
365
         CALL STDWRI(IMNOC,'NAXIS',NAXISA,1,1,UNI,STAT)
 
366
         CALL STDWRI(IMNOC,'NPIX',NPIXC,1,NAXISA,UNI,STAT)
 
367
         STARTC(1) = STARTA(1) - (INBUF(1)*STEPA(1))
 
368
         IF (NAXISA.GT.1)
 
369
     +      STARTC(2) = STARTA(2) - (INBUF(3)*STEPA(2))
 
370
         CALL STDWRD(IMNOC,'START',STARTC,1,NAXISA,UNI,STAT)
 
371
         CALL STDWRD(IMNOC,'STEP',STEPA,1,NAXISA,UNI,STAT)
 
372
         CALL STDWRC(IMNOC,'IDENT',1,IDENTA,1,72,UNI,STAT)
 
373
         N = (NAXISA+1) * 16
 
374
         CALL STDWRC(IMNOC,'CUNIT',1,CUNITA,1,N,UNI,STAT)
 
375
         CALL STDWRR(IMNOC,'LHCUTS',CUTVAL,1,4,UNI,STAT)
 
376
      ELSE
 
377
         CALL STFOPN(OUTFR,D_R4_FORMAT,0,F_IMA_TYPE,IMNOC,STAT)
 
378
      ENDIF
 
379
C 
 
380
C  get buffersize of at least one line but no more than total size
 
381
      NN = 40000                                !first guess ...
 
382
      IF (NN.LT.NPIXA(1)) NN = NPIXA(1)
 
383
      IF (NN.GT.SIZE) NN = SIZE
 
384
      NOLIN = NN / NPIXA(1)                     !buffer = NOLIN * NPIXA(1)
 
385
      CHUNKI = NOLIN * NPIXA(1)
 
386
      CHUNKO = NOLIN * NPIXC(1)
 
387
      IF (NAXISA.GT.1) THEN
 
388
         LOOPLM = NPIXA(2) / NOLIN              !get no. of passes through loop
 
389
         RMAIND = NPIXA(2) - (LOOPLM*NOLIN)     !and remainder
 
390
      ELSE
 
391
         LOOPLM = 1
 
392
         RMAIND = 0
 
393
      ENDIF
 
394
C 
 
395
      CALL STFCRE('expaxi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
396
     +            CHUNKI,IMNOX,STAT)
 
397
      CALL STFMAP(IMNOX,F_X_MODE,1,CHUNKI,IAV,PNTRX,STAT)
 
398
      CALL STFCRE('expaxo',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
399
     +            CHUNKO,IMNOY,STAT)
 
400
      CALL STFMAP(IMNOY,F_X_MODE,1,CHUNKO,IAV,PNTRY,STAT)
 
401
C 
 
402
C  now we loop through input frame
 
403
      MOFF = 1 + OFFPIX                        !point to 1. pixel in input fr.
 
404
      IF (NAXISA.GT.1) THEN                    !point to 1. pixel in output fr.
 
405
         LOFF = (INBUF(3)*NPIXC(1)) + 1       
 
406
      ELSE
 
407
         LOFF = 1
 
408
      ENDIF
 
409
C 
 
410
      DO 2000, N=1,LOOPLM
 
411
         CALL STFGET(IMNOA,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT)    !data in
 
412
         CALL EXPA1(MADRID(PNTRX),MADRID(PNTRY),NOLIN,NPIXA(1),
 
413
     +              NPIXC(1),INBUF)
 
414
         CALL STFPUT(IMNOC,LOFF,CHUNKO,MADRID(PNTRY),STAT)          !data out
 
415
         MOFF = MOFF + CHUNKI
 
416
         LOFF = LOFF + CHUNKO
 
417
2000  CONTINUE
 
418
C 
 
419
      IF (RMAIND.GT.0) THEN
 
420
         CHUNKI = RMAIND * NPIXA(1)
 
421
         CHUNKO = RMAIND * NPIXC(1)
 
422
         CALL STFGET(IMNOA,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT) 
 
423
         CALL EXPA1(MADRID(PNTRX),MADRID(PNTRY),RMAIND,NPIXA(1),
 
424
     +              NPIXC(1),INBUF)
 
425
         CALL STFPUT(IMNOC,LOFF,CHUNKO,MADRID(PNTRY),STAT)  
 
426
      ENDIF
 
427
C 
 
428
C  for 1-dim frames we are finished
 
429
      CALL STFCLO(IMNOX,STAT)
 
430
      CALL STFCLO(IMNOY,STAT)
 
431
      IF (NAXISA.EQ.1) GOTO 8000
 
432
C 
 
433
C  get buffersize of INBUF(3/4)*2 lines
 
434
      IF (INBUF(3).GE.INBUF(4)) THEN
 
435
         NN = INBUF(3)
 
436
      ELSE
 
437
         NN = INBUF(4)
 
438
      ENDIF
 
439
      IF (NN.EQ.0) GOTO 8000                !nothing to do in Y
 
440
C 
 
441
      NN = (2*NN + 1) * NPIXC(1)
 
442
      CALL STFCRE('expaxz',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
 
443
     +            NN,IMNOX,STAT)
 
444
      CALL STFMAP(IMNOX,F_X_MODE,1,NN,IAV,PNTRX,STAT)
 
445
C  
 
446
      NN = 2*INBUF(3) + 1
 
447
      CHUNKI = NN * NPIXC(1)
 
448
      IF (CHUNKI.GT.0) THEN
 
449
         CALL STFGET(IMNOC,1,CHUNKI,IAV,MADRID(PNTRX),STAT) 
 
450
         CALL EXPA2(1,MADRID(PNTRX),NPIXC,INBUF(3))
 
451
         CALL STFPUT(IMNOC,1,CHUNKI,MADRID(PNTRX),STAT)  
 
452
      ENDIF
 
453
C  
 
454
      NN = 2*INBUF(4) + 1
 
455
      CHUNKI = NN * NPIXC(1)
 
456
      IF (CHUNKI.GT.0) THEN
 
457
         MOFF = ((NPIXC(2)-NN)*NPIXC(1)) + 1
 
458
         CALL STFGET(IMNOC,MOFF,CHUNKI,IAV,MADRID(PNTRX),STAT) 
 
459
         CALL EXPA2(2,MADRID(PNTRX),NPIXC,INBUF(4))
 
460
         CALL STFPUT(IMNOC,MOFF,CHUNKI,MADRID(PNTRX),STAT)  
 
461
      ENDIF
 
462
      CALL STFCLO(IMNOX,STAT)
 
463
C 
 
464
8000  CALL STFCLO(IMNOC,STAT)
 
465
      RETURN
 
466
C 
 
467
      END
 
468
 
 
469
      SUBROUTINE EXPA1(X,Y,NOLIN,NPIXA,NPIXC,INBUF)
 
470
C 
 
471
      IMPLICIT NONE
 
472
C 
 
473
      REAL     X(*),Y(*)
 
474
C 
 
475
      INTEGER  NOLIN,NPIXA,NPIXC,INBUF(*)
 
476
      INTEGER  N,NN,IDX,KDX,IOFF,KOFF
 
477
C 
 
478
C  loop over all lines
 
479
      IOFF = 0
 
480
      KOFF = 0
 
481
      DO 5000, NN=1,NOLIN
 
482
C 
 
483
C fill first Y pixels with first flipped X pixels
 
484
         IF (INBUF(1).GT.0) THEN
 
485
            IDX = 1 + IOFF
 
486
            KDX = INBUF(1) + KOFF + 1
 
487
            DO 1000, N=1,INBUF(1)
 
488
               Y(IDX) = X(KDX)
 
489
               IDX = IDX + 1
 
490
               KDX = KDX - 1
 
491
1000        CONTINUE
 
492
         ELSE
 
493
            IDX = 1 + IOFF
 
494
            KDX = 1 + KOFF
 
495
         ENDIF
 
496
C 
 
497
C copy center Y pixels from all X pixels
 
498
         DO 2000, N=1,NPIXA
 
499
            Y(IDX) = X(KDX)
 
500
            IDX = IDX + 1
 
501
            KDX = KDX + 1
 
502
2000     CONTINUE
 
503
C 
 
504
C fill last Y pixels with last flipped X pixels
 
505
         KDX = KDX - 1
 
506
         IF (INBUF(2).GT.0) THEN
 
507
            DO 3000, N=1,INBUF(2)
 
508
               KDX = KDX - 1
 
509
               Y(IDX) = X(KDX)
 
510
               IDX = IDX + 1
 
511
3000        CONTINUE
 
512
         ENDIF
 
513
C 
 
514
         IOFF = IOFF + NPIXC
 
515
         KOFF = KOFF + NPIXA
 
516
5000  CONTINUE
 
517
C 
 
518
      RETURN
 
519
      END
 
520
 
 
521
      SUBROUTINE EXPA2(FLAG,X,NPIXC,INB)
 
522
C 
 
523
      IMPLICIT NONE
 
524
C 
 
525
      REAL     X(*)
 
526
C 
 
527
      INTEGER  FLAG,NPIXC(*),INB
 
528
      INTEGER  N,NN,IDX,KDX,IOFF,KOFF,MX,MMX
 
529
C 
 
530
      MX = NPIXC(1)
 
531
C 
 
532
      IF (FLAG.EQ.1) THEN
 
533
         IOFF = 0
 
534
         KOFF = 2*INB * MX
 
535
         MMX = MX
 
536
      ELSE
 
537
         IOFF = 2*INB * MX
 
538
         KOFF = 0
 
539
         MMX = -MX
 
540
      ENDIF
 
541
C 
 
542
      DO 2000, NN=1,INB
 
543
         IDX = IOFF + 1
 
544
         KDX = KOFF + 1
 
545
         DO 1000, N=1,MX
 
546
            X(IDX) = X(KDX)
 
547
            IDX = IDX + 1
 
548
            KDX = KDX + 1
 
549
1000     CONTINUE
 
550
C 
 
551
         IOFF = IOFF + MMX
 
552
         KOFF = KOFF - MMX
 
553
2000  CONTINUE
 
554
C 
 
555
      RETURN
 
556
      END
 
557
 
 
558
      SUBROUTINE SORTIT(RA,N,NH,LL)
 
559
C 
 
560
C  N  = full dimension of RA
 
561
C  NH = index of median - we only sort fully until NH
 
562
C  then we insert
 
563
C 
 
564
C  this is the Heapsort algorithm from "Numerical Recipes", page 231
 
565
C 
 
566
      IMPLICIT NONE
 
567
C 
 
568
      INTEGER   N,NH,LL
 
569
      INTEGER   L,IR,J,I
 
570
      INTEGER   JL,JU,JM,M
 
571
C 
 
572
      REAL      RA(N),RRA
 
573
C 
 
574
      L = LL                      !that is L = NH/2 + 1  
 
575
      IR = NH
 
576
C 
 
577
10    IF (L.GT.1) THEN            !still hiring
 
578
         L = L - 1
 
579
         RRA = RA(L)
 
580
      ELSE                        !in retirement + promotion phase
 
581
         RRA = RA(IR)             !clear a space at end of array
 
582
         RA(IR) = RA(1)           !retire the top of the heap into it
 
583
         IR = IR - 1           !decrease the size of the corporation
 
584
         IF (IR.EQ.1) THEN        !done with last promotion
 
585
            RA(1) = RRA           !the least competent guy of all...
 
586
            GOTO 1000             !now move to insertion part
 
587
         ELSE
 
588
         ENDIF
 
589
      ENDIF
 
590
C 
 
591
      I = L                       !in hiring as well as in promotion phase
 
592
      J = L + L              !we here set up to sift down RRA to its right level
 
593
C 
 
594
20    IF (J.LE.IR) THEN           !"Do while J<=IR"
 
595
         IF (J.LT.IR) THEN
 
596
            IF (RA(J).LT.RA(J+1)) J = J + 1   !compare to the better underling
 
597
         ENDIF
 
598
         IF (RRA.LT.RA(J)) THEN   !demote RRA
 
599
            RA(I) = RA(J)
 
600
            I = J
 
601
            J = J + J
 
602
         ELSE
 
603
            J = IR + 1
 
604
         ENDIF
 
605
         GOTO 20
 
606
      ENDIF
 
607
C 
 
608
      RA(I) = RRA                 !put RRA into its slot
 
609
      GOTO 10
 
610
C 
 
611
C 
 
612
1000  I = NH + 1
 
613
      DO 2000 L=I,N
 
614
         IF (RA(L).LT.RA(NH)) THEN
 
615
            RRA = RA(L)
 
616
            JL = 0
 
617
            JU = I
 
618
C 
 
619
1010        JM = (JU+JL) / 2
 
620
            IF (RRA.LT.RA(JM)) THEN
 
621
               JU = JM
 
622
            ELSE
 
623
               JL = JM
 
624
            ENDIF
 
625
            IF ((JU-JL).GT.1) GOTO 1010               !loop
 
626
C 
 
627
CC            DO 1100 M=L,JL+2,-1          !RA kept intact
 
628
            DO 1100 M=NH,JL+2,-1           !ojo - RA not kept intact that way!
 
629
               RA(M) = RA(M-1)             !but faster
 
630
1100        CONTINUE
 
631
            RA(JL+1) = RRA
 
632
         ENDIF
 
633
2000  CONTINUE
 
634
C 
 
635
      RETURN
 
636
      END
 
637
 
 
638
      SUBROUTINE XSORT1(RA,NDIM,NH,RINDX)
 
639
C 
 
640
C  RA  = real array of which the elements were used for sorting INDARR
 
641
C  NDIM  = full dimension of RA
 
642
C  NH = index of median - we only sort fully until NH
 
643
C  then we insert
 
644
C 
 
645
C  this is the Heapsort algorithm from "Numerical Recipes", page 233
 
646
C 
 
647
      IMPLICIT NONE
 
648
C 
 
649
      PARAMETER (IXXSIZ = 25000)
 
650
C 
 
651
      INTEGER   NDIM,NH,RINDX
 
652
      INTEGER   L,IR,J,I
 
653
      INTEGER   INDXT
 
654
C 
 
655
      INTEGER   INDARR(IXXSIZ),BACKAR(IXXSIZ),SAVARR(IXXSIZ)
 
656
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
 
657
C 
 
658
      REAL      RA(NDIM),RRA
 
659
C 
 
660
      DO 5, J=1,NDIM
 
661
         INDARR(J) = J
 
662
5     CONTINUE
 
663
C 
 
664
      L = NDIM/2 + 1
 
665
      IR = NDIM
 
666
C 
 
667
10    IF (L.GT.1) THEN            !still hiring
 
668
         L = L - 1
 
669
         INDXT = INDARR(L)
 
670
         RRA = RA(INDXT)
 
671
      ELSE                        !in retirement + promotion phase
 
672
         INDXT = INDARR(IR)
 
673
         RRA = RA(INDXT)          !clear a space at end of array
 
674
         INDARR(IR) = INDARR(1)   !retire the top of the heap into it
 
675
         IR = IR - 1              !decrease the size of the corporation
 
676
         IF (IR.EQ.1) THEN        !done with last promotion
 
677
            INDARR(1) = INDXT     !the least competent guy of all...
 
678
            GOTO 1000
 
679
         ENDIF
 
680
      ENDIF
 
681
C 
 
682
      I = L                 !in hiring as well as in promotion phase
 
683
      J = L + L             !we here set up to sift down RRA to its right level
 
684
C 
 
685
20    IF (J.LE.IR) THEN                    !"Do while J<=IR"
 
686
         IF (J.LT.IR) THEN
 
687
            IF (RA(INDARR(J)).LT.RA(INDARR(J+1)))
 
688
     +         J = J + 1                   !compare to the better underling
 
689
         ENDIF
 
690
         IF (RRA.LT.RA(INDARR(J))) THEN    !demote RRA
 
691
            INDARR(I) = INDARR(J)
 
692
            I = J
 
693
            J = J + J
 
694
         ELSE
 
695
            J = IR + 1
 
696
         ENDIF
 
697
         GOTO 20
 
698
      ENDIF
 
699
C 
 
700
      INDARR(I) = INDXT                 !put RRA into its slot
 
701
      GOTO 10
 
702
C
 
703
1000  DO 4000, J=1,NDIM                 !set up the back index array
 
704
         BACKAR(INDARR(J)) = J
 
705
4000  CONTINUE
 
706
C 
 
707
      RINDX = INDARR(NH)
 
708
C 
 
709
      RETURN
 
710
      END
 
711
 
 
712
      SUBROUTINE XSORT2(RA,NDIM,NH,JWOFF,NEWRA,XSIZE,NCOUNT,RINDX)
 
713
C 
 
714
C  RA  = real array of which the elements were used for sorting INDARR
 
715
C  NDIM  = full dimension of RA
 
716
C  NH = index of median - we only sorted fully until NH then we inserted
 
717
C  IOFF = index of RA element in question
 
718
C  INDARR = semi sorted index array from last sorting of RA
 
719
C 
 
720
C 
 
721
      IMPLICIT NONE
 
722
C 
 
723
      PARAMETER (IXXSIZ = 25000)
 
724
C 
 
725
      INTEGER   NDIM,NH,JWOFF,XSIZE,NCOUNT,RINDX
 
726
      INTEGER   IOFF,INDXT,JL,JU,JM,M,N
 
727
C 
 
728
      INTEGER   INDARR(IXXSIZ),BACKAR(IXXSIZ),SAVARR(IXXSIZ)
 
729
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
 
730
C 
 
731
      REAL      RA(NDIM),NEWRA(NCOUNT)
 
732
      REAL      NEWR,QR
 
733
C 
 
734
      IOFF = JWOFF
 
735
C 
 
736
      DO 5000, N=1,NCOUNT
 
737
C 
 
738
         NEWR = NEWRA(N) 
 
739
         RA(IOFF) = NEWR               !new value
 
740
         INDXT = BACKAR(IOFF)
 
741
C 
 
742
         IF (INDXT.EQ.1) THEN
 
743
            IF (NEWR.LE.RA(INDARR(2)))             !it stays very first
 
744
     +         GOTO 4900
 
745
         ELSE
 
746
            IF (NEWR.LE.RA(INDARR(1))) THEN        !it's very first
 
747
               DO 1000, M=INDXT,2,-1               !remove old entry
 
748
                  INDARR(M) = INDARR(M-1)          !by shifting up
 
749
                  BACKAR(INDARR(M)) = M
 
750
1000           CONTINUE
 
751
               INDARR(1) = IOFF 
 
752
               BACKAR(IOFF) = 1
 
753
               GOTO 4900
 
754
            ENDIF
 
755
         ENDIF
 
756
C 
 
757
         IF (INDXT.EQ.NDIM) THEN
 
758
            IF (NEWR.GE.RA(INDARR(NDIM-1)))        !it stays very last
 
759
     +         GOTO 4900
 
760
         ELSE
 
761
            IF (NEWR.GE.RA(INDARR(NDIM))) THEN     !it very last
 
762
               DO 2000, M=INDXT,NDIM-1             !remove old entry
 
763
                  INDARR(M) = INDARR(M+1)          !by shifting down
 
764
                  BACKAR(INDARR(M)) = M
 
765
2000           CONTINUE
 
766
               INDARR(NDIM) = IOFF
 
767
               BACKAR(IOFF) = NDIM
 
768
               GOTO 4900
 
769
            ENDIF
 
770
         ENDIF
 
771
C 
 
772
         IF (INDXT.EQ.1) THEN
 
773
            QR = NEWR - 1.0                        !so we satisy following test
 
774
         ELSE
 
775
            QR = RA(INDARR(INDXT-1))
 
776
         ENDIF
 
777
C 
 
778
         IF (NEWR.GT.QR) THEN                      !INDXT = 1 also done here...
 
779
            IF (NEWR.LE.RA(INDARR(INDXT+1))) GOTO 4900
 
780
            JL = INDXT
 
781
            JU = NDIM + 1
 
782
C
 
783
2010        JM = (JU+JL) / 2
 
784
            IF (NEWR.LT.RA(INDARR(JM))) THEN
 
785
               JU = JM
 
786
            ELSE
 
787
               JL = JM
 
788
            ENDIF
 
789
            IF ((JU-JL).GT.1) GOTO 2010          !loop
 
790
C
 
791
            DO 2100, M=INDXT,JL-1
 
792
               INDARR(M) = INDARR(M+1)
 
793
               BACKAR(INDARR(M)) = M
 
794
2100        CONTINUE
 
795
            INDARR(JL) = IOFF
 
796
            BACKAR(IOFF) = JL 
 
797
C            
 
798
         ELSE IF (NEWR.LT.QR) THEN
 
799
            JL = 0
 
800
            JU = INDXT
 
801
C
 
802
3010        JM = (JU+JL) / 2
 
803
            IF (NEWR.LT.RA(INDARR(JM))) THEN
 
804
               JU = JM
 
805
            ELSE
 
806
               JL = JM
 
807
            ENDIF
 
808
            IF ((JU-JL).GT.1) GOTO 3010          !loop
 
809
C
 
810
            DO 3100, M=INDXT,JL+2,-1
 
811
               INDARR(M) = INDARR(M-1)
 
812
               BACKAR(INDARR(M)) = M
 
813
3100        CONTINUE
 
814
            INDARR(JL+1) = IOFF
 
815
            BACKAR(IOFF) = JL + 1
 
816
         ENDIF
 
817
C 
 
818
4900     IOFF = IOFF + XSIZE
 
819
5000  CONTINUE
 
820
C 
 
821
      RINDX = INDARR(NH)
 
822
C
 
823
      RETURN
 
824
      END
 
825
 
 
826
      SUBROUTINE XSAVX(DIRFLG,NDIM)
 
827
C 
 
828
      INTEGER   DIRFLG, NDIM
 
829
      INTEGER   N,TEST
 
830
C
 
831
      PARAMETER (IXXSIZ = 25000)
 
832
C 
 
833
      INTEGER   INDARR(IXXSIZ),BACKAR(IXXSIZ),SAVARR(IXXSIZ)
 
834
      COMMON    /TMPARR/  INDARR,BACKAR,SAVARR
 
835
C 
 
836
      IF (DIRFLG.EQ.1) THEN
 
837
         DO 100, N=1,NDIM
 
838
            TEST = INDARR(N)
 
839
            SAVARR(N) = TEST
 
840
100      CONTINUE
 
841
      ELSE
 
842
C 
 
843
         DO 500, N=1,NDIM
 
844
            TEST = SAVARR(N)
 
845
            INDARR(N) = TEST
 
846
500      CONTINUE
 
847
         DO 600, N=1,NDIM                       !rebuild back index array
 
848
            BACKAR(INDARR(N)) = N
 
849
600      CONTINUE
 
850
      ENDIF
 
851
C 
 
852
      RETURN
 
853
      END