1
C @(#)rebdec.for 19.1 (ESO-DMD) 02/25/03 13:31:23
2
C===========================================================================
3
C Copyright (C) 1995 European Southern Observatory (ESO)
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.
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.
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,
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
27
C===========================================================================
31
C++++++++++++++++++++++++++++++++++++++++++++++++++
34
C program REBDEC version 1.00 880930
35
C Basic algorithm: L. Lucy ESO - Garching
36
C MIDAS implementation: D. Baade ESO - Garching
37
C K. Banse 1.10 890517
41
C (nonlinear) rebinning, deconvolution, point spread function, undersampling
44
C Rebin an image to smaller stepsize to remove effects of coarse PSF sampling
45
C ("pixelation"). Optionally deconvolve with PSD.
48
C Get the names of input + output frames from IN_A + OUT_A (both 2-D)
49
C Do the rebinning (which is nonlinear in the way the flux is redistributed
50
C whereas the spatial coordinate transformation is strictly linear) while
51
C simultaneously deconvolving with user-supplied point spread function.
54
C the following keys are used:
56
C IN_A/C/1/60 input frame
57
C IN_B/C/1/60 point spread function
58
C OUT_A/C/1/60 output frame
59
C INPUTI/I/1/2 zoom factors in x and y
60
C INPUTI/I/3/1 number of iterations
63
C 1.10 convert to Portable Midas
64
C 1.20 fix bug in STIPUT (missing NAXIS..)
68
C --------------------------------------------------------------------
72
INTEGER STAT,ZOOM(2),NITER,NAXIN,NPXIN(2)
73
INTEGER INNO,PSFNO,OUTNO,DUMNO,D0NO,D1NO,D2NO,D3NO,D4NO
74
INTEGER NAXPSF,NPXPSF(2),IAV,OFF(2)
75
INTEGER NPXPSR(4),NAXPSR,I,SUBLO(2),GS
76
INTEGER NPXOUT(2),RESPIX(2),SIZE,NAXOUT,NPIXL(2)
77
INTEGER NPIXZ(2),MSIZ,INPUTI(4)
78
INTEGER UNI(1),NULO,MADRID(1)
79
INTEGER*8 INPTR,PSFPTR,OUTPTR,DUMPTR,D0PTR,D1PTR,D2PTR,
82
CHARACTER INFRM*60,OUTFRM*60,PSF*60
83
CHARACTER CUNIN*64,CUNPSF*64,IDENT*72
84
CHARACTER IDENTP*72,CUNOUT*64,IDENTO*72
86
DOUBLE PRECISION STAIN(2),STPIN(2),STAPSF(2)
87
DOUBLE PRECISION STPPSF(2),STPOUT(2),STPPSR(4)
88
DOUBLE PRECISION STAPSR(4),STAOUT(2),TEST,EPS
90
REAL CUTS(2),CUTVALS(4)
92
INCLUDE 'MID_INCLUDE:ST_DEF.INC'
96
DATA CUTVALS /4*0./, SUBLO /1,1/
98
INCLUDE 'MID_INCLUDE:ST_DAT.INC'
100
C set up MIDAS environment + enable automatic error abort
101
CALL STSPRO('REBDEC')
103
C get name of input image + map frame
104
CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT)
105
CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
106
+ 2,NAXIN,NPXIN,STAIN,STPIN,
107
+ IDENT,CUNIN,INPTR,INNO,STAT)
108
IF (NAXIN.NE.2) CALL STETER
109
+ (1,'ERROR: Currently only 2-D frames supported!')
111
C get and check zoom factors in x + y ... and the number of iterations
112
CALL STKRDI('INPUTI',1,4,IAV,INPUTI,UNI,NULO,STAT)
118
IF ((((ZOOM(1)*ZOOM(2))/2)*2).EQ.(ZOOM(1)*ZOOM(2)))
120
+ (2,'ERROR: Both zoom factors need to be odd!')
121
IF ((ZOOM(1).LT.1).OR.(ZOOM(2).LT.1))
123
+ (3,'ERROR: Zoom factors less than unity are illegal!')
125
IF (NITER.LT.1) CALL STETER
126
+ (4,'ERROR: Number of iterations must at least be 1!')
128
C get name of point spread function frame + map it
129
CALL STKRDC('IN_B',1,1,60,IAV,PSF,UNI,NULO,STAT)
130
CALL STIGET(PSF,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
131
+ 2,NAXPSF,NPXPSF,STAPSF,
132
+ STPPSF,IDENTP,CUNPSF,PSFPTR,PSFNO,STAT)
134
C Further error traps: make sure dimensions are the same ...
137
+ (5,'Dimensions of inframe and PSF must be the same!')
139
C ... steps have the right relative proportions,
140
IF ((NINT(STPIN(1)/STPPSF(1)).NE.ZOOM(1)).OR.
141
+ (NINT(STPIN(2)/STPPSF(2)).NE.ZOOM(2)))
143
+ (6,'Steps of inframe and PSF must differ by zoom fact!')
145
C ... the zoom factor is not even when the number of pixels in the input PSF
146
C is odd (check both X and Y)
148
IF ((NPXPSF(I)/2)*2.EQ.NPXPSF(I))
150
+ (7,'NPIX(PSF) must be odd in both X and Y!')
153
C Get name of output frame, ...
154
CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT)
155
IF (OUTFRM.EQ.INFRM) ! in- and output frame must be different
157
+ (8,'Names of inframe + outframe must be different!')
159
C ... depending on its existence or not,
161
C ... open it for reading and writing or create it for writing to it
163
IF (GS.EQ.0) THEN ! we start without an initial guess
165
C Calculate further defining parameters of output frame ...
166
STPOUT(1)=STPIN(1)/ZOOM(1)
167
STPOUT(2)=STPIN(2)/ZOOM(2)
168
NPXOUT(1)=NPXIN(1)*ZOOM(1)
169
NPXOUT(2)=NPXIN(2)*ZOOM(2)
172
STAOUT(1)=STAIN(1)-STPOUT(1)*(ZOOM(1)/2)
173
STAOUT(2)=STAIN(2)-STPOUT(2)*(ZOOM(2)/2)
176
IDENTO(1:) = IDENT(1:)
177
CUNOUT(1:) = CUNIN(1:)
179
CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
180
+ NAXOUT,NPXOUT,STAOUT,
181
+ STPOUT,IDENT,CUNIN,OUTPTR,OUTNO,STAT)
183
ELSE ! there is a guess
185
C If there is a guess, we first map the frame ...
186
CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
187
+ 2,NAXOUT,NPXOUT,STAOUT,
188
+ STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT)
190
C ... and then check if it is technically compatible with the
191
C input frame and the zoom factors:
193
IF (NAXOUT.NE.2) CALL STETER
194
+ (9,'Output frame, too, must be 2-dimensional!')
195
IF ((NPXIN(1)*NPXIN(2)*ZOOM(1)*ZOOM(2)).NE.
196
+ (NPXOUT(1)*NPXOUT(2)))
198
+ (10,'NPIX(OUT) incompatible with NPIX(IN) + zoom factors!')
199
TEST = (STAOUT(1)+STPOUT(1)*(ZOOM(1)/2))*
200
+ (STAOUT(2)+STPOUT(2)*(ZOOM(2)/2))
201
EPS = STPIN(1)*STPIN(2)*0.01
202
IF (ABS(TEST-(STAIN(1)*STAIN(2))).GT.EPS)
204
+ (11,'START values of in- and outframe differ!')
205
EPS = STPIN(1)*STPIN(2)*0.0001
206
IF (ABS((STPOUT(1)*ZOOM(1)*STPOUT(2)*ZOOM(2))-
207
+ (STPIN(1)*STPIN(2))).GT.EPS)
209
+ (12,'steps of in- and outframe must differ by factor ZOOM!')
210
ENDIF ! guess or no guess?
212
C We also need six working buffers for ...
214
C a) the output frame with extrapolation beyond its edges
217
C The radius of the PSF determines the width of the margin to be added
218
C around inframe, for outframe additionally the zoom factor must be considered:
219
OFF(I) = NPXPSF(I)/(2*ZOOM(I)) + 2
220
NPIXL(I) = NPXIN(I)+2*OFF(I)
221
NPIXZ(I) = NPIXL(I)*ZOOM(I)
224
MSIZ = NPIXZ(1) * NPIXZ(2)
225
CALL STFCRE('DUMMY',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
227
CALL STFMAP(DUMNO,F_X_MODE,1,MSIZ,IAV,DUMPTR,STAT)
229
C If there is a guess, we transfer it right away to the enlarged frame:
231
+ CALL EXTND(MADRID(OUTPTR),MADRID(DUMPTR),
232
+ NPXOUT(1),NPXOUT(2),NPIXZ(1),NPIXZ(2),OFF,
235
C In order not to waste virtual memory space, we temporarily close OUTFRM
236
CALL STFCLO(OUTNO,STAT)
238
C b) the input frame after extrapolation beyond its edges
239
MSIZ = NPIXL(1) * NPIXL(2)
240
CALL STFCRE('DUMMY0',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
242
CALL STFMAP(D0NO,F_X_MODE,1,MSIZ,IAV,D0PTR,STAT)
244
C We do the enlargement here and then close inframe (no longer needed)
246
CALL EXTND (MADRID(INPTR),MADRID(D0PTR),NPXIN(1),NPXIN(2),
247
+ NPIXL(1),NPIXL(2),OFF,1,1)
248
CALL STFCLO(INNO,STAT)
250
C c) ratio of observed to calculated image
251
MSIZ = NPIXL(1) * NPIXL(2)
252
CALL STFCRE('DUMMY1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
254
CALL STFMAP(D1NO,F_X_MODE,1,MSIZ,IAV,D1PTR,STAT)
256
C d) correction factors to be applied to previous intermediate result
257
MSIZ = NPIXZ(1) * NPIXZ(2)
258
CALL STFCRE('DUMMY2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
260
CALL STFMAP(D2NO,F_X_MODE,1,MSIZ,IAV,D2PTR,STAT)
262
C e) re-calculated frame on coarse grid
263
MSIZ = NPIXL(1) * NPIXL(2)
264
CALL STFCRE('DUMMY3',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
266
CALL STFMAP(D3NO,F_X_MODE,1,MSIZ,IAV,D3PTR,STAT)
268
C f) PSF, but re-sampled on coarse grid for various centerings.
270
C First, the dimension of this new PSF need to be properly determined:
274
C The number of pixels of the coarse PSF (see subroutine PSFRESAMP) must be
275
C odd. If ZOOM divides NPIX(PSF), this is automatically fulfilled because
276
C both numers are odd.
277
NPXPSR(2*I-1) = NPXPSF(I) / ZOOM(I)
279
C Otherwise, we just make it odd:
280
IF (((NPXPSF(I)/ZOOM(I))*ZOOM(I)).NE.NPXPSF(I))
281
+ NPXPSR(2*I-1) = (NPXPSR(2*I-1)/2)*2 + 1
283
C Finally, an additional margin of 1 is added at either limit in order
284
C to avoid truncations of input-PSF:
285
NPXPSR(2*I-1) = NPXPSR(2*I-1)+2
286
NPXPSR(2*I) = ZOOM(I)
287
STPPSR(2*I-1) = STPPSF(I)
288
STPPSR(2*I) = 1.D0/ZOOM(I)
289
STAPSR(2*I-1) = STAPSF(I)
290
STAPSR(2*I) = 1.D0/ZOOM(I)
293
MSIZ = NPXPSR(1) * NPXPSR(2) * NPXPSR(3) * NPXPSR(4)
294
CALL STFCRE('DUMMY4',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
296
CALL STFMAP(D4NO,F_X_MODE,1,MSIZ,IAV,D4PTR,STAT)
298
C Do the resampling of the PSF here and then close the input PSF
299
C because it's no longer used
300
CALL PSFRES(NPXPSF(1),NPXPSF(2),NPXPSR(1),NPXPSR(3),
301
+ ZOOM(1),ZOOM(2),MADRID(PSFPTR),MADRID(D4PTR))
302
CALL STFCLO(PSFNO,STAT)
304
C Now comes the real thing:
305
CALL REBDC(NPIXL(1),NPIXL(2),NPIXZ(1),NPIXZ(2),NITER,
306
+ NPXPSR(1),NPXPSR(3),ZOOM(1),ZOOM(2),GS,
307
+ MADRID(D0PTR),MADRID(D1PTR),MADRID(D2PTR),
308
+ MADRID(D3PTR),MADRID(D4PTR),MADRID(DUMPTR))
310
C Close dispensible intermediate files and re-open outframe . Cut off
311
C extrapolated regions of intermediate result as to match size of OUTFRM,
312
C insert into outframe:
314
CALL STFCLO(D0NO,STAT)
315
CALL STFCLO(D1NO,STAT)
316
CALL STFCLO(D2NO,STAT)
317
CALL STFCLO(D3NO,STAT)
318
CALL STFCLO(D4NO,STAT)
320
CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
321
+ 2,NAXOUT,NPXOUT,STAOUT,
322
+ STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT)
323
CALL TRUNC(MADRID(DUMPTR),MADRID(OUTPTR),NPIXZ(1),NPIXZ(2),
324
+ NPXOUT(1),NPXOUT(2),OFF,ZOOM(1),ZOOM(2))
326
C Last, we determine physical minimum and maximum of OUTFRM
327
CALL STVALS('MIN',MADRID(OUTPTR),2,NPXOUT,SUBLO,
328
+ NPXOUT,CUTVALS,CUTS,RESPIX,SIZE,STAT)
330
C Copy descriptors and update LHCUTS + HISTORY
331
CALL DSCUPT(INNO,OUTNO,' ',STAT)
332
CALL STDWRR(OUTNO,'LHCUTS',CUTS,3,2,UNI,STAT)
338
SUBROUTINE REBDC(ICX,ICY,I1X,I1Y,NITER,OPSFX,OPSFY,ZX,ZY,GS,
339
1 PHT,RT,CF,PHI,PSF,PSI)
341
C Removal of pixel grid from CCD data plus seeing by appropriate
342
C rebinning (including 2-D deconvolution) to smaller step size.
343
C The deconvolution algorithm follows the description given in
344
C The Astronomical Journal, Vol. 79, p. 745-754 (1974) but the
345
C rebinning to a smaller step is an extra feature.
347
C L.B. Lucy ESO - Garching 1988 August
348
C Modifications applied for MIDAS version (incl. generalisation
349
C to non-square images and PSF's):
351
C D. Baade ESO - Garching 880930
355
INTEGER ICX,ICY,GS,ZX,ZY,NITER,STAT,I1X,I1Y,IT,I,J
358
REAL PHT(ICX,ICY),PHI(ICX,ICY),CD,CD1,RT(ICX,ICY)
359
REAL PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),PSI(I1X,I1Y)
361
CHARACTER TIME*25,CBUF*80
363
C PHT inframe, but extrapolated beyond edges, dims: ICX x ICY
364
C PHI calculated frame in coarse grid
365
C RT basically, the ratio observed/calculated in coarse grid
366
C CF correction factor to be applied to current fine-grid result
367
C PSI fine-grid output frame, STAINg from some trivial guess
368
C IC dimensions of input frame: ICX x ICY
369
C Z zoom factors coarse -> fine grid: ZX, ZY
370
C I1 dimensions of output frame: ZX*ICX x ZY*ICY
373
C ... If an initial guess has not been provided, we start with a
374
C ... smooth trial-image
381
PSI(I,J) = CD1 ! initial smooth trial - image
384
ENDIF ! start frame prepared
386
DO 1000 IT=1,NITER ! iteration loop
388
C Form PSF-mapped coarse-grid image of guess PSI (result in frame PHI):
389
CALL INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
392
C Compute ratio of observed to calculated (guessed) image
396
RT(I,J) = 1.0 ! initialize
398
+ RT(I,J) = CD*PHT(I,J)/PHI(I,J) ! obs/calc
402
C Determine correction factor to current intermediate result
403
CALL INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
406
C Apply the correction, thereby completing this iteration
409
PSI(I,J) = PSI(I,J)*CF(I,J) ! apply correc factor
413
C The code is slow. Progress messages seem a good idea:
415
WRITE (CBUF,10000) TIME,IT
416
CALL STTPUT(CBUF,STAT)
422
10000 FORMAT (A,' ',I2,'. iteration completed')
425
SUBROUTINE EXTND (RAW,PHT,IRX,IRY,ICX,ICY,OFF,ZX,ZY)
427
C To avoid ugly effects at the frame limits, repeat the outermost
428
C columns (rows) RX (RY) times, i.e. by the halfwidth of the
429
C PSF resampled on the coarse grid.
433
INTEGER IRX,IRY,ICX,ICY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY
435
REAL RAW(IRX,IRY),PHT(ICX,ICY)
440
C Insert raw frame at center of output frame
446
PHT(IQ,JQ) = RAW(I,J)
450
C ... Do the extrapolation, first in X:
454
PHT(I,J+OFY) = RAW(1,J)
455
PHT(IRX+OFX+I,J+OFY) = RAW(IRX,J)
458
C Here we also deal with the four corners
462
PHT(IRX+OFX+I,J) = RAW(IRX,1)
463
PHT(I,IRY+OFY+J) = RAW(1,IRY)
464
PHT(IRX+OFX+I,IRY+OFY+J) = RAW(IRX,IRY)
468
C The continuation in Y:
472
PHT(I+OFX,J) = RAW(I,1)
473
PHT(I+OFX,IRY+OFY+J) = RAW(I,IRY)
480
SUBROUTINE PSFRES (IPSFX,IPSFY,OPSFX,OPSFY,ZX,ZY,PSFIN,PSF)
482
C Author: L. Lucy ESO - Garching 1988 August
483
C Modifications: D. Baade ESO - Garching 880930
485
C The subroutine resamples the user-supplied high-resolution PSF
486
C on the coarser grid of the image to be deconvolved. The array
487
C holding the resulting PSF is four-dimensional because the resampling
488
C is done ZX x ZY times where each time the input PSF is
489
C shifted by 1/ZX and/or 1/ZY as to cover all cases of the
490
C relative positioning of the two PSF grids.
494
INTEGER II,JJ,IP,JP,K,L,ZX,ZY,INDX,INDY,IPSFX,IPSFY
495
INTEGER OPSFX,OPSFY,DX,DY,OFFX,OFFY,CX,CY
497
REAL PSFIN(IPSFX,IPSFY),PSFBUF,CD,PSF(OPSFX,ZX,OPSFY,ZY)
500
DX = (OPSFX*ZX-IPSFX)/2 ! the excess in width of the resampled PSF
501
DY = (OPSFY*ZY-IPSFY)/2 ! is always an even number of small pixels.
503
C Shift by one small pixel in X ...
506
OFFX = II-(ZX+1)/2+DX
510
OFFY = JJ-(ZY+1)/2+DY
512
C For all coarse-grid pixels ...
517
PSFBUF = 0. ! initialize
519
C co-add all fine-grid pixels covered by them at the current offsets
523
C Skip fine-grid pixels for which the input-PSF is undefined:
524
IF ( (INDX.LT.1) .OR. (INDX.GT.IPSFX) ) GOTO 400
527
IF ( (INDY.GE.1) .AND. (INDY.LE.IPSFY) )
528
+ PSFBUF = PSFBUF + PSFIN(INDX,INDY)
532
PSF(IP,II,JP,JJ) = PSFBUF*CD
542
SUBROUTINE INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,
545
C Author: L. Lucy ESO - Garching 1988 August
546
C Modifications: D. Baade ESO - Garching 880930
548
C This subroutine maps the current fine-grid guess, PSI,
549
C on to the coarse grid (result in array PHI). The mapping
550
C uses the stack of PSF's that have been constructed by
551
C re-sampling the input PSF onto the coarse grid, i.e. basically
552
C a convolution with the psf is done.
556
INTEGER ICX,ICY,I1X,I1Y,OPSFX,OPSFY,IMIN,IMAX,JMIN,JMAX
557
INTEGER ID,ID0,ID1,ID2,II,IP,JD,JD0,JD1,JD2,JP,JJ,ZX,ZY
558
INTEGER I,J,CX(10000),CY(10000),RX,RY,RX1,RY1,RX1I,RY1J
559
INTEGER AUXX(2000),AUXY(2000)
561
REAL PSF(OPSFX,ZX,OPSFY,ZY),PSI(I1X,I1Y),PHIBUF,PHI(ICX,ICY)
563
RX = OPSFX/2 ! 'radius' of re-sampled PSF,
564
RY = OPSFY/2 ! full diameter is 2*R+1
568
C For all fine pixels, determine X and Y coordinates on the coarse grid.
569
C Also, introduce auxiliary array to further save computing time.
573
AUXX(I) = (CX(I)-1)*ZX
577
AUXY(J) = (CY(J)-1)*ZY
580
C Proceed pixel by pixel on the coarse grid
589
DO 2000 I=IMIN,IMAX ! pixel of coarse grid
593
DO 1800 J=JMIN,JMAX ! pixel of coarse grid
598
C Then sum over all relevant pixels in fine grid. the PSF serves as
599
C the weighting function. For a small pixel (II,JJ) within a given large
600
C pixel (I,J) the value PSF(IP,II,JP,JJ) is taken. IP and JP point
601
C to the coarse pixels contributing (via the PSF) to the currently
602
C treated coarse pixel. IP and JP are no absolute coordinates but
603
C rather are the distances of the contributing large pixels from the
604
C large pixel whose flux is just being computed. II and JJ are also
605
C relative indices and only range from 1 to ZX and ZY, respectively.
607
PHIBUF = 0. ! initialize
609
DO 1600 ID=ID1,ID2 ! relevant range in fine grid
610
II=ID-AUXX(ID) ! position within large pixel
611
IP=CX(ID)+RX1I ! position of large pixel in PSF
612
DO 1500 JD=JD1,JD2 ! relevant range in fine grid
613
JJ = JD-AUXY(JD) ! position within large pixel
614
JP = CY(JD)+RY1J ! position of large pixel in PSF
615
PHIBUF = PHIBUF+PSF(IP,II,JP,JJ)*PSI(ID,JD) ! off we go
618
PHI(I,J) = PHIBUF ! this coarse bin finished
623
C Now extrapolate solution achieved so far to edges of frame by
624
C constant extrapolation
629
PHI(IMAX+I,J)=PHI(IMAX,J)
632
PHI(I,J)=PHI(IMIN,JMIN)
633
PHI(IMAX+I,J)=PHI(IMAX,JMIN)
634
PHI(I,JMAX+J)=PHI(IMIN,JMAX)
635
PHI(IMAX+I,JMAX+J)=PHI(IMAX,JMAX)
641
PHI(I,JMAX+J)=PHI(I,JMAX)
648
SUBROUTINE INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,PSF,RT,CF)
650
C Author: L. Lucy ESO - Garching 1988 August
651
C Modifications: D. Baade ESO - Garching 880930
653
C Here the correction factors (pixelwise on the fine grid, array CF)
654
C are computed as the PSF-mapped ratio (pixelwise on the coarse grid,
655
C array RT) of the observed image to the result of the last iteration
656
C after its re-sampling on the coarse grid.
660
INTEGER I,J,I1X,I1Y,ICX,ICY,II,ID1,ID2,IMIN,IMAX,JMIN,JMAX
661
INTEGER JJ,JD1,JD2,ID,IP,JD,JP,OPSFX,OPSFY,ZX,ZY,RX,RY
662
INTEGER CX(10000),CY(10000),AUXX1(10000),AUXY1(10000)
663
INTEGER AUXX2(10000),AUXY2(10000)
665
REAL PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),RT(ICX,ICY)
668
RX=OPSFX/2 ! 'radius' of coarse PSF,
669
RY=OPSFY/2 ! full diameter is 2*R+1
676
C Determine for all fine pixels coordinates on coarse grid (taking
677
C into account that the input frame has been enlarged).
678
C Also, introduce auxiliary arrays to further reduce computing time.
681
CX(I)=(I-1)/ZX+1 ! X-coordinates on coarse grid
683
AUXX2(I)=(CX(I)-1)*ZX
686
CY(J)=(J-1)/ZY+1 ! Y-coordinates on coarse grid
688
AUXY2(J)=(CY(J)-1)*ZY
691
DO 1000 I=IMIN,IMAX ! pixel on fine grid
692
II=I-AUXX2(I) ! pixel-# within coarse pixel
696
DO 900 J=JMIN,JMAX ! pixel on fine grid
697
JJ=J-AUXY2(J) ! pixel-# within coarse pixel
701
CFBUF = 0. ! initialize
703
DO 800 ID=ID1,ID2 ! relevant coarse-grid pixels
704
IP=AUXX1(I)+ID ! X-position within coarse PSF
705
DO 700 JD=JD1,JD2 ! relevant coarse-grid pixels
706
JP=AUXY1(J)+JD ! Y-position within coarse PSF
707
CFBUF=CFBUF+PSF(IP,II,JP,JJ)*RT(ID,JD) ! the correction
712
C At this point the corrections for fine-grid pixel (I,J) to the
713
C ratios at positions (ID,JD) on the coarse grid have been fully
714
C evaluated for all pixels (ID,JD) concerned by pixel (I,J).
719
C Fill up rest of frame by constant extrapolation of the edges
724
CF(IMAX+I,J)=CF(IMAX,J)
727
CF(I,J)=CF(IMIN,JMIN)
728
CF(IMAX+I,J)=CF(IMAX,JMIN)
729
CF(I,JMAX+J)=CF(IMIN,JMAX)
730
CF(IMAX+I,JMAX+J)=CF(IMAX,JMAX)
736
CF(I,JMAX+J)=CF(I,JMAX)
743
SUBROUTINE TRUNC(PSI,OUT,I1X,I1Y,IOX,IOY,OFF,ZX,ZY)
745
C This subroutine truncates the result frame such as to correspond
746
C to the true dimensions (in world coordinates) of the input frame.
750
INTEGER I1X,I1Y,IOX,IOY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY
752
REAL PSI(I1X,I1Y),OUT(IOX,IOY)
761
OUT(I,J) = PSI(IQ,JQ)
added
removed
contrib/surfphot/src/rebdec.for
