Skip to content

Commit

Permalink
more updates to DIC programs
Browse files Browse the repository at this point in the history
  • Loading branch information
@marcdegraef
marcdegraef committed Jan 15, 2025
1 parent 40cdb3e commit ea57df3
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 162 deletions.
11 changes: 8 additions & 3 deletions Source/EMsoftOOLib/mod_DIC.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -609,9 +609,12 @@ recursive subroutine applyHomography_(self, h, PCx, PCy, dotarget)
write (*,*) ' applyHomography_::range(defpat) : ',minval(self%defpat), maxval(self%defpat)
end if


call self%sdef%initialize(self%x,self%y,self%defpat,self%kx,self%ky,iflag)


! finally get the b-spline coefficients for the deformed pattern
if (.not.tgt) then
call self%sdef%initialize(self%x,self%y,self%defpat,self%kx,self%ky,iflag)
if (self%verbose) call Message%printMessage(' applyHomography_::sdef class initialized')
call self%sdef%set_extrap_flag(.TRUE.)
end if
Expand Down Expand Up @@ -757,6 +760,7 @@ recursive subroutine solveHessian_(self, SOL, normdp)
xi2max = maxval( self%XiPrime(1,:) )
normdp = sqrt(xi1max * (SL(1,1)**2+SL(4,1)**2+SL(7,1)**2) + &
xi2max * (SL(2,1)**2+SL(5,1)**2+SL(8,1)**2) + SL(3,1)**2 + SL(6,1)**2)
! normdp = sum(SOL*SOL)
if (self%verbose) write (*,*) ' solveHessian_::normdp : ', normdp

end subroutine solveHessian_
Expand Down Expand Up @@ -788,7 +792,8 @@ recursive function getShapeFunction_(self, h, store) result(W)
W(1,2:3) = -h(2:3)
W(2,1) = -h(4)
W(2,3) = -h(6)
W(3,1:2) = -h(7:8)
W(3,1) = h(7)
W(3,2) = -h(8)

if (present(store)) then
if (store.eqv..TRUE.) then
Expand Down Expand Up @@ -819,7 +824,7 @@ recursive function getHomography_(self, W, store) result(h)
logical, INTENT(IN), OPTIONAL :: store
real(wp) :: h(8)

h = (/ 1.0_wp/W(1,1)-1.0_wp, -W(1,2), -W(1,3), -W(2,1), 1.0_wp/W(2,2)-1.0_wp, -W(2,3), -W(3,1), W(3,2) /)
h = (/ 1.0_wp/W(1,1)-1.0_wp, -W(1,2), -W(1,3), -W(2,1), 1.0_wp/W(2,2)-1.0_wp, -W(2,3), W(3,1), -W(3,2) /)

if (present(store)) then
if (store.eqv..TRUE.) then
Expand Down
38 changes: 22 additions & 16 deletions Source/EMsoftOOLib/program_mods/mod_HREBSDDIC.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -442,7 +442,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, PCx, PCy
real(wp) :: hg(8), W(3,3), ndp, oldnorm, SOL(8), Hessian(8,8), CIC, PCx, PCy, hpartial(8)
character(11) :: dstr
character(15) :: tstrb
character(15) :: tstre
Expand Down Expand Up @@ -502,8 +502,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
PCx = 0.5_wp ! real(binx,wp)/2.0_wp
PCy = 0.5_wp ! real(biny,wp)/2.0_wp

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

Expand All @@ -529,8 +529,11 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
read(dataunit,rec=offset) exptpattern
close(dataunit,status='keep')

exptpattern = exptpattern - minval(exptpattern)
exptpattern = exptpattern/maxval(exptpattern)

write (*,*) ' reference entries : ', exptpattern(1:2,1:2)
write (*,*) 'reference pattern : ', minval(exptpattern), maxval(exptpattern)
write (*,*) ' reference pattern : ', minval(exptpattern), maxval(exptpattern)

