Skip to content

Commit

Permalink
further updates to EMHREBSDDIC program
Browse files Browse the repository at this point in the history
  • Loading branch information
@marcdegraef
marcdegraef committed Dec 19, 2024
1 parent b41ddeb commit b4ba7c1
Show file tree
Hide file tree
Showing 3 changed files with 163 additions and 69 deletions.
26 changes: 21 additions & 5 deletions Source/EMsoftOOLib/mod_DIC.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ recursive type(DIC_T) function DIC_constructor( nx, ny ) result(DIC)
integer(kind=irg), INTENT(IN) :: ny

integer(kind=irg) :: i, j
real(wp) :: arx, ary

! set pattern dimensions
DIC%nx = nx
Expand All @@ -161,14 +162,28 @@ recursive type(DIC_T) function DIC_constructor( nx, ny ) result(DIC)
allocate( DIC%x(0:nx-1), DIC%y(0:ny-1) )
DIC%rnxi = 1.0_wp/real(nx-1,wp)
DIC%rnyi = 1.0_wp/real(ny-1,wp)

! account for pattern aspect ratio
arx = 1.0_wp
ary = 1.0_wp
! if (nx.ne.ny) then
! if (nx.lt.ny) then
! arx = DIC%rnyi / DIC%rnxi
! else
! ary = DIC%rnxi / DIC%rnyi
! end if
! end if

do i=0,nx-1
DIC%x(i) = real(i,wp)
end do
DIC%x = DIC%x * DIC%rnxi
DIC%x = DIC%x * arx
! DIC%x = DIC%x * DIC%rnxi
do j=0,ny-1
DIC%y(j) = real(j,wp)
end do
DIC%y = DIC%y * DIC%rnyi
DIC%y = DIC%y * ary
! DIC%y = DIC%y * DIC%rnyi

end function DIC_constructor

Expand Down Expand Up @@ -458,8 +473,8 @@ recursive subroutine defineSR_(self, nbx, nby, PCx, PCy)

! define the coordinate arrays for the sub-region
allocate( self%xiX(0:self%nxSR-1), self%xiY(0:self%nySR-1) )
self%xiX = self%x(nbx:self%nx-nbx-1) - PCx
self%xiY = self%y(nby:self%ny-nby-1) - PCy
self%xiX = self%x(nbx:self%nx-nbx-1) - self%nxSR/2 ! PCx
self%xiY = self%y(nby:self%ny-nby-1) - self%nySR/2 ! PCy
if (self%verbose) call Message%printMessage(' defineSR_::xiX, xiY arrays allocated')

! allocate array for the product of the gradient and the Jacobian
Expand Down Expand Up @@ -530,11 +545,12 @@ recursive subroutine applyHomography_(self, h, PCx, PCy, dotarget)
!! apply a homography to the reference or target pattern and store in defpat

use mod_IO
use mod_math

IMPLICIT NONE

class(DIC_T),INTENT(INOUT) :: self
real(wp),INTENT(IN) :: h(0:7)
real(wp),INTENT(IN) :: h(8)
real(wp),INTENT(IN) :: PCx
real(wp),INTENT(IN) :: PCy
logical,INTENT(IN),OPTIONAL :: dotarget
Expand Down
24 changes: 15 additions & 9 deletions Source/EMsoftOOLib/program_mods/mod_HREBSDDIC.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -429,7 +429,7 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
numpats, TID
real(kind=dbl) :: interp_step, std
real(kind=dbl) :: Wnew(3,3), Winv(3,3), oldW(3,3)
real(wp) :: hg(8), W(3,3), ndp, oldnorm, SOL(8), Hessian(8,8), CIC
real(wp) :: hg(8), W(3,3), ndp, oldnorm, SOL(8), Hessian(8,8), CIC, PCx, PCy
character(11) :: dstr
character(15) :: tstrb
character(15) :: tstre
Expand Down Expand Up @@ -489,6 +489,8 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
L = binx * biny
recordsize = 4 * L
patsz = L
PCx = real(binx,wp)/2.0_wp
PCy = real(biny,wp)/2.0_wp

write (*,*) ' pattern dimensions : ', binx, biny, L, self%ppEBSDnml%ROI

Expand Down Expand Up @@ -547,7 +549,7 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
! define the border widths nbx and nby for the subregion
nbx = enl%nbx
nby = enl%nby
call DIC%defineSR( nbx, nby, 0.5_wp, 0.5_wp)
call DIC%defineSR( nbx, nby, PCx, PCy) ! 0.5_wp, 0.5_wp)

! generate the b-splines for the reference pattern and verify the accuracy
! also get the derivatives of the reference pattern to determine the Hessian
Expand All @@ -574,6 +576,8 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)

call DIC%setpattern( 'd', dble(targetpattern) )
call DIC%getbsplines(defp=.TRUE.)
! ! zero-mean and normalize the referenceSR and targetSR arrays
call DIC%applyZMN(dotarget=.TRUE.)

