From 017a31d3767804c0414c5eaa4e52df13c5ae239c Mon Sep 17 00:00:00 2001 From: Marc DeGraef Date: Tue, 3 Dec 2024 09:25:21 -0500 Subject: [PATCH] functional implementation of HR-EBSD-DIC algorithm code compiles and runs correctly on a test case; needs more integration with the rest of the package --- Source/EMsoftOOLib/CMakeLists.txt | 4 + Source/EMsoftOOLib/mod_DIC.f90 | 776 ++++++++++++++++++++++++++++++ Source/EMsoftOOLib/mod_math.f90 | 43 ++ Source/TestPrograms/play.f90 | 533 ++------------------ 4 files changed, 872 insertions(+), 484 deletions(-) create mode 100644 Source/EMsoftOOLib/mod_DIC.f90 diff --git a/Source/EMsoftOOLib/CMakeLists.txt b/Source/EMsoftOOLib/CMakeLists.txt index d5f4cec..546c3a8 100644 --- a/Source/EMsoftOOLib/CMakeLists.txt +++ b/Source/EMsoftOOLib/CMakeLists.txt @@ -108,6 +108,8 @@ set(EMsoftOOLib_SRCS ${EMsoftOOLib_SOURCE_DIR}/mod_Detector.f90 ${EMsoftOOLib_SOURCE_DIR}/mod_SphereIndexer.f90 ${EMsoftOOLib_SOURCE_DIR}/mod_SphInxSupport.f90 + ${EMsoftOOLib_SOURCE_DIR}/mod_DIC.f90 + # ${EMsoftOOLib_SOURCE_DIR}/mod_shapes.f90 ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_sampleRFZ.f90 @@ -266,6 +268,7 @@ target_link_libraries(EMsoftOOLib # EMsoftOOLib_Cpp ${EMSOFTOO_hdf5LinkLibs} jsonfortran + ${BSPLINEFORTRAN_DIR}/libbspline-fortran.a bcls::bcls ) @@ -290,6 +293,7 @@ target_include_directories(EMsoftOOLib ${FFTW3_INCLUDE_DIR} ${HDF5_INCLUDE_DIR} ${HDF5_INCLUDE_DIR_FORTRAN} + ${BSPLINEFORTRAN_DIR} PRIVATE "${EMsoftOO_SOURCE_DIR}/Source" ) diff --git a/Source/EMsoftOOLib/mod_DIC.f90 b/Source/EMsoftOOLib/mod_DIC.f90 new file mode 100644 index 0000000..59df300 --- /dev/null +++ b/Source/EMsoftOOLib/mod_DIC.f90 @@ -0,0 +1,776 @@ +! ################################################################### +! Copyright (c) 2014-2024, Marc De Graef Research Group/Carnegie Mellon University +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without modification, are +! permitted provided that the following conditions are met: +! +! - Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! - Redistributions in binary form must reproduce the above copyright notice, this +! list of conditions and the following disclaimer in the documentation and/or +! other materials provided with the distribution. +! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names +! of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE +! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! ################################################################### + +module mod_DIC + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! this class implements the DIC algorithm as described in the series of + !! papers by Ernould et al. + !! + +use mod_kinds +use mod_global +use bspline_module +use bspline_kinds_module, only: wp, ip + +IMPLICIT NONE + +type, public :: DIC_T + private + real(wp) :: W(3,3) + real(wp) :: h(8) + real(wp) :: CIC + integer(ip) :: kx = 5 ! spline order + integer(ip) :: ky = 5 ! spline order + integer(kind=irg) :: nx ! pattern x-size + integer(kind=irg) :: ny ! pattern y-size + real(wp) :: rnxi ! pattern x scale factor + real(wp) :: rnyi ! pattern y scale factor + integer(kind=irg) :: nxSR ! sub-region x-size + integer(kind=irg) :: nySR ! sub-region y-size + integer(kind=irg) :: nSR ! sub-region total size + integer(kind=irg) :: nbx ! sub-region border x-width + integer(kind=irg) :: nby ! sub-region border y-width + type(bspline_2d) :: sref + type(bspline_2d) :: sdef + real(wp) :: tol = 100.0_wp * epsilon(1.0_wp) + logical :: verbose = .FALSE. ! useful for debugging + + ! the arrays that are tied to a given reference/target pattern are defined here + real(wp),allocatable :: x(:) + real(wp),allocatable :: y(:) + real(wp),allocatable :: xiX(:) + real(wp),allocatable :: xiY(:) + real(wp),allocatable :: XiPrime(:,:) + real(wp),allocatable :: refpat(:,:) + real(wp),allocatable :: defpat(:,:) + real(wp),allocatable :: wtarget(:) + real(wp),allocatable :: gradx(:,:) + real(wp),allocatable :: grady(:,:) + real(wp),allocatable :: gxSR(:) + real(wp),allocatable :: gySR(:) + real(wp),allocatable :: GradJac(:,:) + real(wp),allocatable :: referenceSR(:) + real(wp),allocatable :: targetSR(:) + real(wp),allocatable :: refzmn(:) + real(wp),allocatable :: tarzmn(:) + real(wp),allocatable :: residuals(:) + real(wp) :: Hessian(8,8) + real(wp) :: gradCIC(8) + real(wp) :: refmean + real(wp) :: tarmean + real(wp) :: refstdev + real(wp) :: tarstdev + + contains + private + + procedure, pass(self) :: getShapeFunction_ + procedure, pass(self) :: getHomography_ + procedure, pass(self) :: homography2Fe_ + procedure, pass(self) :: Fe2homography_ + procedure, pass(self) :: getbsplines_ + procedure, pass(self) :: getHessian_ + procedure, pass(self) :: solveHessian_ + procedure, pass(self) :: defineSR_ + procedure, pass(self) :: applyZMN_ + procedure, pass(self) :: applyHomography_ + procedure, pass(self) :: getresiduals_ + procedure, pass(self) :: setverbose_ + procedure, pass(self) :: setpattern_ + + final :: DIC_destructor + + generic, public :: getShapeFunction => getShapeFunction_ + generic, public :: getHomography => getHomography_ + generic, public :: Fe2homography=> Fe2homography_ + generic, public :: homography2Fe=> homography2Fe_ + generic, public :: getbsplines => getbsplines_ + generic, public :: getHessian => getHessian_ + generic, public :: solveHessian => solveHessian_ + generic, public :: defineSR => defineSR_ + generic, public :: applyZMN => applyZMN_ + generic, public :: applyHomography => applyHomography_ + generic, public :: getresiduals => getresiduals_ + generic, public :: setverbose => setverbose_ + generic, public :: setpattern => setpattern_ + +end type DIC_T + +! the constructor routine for this class +interface DIC_T + module procedure DIC_constructor +end interface DIC_T + +contains + +!-------------------------------------------------------------------------- +type(DIC_T) function DIC_constructor( nx, ny ) result(DIC) +!DEC$ ATTRIBUTES DLLEXPORT :: DIC_constructor + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! constructor for the DIC_T Class; + +use bspline_kinds_module, only: wp, ip + +IMPLICIT NONE + +integer(kind=irg), INTENT(IN) :: nx +integer(kind=irg), INTENT(IN) :: ny + +integer(kind=irg) :: i, j + +! set pattern dimensions +DIC%nx = nx +DIC%ny = ny + +! allocate and initialize the normalized coordinate arrays +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) +do i=0,nx-1 + DIC%x(i) = real(i,wp) +end do +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 + +end function DIC_constructor + +!-------------------------------------------------------------------------- +subroutine DIC_destructor( self ) +!DEC$ ATTRIBUTES DLLEXPORT :: DIC_destructor + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! deallocate all arrays + +IMPLICIT NONE + +type(DIC_T), INTENT(INOUT) :: self + +call self%sref%destroy() +call self%sdef%destroy() + +call reportDestructor('DIC_T') + +end subroutine DIC_destructor + +!-------------------------------------------------------------------------- +recursive subroutine setverbose_(self, v) +!DEC$ ATTRIBUTES DLLEXPORT :: setverbose_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! set the verbose parameter + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +logical, INTENT(IN) :: v + +type(IO_T) :: Message + +self%verbose = v + +if (v.eqv..TRUE.) then + call Message%printMessage(' DIC_T class verbosity turned on') +else + call Message%printMessage(' DIC_T class verbosity turned off') +end if + +end subroutine setverbose_ + +!-------------------------------------------------------------------------- +recursive subroutine setpattern_(self, rp, pattern) +!DEC$ ATTRIBUTES DLLEXPORT :: setpattern_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! set the reference or target pattern + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +character(1),INTENT(IN) :: rp +real(wp),INTENT(IN) :: pattern(:,:) + +integer(kind=irg) :: sz(2) + +sz = shape(pattern) + +if (rp.eq.'r') then + allocate( self%refpat(0:sz(1)-1,0:sz(2)-1) ) + self%refpat = pattern +end if + +if (rp.eq.'p') then + allocate( self%defpat(0:sz(1)-1,0:sz(2)-1) ) + self%defpat = pattern +end if + +end subroutine setpattern_ + +!-------------------------------------------------------------------------- +recursive subroutine getbsplines_(self, verify, refp, defp, grads) +!DEC$ ATTRIBUTES DLLEXPORT :: getbsplines_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! compute the bspline coefficients for the reference pattern and, optionally, + !! also for the deformed pattern; optionally also compute the gradients of the + !! reference pattern + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +logical, INTENT(IN), OPTIONAL :: verify ! reconstruct the reference pattern from splines ? +logical, INTENT(IN), OPTIONAL :: refp ! get bsplines for deformed pattern ? +logical, INTENT(IN), OPTIONAL :: defp ! get bsplines for deformed pattern ? +logical, INTENT(IN), OPTIONAL :: grads ! compute gradients of the reference pattern ? + +type(IO_T) :: Message +integer(ip) :: iflag +integer(kind=irg) :: i, j +real(wp) :: val, err, errmax +real(kind=dbl) :: io_real(1) + +if (present(refp)) then + if (refp.eqv..TRUE.) then + ! initialize the bsplines for the reference pattern + call self%sref%initialize(self%x,self%y,self%refpat,self%kx,self%ky,iflag) + ! for potential later calls, allow for extrapolation + call self%sref%set_extrap_flag(.TRUE.) + if (self%verbose) call Message%printMessage(' b-splines for reference pattern completed') + end if +end if + +if (present(defp)) then + if (defp.eqv..TRUE.) then +! initialize the bsplines for the deformed pattern + call self%sdef%initialize(self%x,self%y,self%defpat,self%kx,self%ky,iflag) +! for potential later calls, allow for extrapolation + call self%sdef%set_extrap_flag(.TRUE.) + if (self%verbose) call Message%printMessage(' b-splines for target pattern completed') + end if +end if + +if (present(verify)) then + if (verify.eqv..TRUE.) then + do i=0,self%nx-1 + do j=0,self%ny-1 + ! determine how well the spline coefficients represent the original pattern + call self%sref%evaluate(self%x(i),self%y(j),0,0,val,iflag) + err = abs(self%refpat(i,j)-val) + errmax = max(err,errmax) + end do + end do + ! check max error against tolerance + io_real(1) = dble(errmax) + call Message%WriteValue(' max b-spline reconstruction error of reference pattern:', io_real, 1) + if (errmax >= self%tol) then + call Message%printMessage(' ** reference reconstruction test failed ** ') + else + call Message%printMessage(' ** reference reconstruction test passed ** ') + end if + end if +end if + +if (present(grads)) then + if (grads.eqv..TRUE.) then + allocate( self%gradx(0:self%nx-1,0:self%ny-1), self%grady(0:self%nx-1,0:self%ny-1) ) + do i=0,self%nx-1 + do j=0,self%ny-1 + ! extract the gradient arrays + call self%sref%evaluate(self%x(i),self%y(j),1,0,self%gradx(i,j),iflag) + call self%sref%evaluate(self%x(i),self%y(j),0,1,self%grady(i,j),iflag) + end do + end do + if (self%verbose) call Message%printMessage(' gradient of reference pattern completed') + self%gxSR = reshape( self%gradx(self%nbx:self%nx-self%nbx-1,self%nby:self%ny-self%nby-1), (/ self%NSR /) ) + self%gySR = reshape( self%grady(self%nbx:self%nx-self%nbx-1,self%nby:self%ny-self%nby-1), (/ self%NSR /) ) + end if +end if + +end subroutine getbsplines_ + +!-------------------------------------------------------------------------- +recursive subroutine defineSR_(self, nbx, nby, PCx, PCy) +!DEC$ ATTRIBUTES DLLEXPORT :: defineSR_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! set dimensions and allocate all SR-related arrays + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +integer(kind=irg),INTENT(IN) :: nbx +integer(kind=irg),INTENT(IN) :: nby +real(wp),INTENT(IN) :: PCx +real(wp),INTENT(IN) :: PCy + +type(IO_T) :: Message + +! set the array dimensions for the sub-region +self%nbx = nbx +self%nby = nby +self%nxSR = self%nx-2*nbx +self%nySR = self%ny-2*nby +self%NSR = self%nxSR * self%nySR + +! allocate and populate relevant arrays +allocate( self%referenceSR(0:self%NSR-1), self%gxSR(0:self%NSR-1), self%gySR(0:self%NSR-1), & + self%refzmn(0:self%NSR-1), self%tarzmn(0:self%NSR-1) ) +self%referenceSR = reshape( self%refpat(nbx:self%nx-nbx-1,nby:self%ny-nby-1), (/ self%NSR /) ) + +! 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 + +! allocate array for the product of the gradient and the Jacobian +allocate( self%GradJac(0:7,0:self%NSR-1) ) + +! allocate the residuals array +allocate( self%residuals(0:self%NSR-1) ) + +! and finally the targetdef array +allocate( self%targetSR(0:self%NSR-1) ) + +if (self%verbose) call Message%printMessage(' sub-region arrays allocated') + +end subroutine defineSR_ + +!-------------------------------------------------------------------------- +recursive subroutine applyZMN_(self, doreference, dotarget) +!DEC$ ATTRIBUTES DLLEXPORT :: applyZMN_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! apply zero-mean and normalization to a pattern + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +logical,INTENT(IN), OPTIONAL :: doreference +logical,INTENT(IN), OPTIONAL :: dotarget +real(wp) :: io_real(2) + +type(IO_T) :: Message + +if (present(doreference)) then + if (doreference.eqv..TRUE.) then + self%refmean = sum(self%referenceSR)/dble(self%NSR) + self%refstdev = sqrt( sum( (self%referenceSR - self%refmean)**2 ) ) + self%refzmn = (self%referenceSR-self%refmean)/self%refstdev + io_real(1:2) = (/ self%refmean, self%refstdev /) + if (self%verbose) call Message%WriteValue(' reference mean/stdev : ', io_real, 2) + end if +end if + +if (present(dotarget)) then + if (dotarget.eqv..TRUE.) then + self%tarmean = sum(self%targetSR)/dble(self%NSR) + self%tarstdev = sqrt( sum( (self%targetSR - self%tarmean)**2 ) ) + self%tarzmn = (self%targetSR-self%tarmean)/self%tarstdev + io_real(1:2) = (/ self%tarmean, self%tarstdev /) + if (self%verbose) call Message%WriteValue(' target mean/stdev : ', io_real, 2) + end if +end if + +end subroutine applyZMN_ + +!-------------------------------------------------------------------------- +recursive subroutine applyHomography_(self, h, PCx, PCy, dotarget) +!DEC$ ATTRIBUTES DLLEXPORT :: applyHomography_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! apply a homography to the reference or target pattern and store in defpat + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp),INTENT(IN) :: h(0:7) +real(wp),INTENT(IN) :: PCx +real(wp),INTENT(IN) :: PCy +logical,INTENT(IN),OPTIONAL :: dotarget + +type(IO_T) :: Message + +real(wp) :: W(3,3), p(3) +integer(ip) :: iflag +integer(kind=irg) :: cnt, i, j, lnx, lny +logical :: tgt + +tgt = .FALSE. +if (present(dotarget)) then + if (dotarget.eqv..TRUE.) then + tgt = .TRUE. + if (allocated(self%XiPrime)) deallocate(self%XiPrime) + end if +end if + +lnx = self%nx +lny = self%ny + +! convert to shape function +W = self%getShapeFunction(h) + +! determine the XiPrime coordinates so that we can apply the deformation by +! means of the evaluate method in the bspline class +if (.not.allocated(self%XiPrime)) allocate( self%XiPrime(0:1, 0:lnx*lny-1)) +cnt = 0 +do j = 0, lny-1 + do i = 0, lnx-1 + p = matmul( W, (/ self%x(i)-PCx, self%y(j)-PCy, 1.0_wp /) ) + if (p(3).ne.0.0_wp) then + self%XiPrime(0:1,cnt) = p(1:2)/p(3) + else + self%XiPrime(0:1,cnt) = p(1:2) + end if + cnt = cnt + 1 + end do +end do + +! and interpolate/extrapolate the deformed pattern as needed +if (.not.allocated(self%defpat)) allocate( self%defpat( 0:lnx-1, 0:lny-1 )) +cnt = 0 +do j=0,lny-1 + do i=0,lnx-1 + iflag = 0_wp + if (tgt) then + call self%sdef%evaluate(self%XiPrime(0,cnt)+PCx,self%XiPrime(1,cnt)+PCy,0,0,self%defpat(i,j),iflag) + else + call self%sref%evaluate(self%XiPrime(0,cnt)+PCx,self%XiPrime(1,cnt)+PCy,0,0,self%defpat(i,j),iflag) + end if + cnt = cnt+1 + end do +end do +! reduce to interval [0,1] +self%defpat = self%defpat - minval(self%defpat) +self%defpat = self%defpat/maxval(self%defpat) + +! 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) + call self%sdef%set_extrap_flag(.TRUE.) +end if + +! copy the sub-region of the deformed target pattern into targetSR array +if (.not.allocated(self%targetSR)) then + allocate( self%targetSR( 0:self%NSR-1 ) ) +end if +self%targetSR = reshape( self%defpat( self%nbx:self%nx-self%nbx-1, self%nby:self%ny-self%nby-1 ), (/ self%NSR /) ) + +end subroutine applyHomography_ + +!-------------------------------------------------------------------------- +recursive subroutine getresiduals_(self, CIC) +!DEC$ ATTRIBUTES DLLEXPORT :: getresiduals_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! compute the residuals and optionally CIC + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp),INTENT(OUT),OPTIONAL :: CIC + +type(IO_T) :: Message +real(kind=dbl) :: io_real(1) + +self%residuals = self%refzmn-self%tarzmn + +if (present(CIC)) then + self%CIC = sum( self%residuals**2 ) + CIC = self%CIC + if (self%verbose) then + io_real(1) = dble(CIC) + call Message%WriteValue(' CIC value : ', io_real, 1) + end if +end if + +end subroutine getresiduals_ + +!-------------------------------------------------------------------------- +recursive subroutine getHessian_(self, Hessian) +!DEC$ ATTRIBUTES DLLEXPORT :: getHessian_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! compute the Hessian matrix for the reference pattern + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp),INTENT(OUT),OPTIONAL :: Hessian(8,8) + +type(IO_T) :: Message +integer(kind=irg) :: cnt, i, j +real(wp) :: GJ(8) + +cnt = 0 +self%Hessian = 0.0_wp +do j = 0, self%nySR-1 + do i = 0, self%nxSR-1 + GJ = (/ self%gxSR(cnt)*self%xiX(i), self%gxSR(cnt)*self%xiY(j), self%gxSR(cnt), & + self%gySR(cnt)*self%xiX(i), self%gySR(cnt)*self%xiY(j), self%gySR(cnt), & + self%gxSR(cnt)*self%xiX(i)**2-self%gySR(cnt)*self%xiX(i)*self%xiY(j), & + self%gySR(cnt)*self%xiY(j)**2-self%gxSR(cnt)*self%xiX(i)*self%xiY(j) /) + self%GradJac(0:7,cnt) = GJ + self%Hessian = self%Hessian + matmul( reshape(GJ, (/8,1/)), reshape(GJ, (/1,8/)) ) + cnt = cnt + 1 + end do +end do +self%Hessian = self%Hessian * 2.0_wp / self%refstdev**2 + +if (present(Hessian)) Hessian = self%Hessian + +if (self%verbose) call Message%printMessage(' reference Hessian computed') + +end subroutine getHessian_ + +!-------------------------------------------------------------------------- +recursive subroutine solveHessian_(self, SOL, normdp) +!DEC$ ATTRIBUTES DLLEXPORT :: solveHessian_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! solve the Hessian equation using Cholesky decomposition (LAPACK) + +use mod_IO + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp),INTENT(INOUT) :: SOL(8) +real(wp),INTENT(INOUT) :: normdp + +real(wp) :: Hessiancopy(8,8) ! local copy +real(wp) :: xi1max, xi2max + +! LAPACK variables +character :: UPLO ='U' +integer :: NN = 8 +integer :: NRHS = 1 +integer :: LDA = 8 +real(wp) :: SL(1:8,1) +integer :: LDB = 8 +integer :: INFO + +! compute gradCIC (right-hand side of the equation) +self%gradCIC = matmul(self%GradJac, self%residuals) +self%gradCIC = self%gradCIC * 2.0_wp/self%refstdev + +! solve by Cholesky decomposition +Hessiancopy = self%Hessian +SL(1:8,1) = -self%gradCIC(1:8) +call DPOSV(UPLO, NN, NRHS, Hessiancopy, LDA, SL, LDB, INFO) +SOL(1:8) = SL(1:8,1) + +! get the weighted norm of the solution vector +xi1max = maxval( self%XiPrime(0,:) ) +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) + +end subroutine solveHessian_ + +!-------------------------------------------------------------------------- +recursive function getShapeFunction_(self, h, store) result(W) +!DEC$ ATTRIBUTES DLLEXPORT :: getShapeFunction_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! convert input homography to a 3x3 shape function W and optionally store it + !! as a class variable + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp), INTENT(IN) :: h(8) +logical, INTENT(IN), OPTIONAL :: store +real(wp) :: W(3,3) + +W(1,1) = 1.0_wp + h(1) +W(2,2) = 1.0_wp + h(5) +W(3,3) = 1.0_wp + +W(1,2:3) = h(2:3) +W(2,1) = h(4) +W(2,3) = h(6) +W(3,1:2) = h(7:8) + +if (present(store)) then + if (store.eqv..TRUE.) then + self%W = W + end if +end if + +end function getShapeFunction_ + +!-------------------------------------------------------------------------- +recursive function getHomography_(self, W, store) result(h) +!DEC$ ATTRIBUTES DLLEXPORT :: getHomography_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! convert input 3x3 shape function to homography h and optionally store it + !! as a class variable + +IMPLICIT NONE + +class(DIC_T),INTENT(INOUT) :: self +real(wp), INTENT(IN) :: W(3,3) +logical, INTENT(IN), OPTIONAL :: store +real(wp) :: h(8) + +h = (/ W(1,1)-1.0_wp, W(1,2), W(1,3), W(2,1), W(2,2)-1.0_wp, W(2,3), W(3,1), W(3,2) /) + +if (present(store)) then + if (store.eqv..TRUE.) then + self%h = h + end if +end if + +end function getHomography_ + +!-------------------------------------------------------------------------- +recursive function homography2Fe_(self, h, PC, DD) result(W) +!DEC$ ATTRIBUTES DLLEXPORT :: homography2Fe_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! convert a homography to a reduced deformation tensor (Fe-hat in Ernould papers) + +class(DIC_T),INTENT(INOUT) :: self +real(wp), INTENT(IN) :: h(0:7) +real(wp), INTENT(IN) :: PC(0:1) +real(wp), INTENT(IN) :: DD +real(wp) :: W(3,3) + +real(wp) :: beta0, Fe11, Fe12, Fe13, Fe21, Fe22, Fe23, Fe31, Fe32, Fe33 + +! h11, h12, h13, h21, h22, h23, h31, h32 +! 0 1 2 3 4 5 6 7 + +beta0 = 1.0_wp - h(6) * PC(0) - h(7) * PC(1) +Fe11 = 1.0_wp + h(0) + h(6) * PC(0) +Fe12 = h(1) + h(7) * PC(0) +Fe13 = (h(2) - h(0)*PC(0) - h(1)*PC(1) + PC(0)*(beta0 - 1.0_wp))/DD +Fe21 = h(3) + h(6) * PC(1) +Fe22 = 1.0_wp + h(4) + h(7) * PC(1) +Fe23 = (h(5) - h(3)*PC(0) - h(4)*PC(1) + PC(1)*(beta0 - 1.0_wp))/DD +Fe31 = DD * h(6) +Fe32 = DD * h(7) +Fe33 = beta0 + +W = reshape( (/ Fe11, Fe12, Fe13, Fe21, Fe22, Fe23, Fe31, Fe32, Fe33 /), (/ 3,3 /)) / beta0 + +end function homography2Fe_ + +!-------------------------------------------------------------------------- +recursive function Fe2homography_(self, Fe_input, PC, DD) result(h) +!DEC$ ATTRIBUTES DLLEXPORT :: Fe2homography_ + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! convert a reduced deformation tensor (with (3,3) entry = 1.0_wp) to a homography + +class(DIC_T),INTENT(INOUT) :: self +real(wp), INTENT(IN) :: Fe_input(0:2,0:2) +real(wp), INTENT(IN) :: PC(0:1) +real(wp), INTENT(IN) :: DD +real(wp) :: h(8) + +real(wp) :: Fe(0:2,0:2), g0, g11, g22, g13, g23, h11, h12, h13, h21, h22, h23, h31, h32 + +Fe = Fe_input/Fe_input(2,2) + +g0 = DD + Fe(2,0) * PC(0) + Fe(2,1) * PC(1) +g11 = DD * Fe(0,0) - Fe(2,0) * PC(0) - g0 +g22 = DD * Fe(1,1) - Fe(2,1) * PC(1) - g0 +g13 = DD * ((Fe(0,0) - 1.D0) * PC(0) + Fe(0,1) * PC(1) + Fe(0,2) * DD) + PC(0) * (DD - g0) +g23 = DD * (Fe(1,0) * PC(0) + (Fe(1,1) - 1.D0) * PC(1) + Fe(1,2) * DD) + PC(1) * (DD - g0) +h11 = g11 +h12 = DD * Fe(0,1) - Fe(2,1) * PC(0) +h13 = g13 +h21 = DD * Fe(1,0) - Fe(2,0) * PC(1) +h22 = g22 +h23 = g23 +h31 = Fe(2,0) +h32 = Fe(2,1) + +h = (/ h11, h12, h13, h21, h22, h23, h31, h32 /) / g0 + +end function Fe2homography_ + + + + + +end module mod_DIC diff --git a/Source/EMsoftOOLib/mod_math.f90 b/Source/EMsoftOOLib/mod_math.f90 index 8a6997e..21784f9 100644 --- a/Source/EMsoftOOLib/mod_math.f90 +++ b/Source/EMsoftOOLib/mod_math.f90 @@ -479,6 +479,49 @@ recursive subroutine mInvert(a,b,uni) end subroutine mInvert +!-------------------------------------------------------------------------- +recursive function matrixInvert_wp(a) result(b) +!DEC$ ATTRIBUTES DLLEXPORT :: matrixInvert_wp + !! author: MDG + !! version: 1.0 + !! date: 12/01/24 + !! + !! Invert a double precision 3x3 matrix in the "wp" kind used by bspline-fortran module; + !! this is only used in code that uses the bsplines. + +use bspline_kinds_module, only: wp +use mod_io + +real(wp), INTENT(IN) :: a(3,3) +real(wp) :: b(3,3) + +type(IO_T) :: Message +real(wp) :: d + +d = a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+ & + a(1,3)*a(2,1)*a(3,2)-a(1,3)*a(2,2)*a(3,1)- & + a(1,2)*a(2,1)*a(3,3)-a(1,1)*a(2,3)*a(3,2) +if (d.ne.0.D0) then + b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2) + b(1,2)=a(1,3)*a(3,2)-a(1,2)*a(3,3) + b(1,3)=a(1,2)*a(2,3)-a(1,3)*a(2,2) + b(2,1)=a(2,3)*a(3,1)-a(2,1)*a(3,3) + b(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1) + b(2,3)=a(1,3)*a(2,1)-a(1,1)*a(2,3) + b(3,1)=a(2,1)*a(3,2)-a(2,2)*a(3,1) + b(3,2)=a(1,2)*a(3,1)-a(1,1)*a(3,2) + b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1) + b = b/d +else + do i=1,3 + write (*,*) (a(i,j),j=1,3) + end do + call Message%printMessage('matrixInvert_wp: matrix has zero determinant') + b = a +end if + +end function matrixInvert_wp + !-------------------------------------------------------------------------- ! ! SUBROUTINE:cInvert diff --git a/Source/TestPrograms/play.f90 b/Source/TestPrograms/play.f90 index 952def3..17e6fc7 100644 --- a/Source/TestPrograms/play.f90 +++ b/Source/TestPrograms/play.f90 @@ -44,6 +44,7 @@ program EMplay use mod_io use mod_symmetry use mod_crystallography +use mod_DIC use HDF5 use mod_HDFsupport use bspline_module @@ -51,141 +52,25 @@ program EMplay IMPLICIT NONE -interface - function getShapeFunction(h) result(W) - - use bspline_kinds_module, only: wp - - real(wp), INTENT(IN) :: h(8) - real(wp) :: W(3,3) - - end function getShapeFunction - - function getHomography(W) result(h) - - use bspline_kinds_module, only: wp - - real(wp), INTENT(IN) :: W(3,3) - real(wp) :: h(8) - - end function getHomography - - function matrixInvert(a) result(b) - - use bspline_kinds_module, only: wp - use mod_io - - real(wp), INTENT(IN) :: a(3,3) - real(wp) :: b(3,3) - - end function matrixInvert - - function Fe2homography(Fe_input, PC, DD) result(h) - - use bspline_kinds_module, only: wp - - real(wp), INTENT(IN) :: Fe_input(0:2,0:2) - real(wp), INTENT(IN) :: PC(0:1) - real(wp), INTENT(IN) :: DD - real(wp) :: h(8) - - end function Fe2homography - - function homography2Fe(p, PC, DD) result(W) - - use bspline_kinds_module, only: wp - - real(wp), INTENT(IN) :: p(0:7) - real(wp), INTENT(IN) :: PC(0:1) - real(wp), INTENT(IN) :: DD - real(wp) :: W(3,3) - - end function homography2Fe -end interface - - - character(fnlen) :: progname = 'EMplay.f90' character(fnlen) :: progdesc = 'test program for HR-EBSD-DIC' type(EMsoft_T) :: EMsoft -type(bspline_2d) :: sref, sdef +type(DIC_T) :: DIC type(q_T) :: qu type(o_T) :: om type(e_T) :: eu -real(kind=wp), allocatable :: refpat(:,:), defpat(:,:), x(:), y(:), gradx(:,:), grady(:,:), & - referenceSR(:), targetSR(:), gxSR(:), gySR(:), refzmn(:), tarzmn(:), & - GradJac(:,:), xiX(:), xiY(:), XiPrime(:,:), targetdef(:,:), residuals(:), & - wtarget(:) -real(kind=wp) :: DD, PCx, PCy, val, err, errmax, rnxi, rnyi, hg(8), W(3,3), gradCIC(8), & - refmean, refstdev, tarmean, tarstdev, CIC, Hessian(8,8), GJ(8), minx, miny, & - xi1max, xi2max, normdp, oldnorm, oldW(3,3), horiginal(8) +real(kind=wp), allocatable :: refpat(:,:), defpat(:,:) +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) 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 -integer(ip),parameter :: kx = 5, ky = 5 ! spline order -integer(ip) :: iflag real(wp) :: tol -! LAPACK variables -character :: UPLO ='U' -integer :: NN = 8 -integer :: NRHS = 1 -real(wp) :: Hessiancopy(8,8) -integer :: LDA = 8 -real(wp) :: SOL(8,1) -integer :: LDB = 8 -integer :: INFO - EMsoft = EMsoft_T( progname, progdesc ) -call setRotationPrecision('d') - -! let's test the second rank tensor transformation routine -alp = 45.D0 * dtor -qu = q_T( qdinp = (/ cos(alp*0.5D0), 0.D0, 0.D0, sin(alp*0.5D0) /) ) -om = qu%qo() -Wnew = om%o_copyd() - -call om%o_print(str = ' rotation matrix : ') - -write(*,*) ' extracted rotation matrix : ' -do i=1,3 - write(*,*) Wnew(i,1:3) -end do - - -! create a symmetric second rank tensor -srt = reshape( (/ 1.D0, 2.D0, 3.D0, 2.D0, 4.D0, 5.D0, 3.D0, 5.D0, 8.D0/), (/ 3,3 /) ) - -write(*,*) ' second rank tensor : ' -do i=1,3 - write(*,*) srt(i,1:3) -end do - -srtrot = TransSecondRankTensor(srt, Wnew, Wnew) - -write(*,*) ' second rank tensor (rotated) : ' -do i=1,3 - write(*,*) srtrot(i,1:3) -end do - -alp = 90.D0 * dtor -qu = q_T( qdinp = (/ cos(alp*0.5D0), 0.D0, 0.D0, sin(alp*0.5D0) /) ) -om = qu%qo() -Wnew = om%o_copyd() -srtrot = TransSecondRankTensor(srt, Wnew) -write(*,*) ' second rank tensor (rotated by single 90°) : ' -do i=1,3 - write(*,*) srtrot(i,1:3) -end do - - - -stop - - ! this program is a first conversion of the IDL HR-EBSD-DIC code into f90; we'll test ! everything here before turning it into a fullblown EMsoftOO program. We will use the IDL ! convention of starting arrays at coordinate 0 for ease of comparison. Also, keep in mind @@ -195,263 +80,84 @@ end function homography2Fe ! pre-processed and scaled between 0.D0 and 1.D0 nx = 640 ny = 480 + +! instantiate the DIC class +! this also initializes the x and y coordinate arrays +DIC = DIC_T( nx, ny ) + nxy = nx * ny allocate(refpat(0:nx-1,0:ny-1), defpat(0:nx-1,0:ny-1)) -! open(20,file='/Users/mdg/playarea/DIC/refpat.data',status='old',form='unformatted') -open(20,file='/Users/mdg/Files/EMPlay/playarea/DIC/refpat.data',status='old',form='unformatted') +open(20,file='/Users/mdg/playarea/DIC/refpat.data',status='old',form='unformatted') +! open(20,file='/Users/mdg/Files/EMPlay/playarea/DIC/refpat.data',status='old',form='unformatted') read(20) refpat close(20,status='keep') +call DIC%setpattern('r', refpat) + ! define some parameters for now... tol = 100 * epsilon(1.0_wp) PCx = real(nx/2-1,wp) PCy = real(ny/2-1,wp) DD = 15000.0_wp/50.0_wp -! generate the reference coordinate grid; these are two separate line arrays at this point -write(*,*) ' creating coordinate arrays x and y' -allocate(x(0:nx-1), y(0:ny-1)) -rnxi = 1.0_wp/real(nx-1,wp) -rnyi = 1.0_wp/real(ny-1,wp) -do i=0,nx-1 - x(i) = real(i,wp) -end do -x = x * rnxi -do j=0,ny-1 - y(j) = real(j,wp) -end do -y = y * rnyi - -! initialize the spline classes stef and sdef -! all interpolations on the deformed region will be done on the whole -! pattern and then the subregion will be extracted for comparison with -! the reference subregion -write(*,*) ' initializing spline coefficients of order ', kx, ky -call sref%initialize(x,y,refpat,kx,ky,iflag) +! define the border widths nbx and nby for the subregion +nbx = 30 +nby = 30 +call DIC%defineSR( nbx, nby, 0.5_wp, 0.5_wp) + + +call DIC%getbsplines(refp=.TRUE., verify=.TRUE.) ! here we deform the reference pattern to obtain a defpat array with known homography ! define the homography hg horiginal = (/ 0.0002_wp, 0.0001_wp, -0.0004_wp, -0.0001_wp, -0.0003_wp, 0.0005_wp, 0.0_wp, 0.0_wp /) -! convert to shape function -W = getShapeFunction(horiginal) -! determine the XiPrime coordinates so that we can apply the deformation by -! means of the evaluate method in the bspline class -allocate( XiPrime(0:1, 0:nx*ny-1)) -cnt = 0 -do j = 0, ny-1 - do i = 0, nx-1 - p2 = matmul( W, (/ x(i)-0.5_wp, y(j)-0.5_wp, 1.0_wp /) ) - if (p2(3).ne.0.0_wp) then - XiPrime(0:1,cnt) = p2(1:2)/p2(3) - else - XiPrime(0:1,cnt) = p2(1:2) - end if - cnt = cnt + 1 - end do -end do -! and interpolate/extrapolate the deformed pattern as needed -call sref%set_extrap_flag(.TRUE.) -cnt = 0 -do j=0,ny-1 - do i=0,nx-1 - iflag = 0_wp - call sref%evaluate(XiPrime(0,cnt)+0.5_wp,XiPrime(1,cnt)+0.5_wp,0,0,defpat(i,j),iflag) - if (iflag.ne.0) write(*,*) 'iflag : ', iflag, i, j - cnt = cnt+1 - end do -end do -! write (*,*) ' defpat: ', minval(defpat), maxval(defpat) -defpat = defpat - minval(defpat) -defpat = defpat/maxval(defpat) -! finally get the b-spline coefficients for the deformed pattern -call sdef%initialize(x,y,defpat,kx,ky,iflag) -deallocate( XiPrime ) - -Woriginal = W +call DIC%applyHomography(horiginal, 0.5_wp, 0.5_wp) + +Woriginal = DIC%getShapeFunction(horiginal) ! get the derivatives of the reference pattern to determine the Hessian -write(*,*) ' computing gradient arrays' -allocate( gradx(0:nx-1,0:ny-1), grady(0:nx-1,0:ny-1) ) -do i=0,nx-1 - do j=0,ny-1 - ! determine how well the spline coefficients represent the original pattern - call sref%evaluate(x(i),y(j),0,0,val,iflag) - err = abs(refpat(i,j)-val) - errmax = max(err,errmax) - ! and extract the gradient arrays - call sref%evaluate(x(i),y(j),1,0,gradx(i,j),iflag) - call sref%evaluate(x(i),y(j),0,1,grady(i,j),iflag) - end do -end do - -! open(20,file='patterns.data',status='unknown',form='unformatted') -! write(20) real(refpat) -! write(20) real(gradx) -! write(20) real(grady) -! close(20) - -! check max error against tolerance -write(*,*) ' max b-spline reconstruction error of reference pattern:', errmax -if (errmax >= tol) then - write(*,*) ' ** test failed ** ' -else - write(*,*) ' ** test passed ** ' -end if -write(*,*) '' -write(*,*) ' range(gradx) : ', minval(gradx), maxval(gradx) -write(*,*) ' range(grady) : ', minval(grady), maxval(grady) +call DIC%getbsplines(defp=.TRUE., grads=.TRUE.) -! define the border widths nbx and nby for the subregion -nbx = 30 -nby = 30 -nxSR = nx-2*nbx -nySR = ny-2*nby -NSR = nxSR * nySR - -! extract the subregions from the reference and gradient arrays -write(*,*) ' creating subregion arrays ' -allocate( referenceSR(0:NSR-1), gxSR(0:NSR-1), gySR(0:NSR-1), refzmn(0:NSR-1), tarzmn(0:NSR-1) ) -referenceSR = reshape( refpat(nbx:nx-nbx-1,nby:ny-nby-1), (/ NSR /) ) -gxSR = reshape( gradx(nbx:nx-nbx-1,nby:ny-nby-1), (/ NSR /) ) -gySR = reshape( grady(nbx:nx-nbx-1,nby:ny-nby-1), (/ NSR /) ) - -! zero-mean and normalize the referenceSR array -write(*,*) ' zero-mean and normalization of subregion arrays ' -refmean = sum(referenceSR)/dble(NSR) -refstdev = sqrt( sum( (referenceSR - refmean)**2 ) ) -refzmn = (referenceSR-refmean)/refstdev -write(*,*) ' reference ', refmean, refstdev, sum(refzmn)/dble(NSR), & - sqrt( sum( (refzmn - sum(refzmn)/dble(NSR))**2 ) ) - -! repeat for the targetSR so that we can compute the initial CIC -tarmean = sum(defpat(nbx:nx-nbx-1,nby:ny-nby-1))/dble(NSR) -tarstdev = sqrt( sum( (defpat(nbx:nx-nbx-1,nby:ny-nby-1) - tarmean)**2 ) ) -tarzmn = reshape( (defpat(nbx:nx-nbx-1,nby:ny-nby-1)-tarmean)/tarstdev, (/NSR /) ) -write(*,*) ' target ', tarmean, tarstdev, sum(tarzmn)/dble(NSR), & - sqrt( sum( (tarzmn - sum(tarzmn)/dble(NSR))**2 ) ) +! zero-mean and normalize the referenceSR and targetSR arrays +call DIC%applyZMN(doreference=.TRUE., dotarget=.TRUE.) ! initial CIC function (to be minimized ...) -CIC = sum( (refzmn-tarzmn)**2 ) +call DIC%getresiduals( CIC ) write (*,*) ' initial CIC value : ', CIC -! next we compute the GradJac array and the Hessian using the standard grid shifted -! by the Pattern Center coordinates -write(*,*) ' computing GradJac and Hessian arrays' -allocate( xiX(0:nxSR-1), xiY(0:nySR-1) ) -xiX = x(nbx:nx-nbx-1) - 0.5_wp ! PCx -xiY = y(nby:ny-nby-1) - 0.5_wp ! PCy -write(*,*) ' range(xiX) : ', minval(xiX), maxval(xiX), shape(xiX), PCx -write(*,*) ' range(xiY) : ', minval(xiY), maxval(xiY), shape(xiY), PCy - -allocate( GradJac(0:7,0:NSR-1) ) -cnt = 0 -Hessian = 0.0_wp -do j = 0, nySR-1 - do i = 0, nxSR-1 - GJ = (/ gxSR(cnt)*xiX(i), gxSR(cnt)*xiY(j), gxSR(cnt), & - gySR(cnt)*xiX(i), gySR(cnt)*xiY(j), gySR(cnt), & - gxSR(cnt)*xiX(i)**2-gySR(cnt)*xiX(i)*xiY(j), & - gySR(cnt)*xiY(j)**2-gxSR(cnt)*xiX(i)*xiY(j) /) - GradJac(0:7,cnt) = GJ - Hessian = Hessian + matmul( reshape(GJ, (/8,1/)), reshape(GJ, (/1,8/)) ) - cnt = cnt + 1 - end do -end do -Hessian = Hessian * 2.0_wp / refstdev**2 - -write(*,*) ' final cnt : ', cnt-1, NSR - -write(*,*) ' Hessian :' -write(*,*) Hessian - -write(*,*) ' range(GradJac) : ', minval(GradJac), maxval(GradJac) +! compute the Hessian matrix +call DIC%getHessian( Hessian ) ! 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 = getShapeFunction(hg) -write (*,*) W - -write (*,*) getHomography(W) - -allocate( targetdef(0:nxSR-1,0:nySR-1), residuals(0:NSR-1), wtarget(0:NSR-1) ) - -! allow for extrapolation in addition to interpolation... -! this required a change to the original bspline-fortran code; the -! modified code is in a personal fork of the package which is -! installed during the EMsoftOO super build process. -call sdef%set_extrap_flag(.TRUE.) - -dx = abs(xiX(2) - xiX(1)) -dy = abs(xiY(2) - xiY(1)) - -allocate( XiPrime(0:1, 0:NSR-1) ) +W = DIC%getShapeFunction(hg) oldnorm = 100.0_wp ! and here we start the loop do ii=1,100 write (*,*) ' iteration # ',ii - ! get XiPrime - cnt = 0 - do j = 0, nySR-1 - do i = 0, nxSR-1 - p2 = matmul( W, (/ xiX(i), xiY(j), 1.0_wp /) ) - XiPrime(0:1,cnt) = p2(1:2)/p2(3) - cnt = cnt + 1 - end do - end do - -! apply this deformation to the target pattern using b-splines - targetdef = 0.D0 - cnt = 0 - do j=0,nySR-1 - do i=0,nxSR-1 - iflag = 0_wp - ! call sdef%evaluate(XiPrime(0,cnt)-minx,XiPrime(1,cnt)-miny,0,0,targetdef(i,j),iflag) - call sdef%evaluate(XiPrime(0,cnt)+0.5_wp,XiPrime(1,cnt)+0.5_wp,0,0,targetdef(i,j),iflag) - if (iflag.ne.0) write(*,*) 'iflag : ', iflag, i, j - cnt = cnt+1 - end do - end do - - ! extract the subregion - wtarget = reshape( targetdef, (/ NSR /) ) - ! zero-mean and normalize - tarmean = sum(wtarget)/dble(NSR) - tarstdev = sqrt( sum( (wtarget - tarmean)**2 ) ) - tarzmn = (wtarget-tarmean)/tarstdev + call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.) + + ! zero-mean and normalize the deformed sub-region array wtarget + call DIC%applyZMN( dotarget = .TRUE. ) + ! residuals - residuals = refzmn - tarzmn - CIC = sum( residuals**2 ) - write (*,*) ' range(residuals) : ', minval(residuals), maxval(residuals) - write (*,*) ' new CIC value : ', CIC - ! gradCIC - gradCIC = matmul(GradJac, residuals) - gradCIC = gradCIC * 2.0_wp/refstdev - write (*,*) ' gradCIC :', shape(GradJac), shape(residuals), shape(gradCIC) - write (*,*) gradCIC - ! solve the Hessian equation using Cholesky decomposition (LAPACK) - Hessiancopy = Hessian - SOL(1:8,1) = -gradCIC(1:8) - call DPOSV(UPLO, NN, NRHS, Hessiancopy, LDA, SOL, LDB, INFO) - ! extract the solution vector from SOL - write (*,*) ' DPOSV INFO = ', INFO - write (*,*) ' Solution vector : ' - write (*,*) SOL - xi1max = maxval( XiPrime(0,:) ) - xi2max = maxval( XiPrime(1,:) ) - normdp = sqrt(xi1max * (SOL(1,1)**2+SOL(4,1)**2+SOL(7,1)**2) + & - xi2max * (SOL(2,1)**2+SOL(5,1)**2+SOL(8,1)**2) + SOL(3,1)**2 + SOL(6,1)**2) + call DIC%getresiduals( CIC ) + + ! gradCIC and solve equation + call DIC%solveHessian(SOL, normdp) ! convert to a shape function and take the inverse - Wnew = getShapeFunction(reshape(SOL, (/8/))) - Winv = matrixInvert( Wnew ) + Wnew = DIC%getShapeFunction(reshape(SOL, (/8/))) + Winv = matrixInvert_wp( Wnew ) W = matmul( W, Winv ) W = W / W(3,3) - write (*,*) getHomography(W) + hg = DIC%getHomography(W) + write (*,*) hg write (*,*) ' norm(deltap) = ', normdp write (*,*) '------------' if (normdp.lt.oldnorm) then @@ -463,158 +169,17 @@ end function homography2Fe end if end do -W = matrixInvert( W ) +W = matrixInvert_wp( W ) W = W / W(3,3) +hg = DIC%getHomography(W) write (*,*) '' write (*,*) ' Final homography : ' -write (*,*) getHomography(W) +write (*,*) DIC%getHomography(W) write (*,*) '' write (*,*) ' Target homography : ' write (*,*) horiginal - -call sdef%destroy() -call sref%destroy() +write (*,*) ' differences : ' +write (*,*) horiginal-hg end program EMplay - - - -! functions and subroutines ... - - - - -function getShapeFunction(h) result(W) - -use bspline_kinds_module, only: wp - -real(wp), INTENT(IN) :: h(8) -real(wp) :: W(3,3) - -W(1,1) = 1.0_wp + h(1) -W(2,2) = 1.0_wp + h(5) -W(3,3) = 1.0_wp - -W(1,2:3) = h(2:3) -W(2,1) = h(4) -W(2,3) = h(6) -W(3,1:2) = h(7:8) - -end function getShapeFunction - - -function getHomography(W) result(h) - -use bspline_kinds_module, only: wp - -real(wp), INTENT(IN) :: W(3,3) -real(wp) :: h(8) - -h = (/ W(1,1)-1.0_wp, W(1,2), W(1,3), W(2,1), W(2,2)-1.0_wp, W(2,3), W(3,1), W(3,2) /) - -end function getHomography - - - - -function matrixInvert(a) result(b) - -use bspline_kinds_module, only: wp -use mod_io - -real(wp), INTENT(IN) :: a(3,3) -real(wp) :: b(3,3) - -type(IO_T) :: Message -real(wp) :: d - -d = a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+ & - a(1,3)*a(2,1)*a(3,2)-a(1,3)*a(2,2)*a(3,1)- & - a(1,2)*a(2,1)*a(3,3)-a(1,1)*a(2,3)*a(3,2) -if (d.ne.0.D0) then - b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2) - b(1,2)=a(1,3)*a(3,2)-a(1,2)*a(3,3) - b(1,3)=a(1,2)*a(2,3)-a(1,3)*a(2,2) - b(2,1)=a(2,3)*a(3,1)-a(2,1)*a(3,3) - b(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1) - b(2,3)=a(1,3)*a(2,1)-a(1,1)*a(2,3) - b(3,1)=a(2,1)*a(3,2)-a(2,2)*a(3,1) - b(3,2)=a(1,2)*a(3,1)-a(1,1)*a(3,2) - b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1) - b = b/d -else - do i=1,3 - write (*,*) (a(i,j),j=1,3) - end do - ! call Message%printError('mInvert','matrix has zero determinant') - call Message%printMessage('mInvert: matrix has zero determinant') - b = a -end if - -end function matrixInvert - - -function homography2Fe(p, PC, DD) result(W) - -use bspline_kinds_module, only: wp - -real(wp), INTENT(IN) :: p(0:7) -real(wp), INTENT(IN) :: PC(0:1) -real(wp), INTENT(IN) :: DD -real(wp) :: W(3,3) - -real(wp) :: beta0, Fe11, Fe12, Fe13, Fe21, Fe22, Fe23, Fe31, Fe32, Fe33 - -! h11, h12, h13, h21, h22, h23, h31, h32 -! 0 1 2 3 4 5 6 7 - -beta0 = 1.0_wp - p(6) * PC(0) - p(7) * PC(1) -Fe11 = 1.0_wp + p(0) + p(6) * PC(0) -Fe12 = p(1) + p(7) * PC(0) -Fe13 = (p(2) - p(0)*PC(0) - p(1)*PC(1) + PC(0)*(beta0 - 1.0_wp))/DD -Fe21 = p(3) + p(6) * PC(1) -Fe22 = 1.0_wp + p(4) + p(7) * PC(1) -Fe23 = (p(5) - p(3)*PC(0) - p(4)*PC(1) + PC(1)*(beta0 - 1.0_wp))/DD -Fe31 = DD * p(6) -Fe32 = DD * p(7) -Fe33 = beta0 - -W = reshape( (/ Fe11, Fe12, Fe13, Fe21, Fe22, Fe23, Fe31, Fe32, Fe33 /), (/ 3,3 /)) / beta0 - -end function homography2Fe - - -function Fe2homography(Fe_input, PC, DD) result(h) - -use bspline_kinds_module, only: wp - -real(wp), INTENT(IN) :: Fe_input(0:2,0:2) -real(wp), INTENT(IN) :: PC(0:1) -real(wp), INTENT(IN) :: DD -real(wp) :: h(8) - -real(wp) :: Fe(0:2,0:2), g0, g11, g22, g13, g23, h11, h12, h13, h21, h22, h23, h31, h32 - -Fe = Fe_input/Fe_input(2,2) - -g0 = DD + Fe(2,0) * PC(0) + Fe(2,1) * PC(1) -g11 = DD * Fe(0,0) - Fe(2,0) * PC(0) - g0 -g22 = DD * Fe(1,1) - Fe(2,1) * PC(1) - g0 -g13 = DD * ((Fe(0,0) - 1.D0) * PC(0) + Fe(0,1) * PC(1) + Fe(0,2) * DD) + PC(0) * (DD - g0) -g23 = DD * (Fe(1,0) * PC(0) + (Fe(1,1) - 1.D0) * PC(1) + Fe(1,2) * DD) + PC(1) * (DD - g0) -h11 = g11 -h12 = DD * Fe(0,1) - Fe(2,1) * PC(0) -h13 = g13 -h21 = DD * Fe(1,0) - Fe(2,0) * PC(1) -h22 = g22 -h23 = g23 -h31 = Fe(2,0) -h32 = Fe(2,1) - -h = (/ h11, h12, h13, h21, h22, h23, h31, h32 /) / g0 - -end function Fe2homography - - -