! allocate global arrays to store the output of the DIC routine.
call mem%alloc(normdp, (/ numpats /), 'normdp', initval=0.D0 )
Expand All @@ -548,7 +551,7 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)

! here we go parallel with OpenMP
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, DIC, nbx, nby, targetpattern, ii, jj, hg, W)&
!$OMP& PRIVATE(oldnorm, oldW, SOL, ndp, Wnew, Winv, io_int, i)
!$OMP& PRIVATE(oldnorm, oldW, SOL, ndp, Wnew, Winv, io_int, i, hpartial)

! set the thread ID
TID = OMP_GET_THREAD_NUM()
Expand All @@ -561,7 +564,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, PCx, PCy) ! 0.5_wp, 0.5_wp)
call DIC%defineSR( nbx, nby, PCx, PCy)

! 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 @@ -573,36 +576,40 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
! compute the Hessian matrix
call DIC%getHessian(Hessian)

hpartial = (/ (0.0_wp, i=1,8) /)
call DIC%applyHomography(hpartial, PCx, PCy)

! allocate the targetpattern array
call memth%alloc(targetpattern, (/ binx,biny /), 'targetpattern', TID=TID )

! in the main loop we iterate over all the experimental patterns in the input file
!$OMP DO SCHEDULE(DYNAMIC)
do ii=1, numpats
! write (*,*) '-----------------'
! write (*,*) 'starting pattern ', ii

!$OMP CRITICAL
read(dataunit,rec=ii) targetpattern
!$OMP END CRITICAL

targetpattern = targetpattern - minval(targetpattern)
targetpattern = targetpattern/maxval(targetpattern)

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.)
call DIC%getresiduals()

! initialize the homography to zero
hg = (/ (0.0_wp, i=1,8) /)
W = DIC%getShapeFunction(hg)
oldnorm = 100.0_wp
oldW = W

! loop until max iterations or convergence
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, PCx, PCy, dotarget=.TRUE.)
if (jj.eq.1) then
hpartial = (/ (0.0_wp, i=1,8) /)
W = DIC%getShapeFunction(hpartial)
end if
call DIC%applyHomography(hpartial, PCx, PCy, dotarget=.TRUE.)
call DIC%applyZMN(dotarget=.TRUE.)
call DIC%getresiduals(CIC)
call DIC%solveHessian(SOL, ndp)
Expand All @@ -612,7 +619,7 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
Winv = matrixInvert_wp( Wnew )
W = matmul( W, Winv )
W = W / W(3,3)
hg = DIC%getHomography(W)
hpartial = reshape(SOL, (/8/))
if (ndp.lt.enl%mindeltap) then
EXIT
end if
Expand All @@ -623,7 +630,6 @@ subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames)
io_int(2) = numpats
call Message%WriteValue(' completed # patterns/total ',io_int,2)
end if
W = matrixInvert_wp( W )
hg = DIC%getHomography(W)
homographies(1:8,ii) = dble(hg)
normdp(ii) = dble(ndp)
Expand Down
159 changes: 29 additions & 130 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,2500)
homographies(8,1600), hpartial(8)
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 Down Expand Up @@ -99,40 +99,17 @@ program EMplay

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/Files/EMPlay/playarea/DIC/refpat.data',status='old',form='unformatted')
if (trim(hostname).eq.'Mac-Studio.fios-router.home') then
fname = EMsoft%generateFilePath('EMdatapathname','DIC/refpat.data')
else
fname = EMsoft%generateFilePath('EMdatapathname','playarea/DIC/refpat.data')
end if
open(20,file=trim(fname),status='old',form='unformatted')
read(20) tmppat
close(20,status='keep')
refpat = dble(tmppat)
!
! open(20,file='/Users/mdg/Files/EMPlay/playarea/UCSB/GaN/pp27238.data',status='old',form='unformatted')
! recordsize = nxy*4
! open(unit=15,file='/Users/mdg/Files/EMPlay/playarea/UCSB/GaN/pp27238.data', &
! status='unknown',form='unformatted',access='direct',action='read',recl=recordsize,iostat=ierr)