! initialize the homography to zero
hg = (/ (0.0_wp, i=1,8) /)
Expand All @@ -585,7 +589,8 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
do jj = 1, enl%maxnumit
! write (*,*) '-----------------'
! write (*,*) ' iteration # ',jj
call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.)
! call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.)
call DIC%applyHomography(hg, PCx, PCy, dotarget=.TRUE.)
call DIC%applyZMN(dotarget=.TRUE.)
call DIC%getresiduals(CIC)
call DIC%solveHessian(SOL, ndp)
Expand All @@ -595,13 +600,14 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
Winv = matrixInvert_wp( Wnew )
W = matmul( W, Winv )
W = W / W(3,3)
hg = DIC%getHomography(W)
if (ndp.lt.0.001D0) then
oldnorm = ndp
oldW = W
hg = DIC%getHomography(W)
else
W = oldW
ndp = oldnorm
! oldnorm = ndp
! oldW = W
! hg = DIC%getHomography(W)
! else
! W = oldW
! ndp = oldnorm
EXIT
end if
end do
Expand Down
182 changes: 127 additions & 55 deletions Source/TestPrograms/play.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ program EMplay
real(kind=sgl),allocatable :: tmppat(:,:)
real(kind=wp) :: DD, PCx, PCy, val, err, errmax, rnxi, rnyi, hg(8), W(3,3), gradCIC(8), Hessian(8,8), &
minx, miny, xi1max, xi2max, normdp, oldnorm, oldW(3,3), horiginal(8), CIC, sol(8), &
homographies(8,1000)
homographies(8,2500)
real(kind=dbl) :: Wnew(3,3), Winv(3,3), dx, dy, p2(3), Woriginal(3,3), alp, srt(3,3), srtrot(3,3)
integer(kind=irg) :: nx, ny, nxy, nbx, nby, i, ii, j, NSR, cnt, nxSR, nySR, jj, recordsize, ierr
real(wp) :: tol
Expand All @@ -85,7 +85,7 @@ program EMplay

! read the reference pattern from a binary file; the patterns have been
! pre-processed and scaled between 0.D0 and 1.D0
nx = 640
nx = 480
ny = 480

! instantiate the DIC class
Expand All @@ -98,10 +98,10 @@ program EMplay
recordsize = nxy * 4

allocate(tmppat(nx,ny), refpat(0:nx-1,0:ny-1), defpat(0:nx-1,0:ny-1))
open(20,file='/Users/mdg/playarea/DIC/verification-sq.data',status='old',form='unformatted')
! open(20,file='/Users/mdg/playarea/DIC/refpat.data',status='old',form='unformatted')
open(20,file='/Users/mdg/playarea/DIC/refpat.data',status='old',form='unformatted')
read(20) refpat
close(20,status='keep')
! read(20) refpat
! close(20,status='keep')
! refpat = dble(tmppat)
!
! open(20,file='/Users/mdg/Files/EMPlay/playarea/UCSB/GaN/pp27238.data',status='old',form='unformatted')
Expand All @@ -111,19 +111,21 @@ program EMplay

! read(15, rec=1, iostat=ierr) tmppat
! refpat = tmppat
if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
fname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographypatterns.data')
else
fname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographypatterns.data')
end if
open(unit=dataunit,file=trim(fname),&
status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr)
! if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
! fname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographypatterns2.data')
! else
! fname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographypatterns.data')
! end if
! open(unit=dataunit,file=trim(fname),&
! status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr)

! read(20) refpat
! read(dataunit, rec=1) tmppat
! refpat = dble(tmppat)
! read(20, rec=1) tmppat
read(20) tmppat
refpat = dble(tmppat)
! close(dataunit,status='keep')
! ! close(20,status='keep')
close(20,status='keep')

! write (*,*) ' tmppat entries : ', tmppat(1:2,1:2)
! write (*,*) ' reference pattern : ', minval(tmppat), maxval(tmppat)
Expand All @@ -150,52 +152,118 @@ program EMplay
! compute the Hessian matrix
call DIC%getHessian( Hessian )

! read 1000 random homographies from a file and run the fitting, then
write (*,*) 'Hessian'
do i=1,8
write (*,*) Hessian(i,1:8)
end do

! read 2500 random homographies from a file and run the fitting, then
! write result to another file for comparison with an IDL routine.
if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
gname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographies.data')
else
gname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographies1.data')
end if
open(unit=dataunit,file=trim(gname),status='unknown',form='unformatted')
read(dataunit) homographies
close(dataunit,status='keep')

if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
gname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographypatterns.data')
else
gname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographypatterns.data')
end if
open(unit=dataunit,file=trim(gname),&
status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr)

open(unit=28,file='output.txt',status='unknown',form='formatted')

do jj=1, 1000
call Message%printMessage(' ---------------------- ')
! if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
! gname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographies2.data')
! else
! gname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographies1.data')
! end if
! open(unit=dataunit,file=trim(gname),status='unknown',form='unformatted')
! read(dataunit) homographies
! close(dataunit,status='keep')

