v3.0.3 - 11 Mar 2016
C$ Procedure
SUBROUTINE FIND_LAMBDA_PIC(NTMP,QSZ,DNX, SP,CP,TUSE,TMPL,
. LAMBDA,PHI,CHI,CHI0)
C$ Abstract
C
C The predicted brightness for a pixel in a maplet is given by
C
C B = LAMBDA*ETA*ILLUM(COS(I),COS(E),PHASE)+PHI
C
C where I and E are incidence and emission angles, ETA is relative
C albedo, normalized to 1 across the maplet. Since no attempt is made
C to account for absolute radiometry, there is a constant lambda that
C serves to relate the topography manifested by the slopes leading to
C variations in I and E to the brightness measured by the imaging data.
C Clearly, a smaller value of lambda requires a more robust topography
C than a larger one, so we must have some initial idea of the topography
C to determine lambda. This is accomplished by making large, low
C resolution maplets to start, that owe most of their relief to the
C curvature of the body from ordinary stereography. Once lambda is
C determined, it sets the scale for the higher frequency topography in
C the maplet. This topography, in turn, sets the scale for the higher
C resolution maplets that overlap the original maplet and their even
C higher frequency topography. The topography scale of any maplet is
C constantly refined by contributions such as differential stereo, limb
C contributions and randomly selected heights from overlapping maplets.
C
C This subroutine estimates the value oflambda in the brightness function
C from the observed brightness distribution and the nominal topography. It
C can also solve for a constant background in the brightness such as was
C observed due to atmospheric scattering in Viking Orbiter images of Mars.
C This latter PHI term is usually ignored by setting the parameter CHI0 = 1
C in INIT_LITHOS.TXT.
C$ Disclaimer
C
C
C$ Required_Reading
C
C R.W. Gaskell, et.al, "Characterizing and navigating small bodies
C with imaging data", Meteoritics & Planetary Science 43,
C Nr 6, 1049-1061 (2008)
C
C$ Declarations
IMPLICIT NONE
INTEGER NTMP
DOUBLE PRECISION SP(3)
DOUBLE PRECISION CP(3)
DOUBLE PRECISION ALPHA
DOUBLE PRECISION RPD
DOUBLE PRECISION VDOT
DOUBLE PRECISION ILLUM
DOUBLE PRECISION GAMMA
DOUBLE PRECISION ETA
DOUBLE PRECISION LAMBDA
DOUBLE PRECISION PHI
DOUBLE PRECISION CHI
DOUBLE PRECISION CHI0
DOUBLE PRECISION X1
DOUBLE PRECISION X2
DOUBLE PRECISION Z0
DOUBLE PRECISION Z1
DOUBLE PRECISION Z2
DOUBLE PRECISION ZI
DOUBLE PRECISION ZB
DOUBLE PRECISION ZBI
DOUBLE PRECISION ZBB
DOUBLE PRECISION ZII
REAL*4 DNX(-NTMP:NTMP,-NTMP:NTMP)
REAL*4 TMPL(-NTMP:NTMP,-NTMP:NTMP,3)
INTEGER QSZ
INTEGER I
INTEGER J
INTEGER N1
LOGICAL TUSE(-NTMP:NTMP,-NTMP:NTMP)
C$ Variable_I/O
C
C Variable I/O Description
C -------- --- --------------------------------------------------
C NTMP
C QSZ
C DNX
C SP
C CP
C TUSE
C TMPL
C LAMBDA
C PHI
C CHI
C CHI0
C
C$ File_I/O
C
C Filename I/O Description
C ---------------------------- --- -------------------------------
C None
C
C$ Restrictions
C
C
C$ Software_Documentation
C
C OSIRIS-REx Stereophotoclinometry Software Design Document
C OSIRIS-REx Stereophotoclinometry Software User's Guide
C
C$ Author_and_Institution
C
C R.W. Gaskell (PSI)
C
C$ Version
C
C
C
C$ SPC_functions_called
C ILLUM
C
C$ SPC_subroutines_called
C None
C
C$ SPICELIB_functions_called
C RPD
C VDOT
C
C$ SPICELIB_subroutines_called
C None
C
C$ Called_by_SPC_Programs
C AUTOREGISTER
C SUBROUTINE FIND_LAMBDA
LAMBDA=1
PHI=0
Z0=1
ZI=0
ZB=0
ZBI=0
ZBB=0
ZII=0
N1=0
% Do per pixel
DO I=-QSZ,QSZ
DO J=-QSZ,QSZ
IF((DNX(I,J).GT.0).AND.TUSE(I,J)) THEN
X1=DNX(I,J)
Z0=MIN(X1,Z0)
GAMMA=SQRT(1+TMPL(I,J,1)**2+TMPL(I,J,2)**2)
Z1=(SP(3)+TMPL(I,J,1)*SP(1)+TMPL(I,J,2)*SP(2))/GAMMA
Z2=(CP(3)+TMPL(I,J,1)*CP(1)+TMPL(I,J,2)*CP(2))/GAMMA
ETA=1+TMPL(I,J,3)
ALPHA=ACOS(VDOT(CP,SP))/RPD() % phase
X2=ETA*ILLUM(Z1,Z2,ALPHA)
ZI=ZI+X1
ZB=ZB+X2
ZBI=ZBI+X1*X2
ZBB=ZBB+X2**2
ZII=ZII+X1**2
N1=N1+1
ENDIF
ENDDO
ENDDO
IF((N1.GT.0).AND.(ZB.GT.0)) THEN
ZI=ZI/N1
ZB=ZB/N1
ZBI=ZBI/N1
ZBB=ZBB/N1
ZII=ZII/N1
LAMBDA=ZI/ZB
PHI=0
CHI=0
IF((ZBB-ZB**2.GT.0.001).AND.(ZII-ZI**2.GT.0.001)) THEN
CHI=(ZBI-ZB*ZI)/SQRT((ZBB-ZB**2)*(ZII-ZI**2))
IF(CHI.GT.CHI0) THEN
Z1=(ZI*ZBB-ZB*ZBI)/(ZBB-ZB**2)
IF(Z1.GT.0) THEN
PHI=MIN(Z0,Z1)
LAMBDA=(ZI-PHI)/ZB
ENDIF
ENDIF
ENDIF
ENDIF
RETURN
END