! 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/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
! read(20, rec=1) tmppat
! read(20) tmppat
! refpat = dble(tmppat)
! close(dataunit,status='keep')
! close(20,status='keep')

! write (*,*) ' tmppat entries : ', tmppat(1:2,1:2)
! write (*,*) ' reference pattern : ', minval(tmppat), maxval(tmppat)

! refpat(20:40,20:22) = 1.0_wp
! refpat(20:22,20:40) = 1.0_wp

refpat = refpat - minval(refpat)
refpat = refpat/maxval(refpat)
call DIC%setpattern('r', refpat)

write (*,*) ' range refpat ', minval(refpat), maxval(refpat)
Expand Down Expand Up @@ -162,112 +139,29 @@ program EMplay
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/homographies2.data')
! else
! gname = EMsoft%generateFilePath('EMdatapathname','playarea/DIC/test/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')
gname = EMsoft%generateFilePath('EMdatapathname','DIC/test/homographypatterns2.data')
else
gname = EMsoft%generateFilePath('EMdatapathname','playarea/DIC/test/homographypatterns.data')
gname = EMsoft%generateFilePath('EMdatapathname','playarea/DIC/test/homographypatterns2.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.2D0, 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.5D0, 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, 50.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.5D0, 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.2D0, 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, 50.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.001D0 /)
! 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')

! ! check

! horiginal = (/ 0.01D0, 0.1D0, 50.0D0, 0.1D0, 0.5D0,-50.D0, 0.002D0, 0.002D0 /)
! write(*,*) horiginal
! W = DIC%getShapeFunction(horiginal)
! hg = DIC%getHomography(W)
! write(*,*) hg

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

! stop
horiginal = (/ (0.0_wp, i=1,8) /)
call DIC%applyHomography(horiginal, PCx, PCy)

do jj=1, 2500
do jj=1, 1600
! call Message%printMessage(' ---------------------- ')
if (mod(jj,100).eq.0) 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, PCx, PCy)

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

! 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
Expand All @@ -283,15 +177,18 @@ program EMplay
! initialize the homography coefficients to zero
! in a more general implementation we might initialize them to the
! values for one of the nearest neighbor patterns
hg = (/ (0.0_wp, i=1,8) /)
W = DIC%getShapeFunction(hg)

oldnorm = 100.0_wp

! and here we start the loop
do ii=1,50
do ii=1,20 ! 50
! write (*,*) ' iteration # ',ii
call DIC%applyHomography(hg, PCx, PCy, dotarget=.TRUE.)
if (ii.eq.1) then ! initialize to identity homography in first cycle
hpartial = (/ (0.0_wp, i=1,8) /)
W = DIC%getShapeFunction(hpartial)
end if

call DIC%applyHomography(hpartial, PCx, PCy, dotarget=.TRUE.)

! zero-mean and normalize the deformed sub-region array wtarget
call DIC%applyZMN( dotarget = .TRUE. )
Expand All @@ -308,11 +205,13 @@ program EMplay
W = matmul( W, Winv )
W = W / W(3,3)
hg = DIC%getHomography(W)
hpartial = reshape(SOL, (/8/))
! write (*,*) '------------'
! write (*,*) hg
! write (*,*) ' norm(deltap) = ', normdp
! write (*,*) '------------'
! if (normdp.lt.oldnorm) then
if (normdp.lt.0.001D0) then
if (normdp.lt.0.0001D0) then
! oldnorm = normdp
! oldW = W
! else
Expand All @@ -322,8 +221,8 @@ program EMplay
! call Message%printMessage(' ---------------------- ')
end do

W = matrixInvert_wp( W )
W = W / W(3,3)
! W = matrixInvert_wp( W )
! W = W / W(3,3)
hg = DIC%getHomography(W)
! write (*,*) ''
! write (*,*) ' Final homography : '
Expand Down
Loading

0 comments on commit ea57df3

Please sign in to comment.