! if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
! gname = EMsoft%generateFilePath('EMdatapathname','UCSB/GaN/homographypatterns2.data')
! else
! gname = EMsoft%generateFilePath('EMdatapathname','playarea/UCSB/GaN/homographypatterns.data')
! end if
! open(unit=dataunit,file=trim(gname),&
! status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr)

open(20,file='/Users/mdg/playarea/DIC/check.data',status='old',form='unformatted')

! open(unit=28,file='output-verify.txt',status='unknown',form='formatted')

PCx = real(nx,wp)/2.0_wp
PCy = real(ny,wp)/2.0_wp

write (20) real(refpat)

horiginal = (/ 0.5D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.0D0, 0.3D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 25.0D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 0.0D0, 0.3D0, 0.D0, 0.D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 0.D0, 0.D0, 0.5D0, 0.D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 0.0D0, 0.D0, 0.D0, 20.0D0, 0.D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 0.0D0, 0.D0, 0.D0, 0.D0, 0.001D0, 0.D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.D0, 0.D0, 0.0D0, 0.D0, 0.D0, 0.D0, 0.0D0, 0.002D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)
call DIC%cleanup()
horiginal = (/ 0.01D0, 0.1D0, 50.0D0, 0.1D0, 0.5D0,-50.D0, 0.002D0, 0.002D0 /)
call DIC%applyHomography(horiginal, PCx, PCy)
defpat = DIC%getpattern('d', nx, ny)
write (20) real(defpat)

close(20,status='keep')

stop

do jj=1, 10 ! 2500
! call Message%printMessage(' ---------------------- ')
write (*,*) 'starting pattern ', jj
! here we deform the reference pattern to obtain a defpat array with known homography
! define the homography hg
horiginal = homographies(1:8,jj) ! (/ 0.002_wp, 0.001_wp, -0.004_wp, -0.001_wp, -0.003_wp, 0.005_wp, 0.0_wp, 0.0_wp /)
! ! horiginal = (/ (0.0_wp, i=1,8) /)
call DIC%applyHomography(horiginal, 0.5_wp, 0.5_wp)
! horiginal = homographies(1:8,jj) ! (/ 0.002_wp, 0.001_wp, -0.004_wp, -0.001_wp, -0.003_wp, 0.005_wp, 0.0_wp, 0.0_wp /)
horiginal = (/ (0.0_wp, i=1,8) /)
! call DIC%applyHomography(horiginal, 0.5_wp, 0.5_wp)

! read(dataunit,rec=jj, iostat=ierr) tmppat
! defpat = dble(tmppat)
! call DIC%setpattern('d', defpat)
! read(20,rec=jj, iostat=ierr) tmppat
read(20) tmppat
write (*,*) 'read error = ', ierr
defpat = dble(tmppat)
call DIC%setpattern('d', defpat)

! write (dataunit,rec=jj) DIC%getpattern('d',nx,ny)

write (dataunit,rec=jj) DIC%getpattern('d',nx,ny)
! if (1.eq.0) then

Woriginal = DIC%getShapeFunction(horiginal)

! get the derivatives of the reference pattern to determine the Hessian
call DIC%getbsplines(defp=.TRUE.)

! ! zero-mean and normalize the referenceSR and targetSR arrays
! call DIC%applyZMN(dotarget=.TRUE.)
call DIC%applyZMN(dotarget=.TRUE.)

! ! initial CIC function (to be minimized ...)
! call DIC%getresiduals( CIC )
call DIC%getresiduals( CIC )
! stop

! initialize the homography coefficients to zero
Expand All @@ -207,7 +275,7 @@ program EMplay
oldnorm = 100.0_wp

! and here we start the loop
do ii=1,500
do ii=1,50
write (*,*) ' iteration # ',ii
call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.)

Expand All @@ -224,16 +292,17 @@ program EMplay
Wnew = DIC%getShapeFunction(reshape(SOL, (/8/)))
Winv = matrixInvert_wp( Wnew )
W = matmul( W, Winv )
! W = W / W(3,3)
W = W / W(3,3)
hg = DIC%getHomography(W)
! write (*,*) hg
! write (*,*) ' norm(deltap) = ', normdp
! write (*,*) '------------'
if (normdp.lt.oldnorm) then
oldnorm = normdp
oldW = W
else
W = oldW
write (*,*) hg
write (*,*) ' norm(deltap) = ', normdp
write (*,*) '------------'
! if (normdp.lt.oldnorm) then
if (normdp.lt.0.001D0) then
! oldnorm = normdp
! oldW = W
! else
! W = oldW
EXIT
end if
! call Message%printMessage(' ---------------------- ')
Expand All @@ -258,9 +327,12 @@ program EMplay
! if (jj.eq.3) stop
call DIC%cleanup()

! end if

end do ! loop over jj

close(28,status='keep')
close(dataunit,status='keep')
close(20,status='keep')
! close(dataunit,status='keep')

end program EMplay

0 comments on commit b4ba7c1

Please sign in to comment.