new program EMppEBSD to pre-process a file of experimental EBSD patterns

README.md

@@ -67,8 +67,10 @@ If you do not need the complete EMsoftOO package, you can compile sections of th
 
 ## Licenses ##
 
+The BSD-3 license text below can be found at the start of every source file in the EMsoftOO package.
+
 	!###################################################################
-	! Copyright (c) 2013-2022, Marc De Graef Research Group/Carnegie Mellon University
+	! Copyright (c) 2013-2025, Marc De Graef Research Group/Carnegie Mellon University
 	! All rights reserved.
 	!
 	! Redistribution and use in source and binary forms, with or without modification, are 

Source/EMsoftOOLib/CMakeLists.txt

@@ -164,7 +164,7 @@ set(EMsoftOOLib_SRCS
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_ISEDI.f90
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_EBSDPCA.f90
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_EBSDdefect.f90
-#  ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_HREBSD.f90
+  ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_HREBSDDIC.f90
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_CBED.f90
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_SRdefect.f90 
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_SRCBED.f90 
@@ -175,6 +175,7 @@ set(EMsoftOOLib_SRCS
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_BWshow.f90 
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_SphInx.f90
   ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_4DEBSD.f90
+  ${EMsoftOOLib_SOURCE_DIR}/program_mods/mod_ppEBSD.f90 
 
   ${EMsoftOOLib_Additional_SRCS}
 )

Source/EMsoftOOLib/mod_DIC.f90

@@ -44,9 +44,9 @@ module mod_DIC
 
 type, public :: DIC_T
   private
-    real(wp)                        :: W(3,3)
-    real(wp)                        :: h(8)
-    real(wp)                        :: CIC
+    real(wp)                        :: W(3,3)     ! shape function
+    real(wp)                        :: h(8)       ! homography
+    real(wp)                        :: CIC        ! current value of function to be minimized
     integer(ip)                     :: kx = 5     ! spline order
     integer(ip)                     :: ky = 5     ! spline order
     integer(kind=irg)               :: nx         ! pattern x-size
@@ -58,12 +58,12 @@ module mod_DIC
     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)
+    type(bspline_2d)                :: sref       ! b-spline class for reference pattern
+    type(bspline_2d)                :: sdef       ! b-spline class for deformed pattern
+    real(wp)                        :: tol = 100.0_wp * epsilon(1.0_wp) ! tolerance
     logical                         :: verbose = .FALSE.   ! useful for debugging
 
-  ! the arrays that are tied to a given reference/target pattern are defined here
+! 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(:)
@@ -105,6 +105,7 @@ module mod_DIC
     procedure, pass(self) :: getresiduals_
     procedure, pass(self) :: setverbose_
     procedure, pass(self) :: setpattern_
+    procedure, pass(self) :: getpattern_
 
     final :: DIC_destructor
 
@@ -121,6 +122,7 @@ module mod_DIC
     generic, public :: getresiduals => getresiduals_
     generic, public :: setverbose => setverbose_
     generic, public :: setpattern => setpattern_
+    generic, public :: getpattern => getpattern_
 
 end type DIC_T
 
@@ -184,6 +186,25 @@ subroutine DIC_destructor( self )
 call self%sref%destroy()
 call self%sdef%destroy()
 
+if (allocated(self%x)) deallocate(self%x)
+if (allocated(self%y)) deallocate(self%y)
+if (allocated(self%xiX)) deallocate(self%xiX)
+if (allocated(self%xiY)) deallocate(self%xiY)
+if (allocated(self%XiPrime)) deallocate(self%XiPrime)
+if (allocated(self%refpat)) deallocate(self%refpat)
+if (allocated(self%defpat)) deallocate(self%defpat)
+if (allocated(self%wtarget)) deallocate(self%wtarget)
+if (allocated(self%gradx)) deallocate(self%gradx)
+if (allocated(self%grady)) deallocate(self%grady)
+if (allocated(self%gxSR)) deallocate(self%gxSR)
+if (allocated(self%gySR)) deallocate(self%gySR)
+if (allocated(self%GradJac)) deallocate(self%GradJac)
+if (allocated(self%referenceSR)) deallocate(self%referenceSR)
+if (allocated(self%targetSR)) deallocate(self%targetSR)
+if (allocated(self%refzmn)) deallocate(self%refzmn)
+if (allocated(self%tarzmn)) deallocate(self%tarzmn)
+if (allocated(self%residuals)) deallocate(self%residuals)
+
 call reportDestructor('DIC_T')
 
 end subroutine DIC_destructor
@@ -242,13 +263,42 @@ recursive subroutine setpattern_(self, rp, pattern)
   self%refpat = pattern 
 end if
 
-if (rp.eq.'p') then 
+if (rp.eq.'d') then 
   allocate( self%defpat(0:sz(1)-1,0:sz(2)-1) )
   self%defpat = pattern 
 end if
 
 end subroutine setpattern_
 
+!--------------------------------------------------------------------------
+recursive function getpattern_(self, rp, nx, ny) result(pattern)
+!DEC$ ATTRIBUTES DLLEXPORT :: getpattern_
+ !! author: MDG
+ !! version: 1.0
+ !! date: 12/01/24
+ !!
+ !! get the reference or target pattern
+
+use mod_IO 
+
+IMPLICIT NONE
+
+class(DIC_T),INTENT(INOUT)      :: self
+character(1),INTENT(IN)         :: rp
+integer(kind=irg),INTENT(IN)    :: nx
+integer(kind=irg),INTENT(IN)    :: ny
+real(kind=sgl)                  :: pattern(nx, ny)
+
+if (rp.eq.'r') then 
+  pattern = real(self%refpat)
+end if
+
+if (rp.eq.'d') then 
+  pattern = real(self%defpat)
+end if
+
+end function getpattern_
+
 !--------------------------------------------------------------------------
 recursive subroutine getbsplines_(self, verify, refp, defp, grads)
 !DEC$ ATTRIBUTES DLLEXPORT :: getbsplines_
@@ -300,7 +350,7 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads)
   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
+            ! determine how well the spline coefficients represent the original reference 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)
@@ -319,7 +369,9 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads)
 
 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) )
+    if (.not.allocated(self%gradx)) then 
+      allocate( self%gradx(0:self%nx-1,0:self%ny-1), self%grady(0:self%nx-1,0:self%ny-1) )
+    end if
     do i=0,self%nx-1
        do j=0,self%ny-1
             ! extract the gradient arrays

Source/EMsoftOOLib/stringconstants.in

@@ -1,5 +1,5 @@
 ! ###################################################################
-! Copyright (c) 2017-2022, Marc De Graef Research Group/Carnegie Mellon University
+! Copyright (c) 2017-2025, Marc De Graef Research Group/Carnegie Mellon University
 ! All rights reserved.
 !
 ! Redistribution and use in source and binary forms, with or without modification, are
@@ -312,6 +312,9 @@ OSM;"OSM"
 OpenCLpathname;"OpenCLpathname"
 Operator;"Operator"
 OrderParameter;"CPLP"
+ppEBSD;"ppEBSD"
+ppEBSDNameList;"ppEBSDNameList"
+ppEBSDNML;"ppEBSDNML"
 PEDZA;"PEDZA"
 PEDZANML;"PEDZANML"
 PEDZANameList;"PEDZANameList"

Source/SEM/CMakeLists.txt

@@ -278,6 +278,15 @@ if(EMsoftOO_ENABLE_HDF5_SUPPORT)
                         INCLUDE_DIRS ${EMsoftOOLib_BINARY_DIR}
                         )
 
+  Add_EMsoftOO_Executable(TARGET EMppEBSD
+                        SOURCES ${APP_DIR}/EMppEBSD.f90
+                        LINK_LIBRARIES ${EXE_LINK_LIBRARIES}
+                        TEMPLATE ${TMPLT_DIR}/EMppEBSD.template
+                        SOLUTION_FOLDER EMsoftOOPublic/SEM
+                        INSTALL_PROGRAM TRUE
+                        INCLUDE_DIRS ${EMsoftOOLib_BINARY_DIR}
+                        )
+
   Add_EMsoftOO_Executable(TARGET EMkinematical
                         SOURCES ${APP_DIR}/EMkinematical.f90
                         LINK_LIBRARIES ${EXE_LINK_LIBRARIES}

Source/SEM/EMppEBSD.f90

@@ -0,0 +1,61 @@
+! ###################################################################
+! Copyright (c) 2016-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.
+! ###################################################################
+
+program EMppEBSD
+  !! author: MDG
+  !! version: 1.0 
+  !! date: 12/06/24
+  !!
+  !! Standard pre-processing of a file containing EBSD patterns; any format is accepted
+  !! output will be a binary file... can be a very large file ...
+  !!
+
+use mod_kinds
+use mod_global
+use mod_EMsoft
+use mod_ppEBSD
+use stringconstants
+
+IMPLICIT NONE
+
+character(fnlen)        :: progname = 'EMppEBSD.f90'
+character(fnlen)        :: progdesc = 'Pre-processing of an EBSD/ECP/TKD pattern data set'
+
+type(EMsoft_T)          :: EMsoft
+type(ppEBSD_T)          :: ppEBSD
+
+! print the EMsoft header and handle any command line arguments  
+EMsoft = EMsoft_T( progname, progdesc, tpl = (/ 104 /) )
+
+! deal with the namelist stuff
+ppEBSD = ppEBSD_T(EMsoft%nmldeffile)
+
+! perform the pattern computations
+call ppEBSD%ppEBSD(EMsoft, progname)
+
+end program EMppEBSD

Source/TestPrograms/play.f90

@@ -31,7 +31,7 @@ program EMplay
   !! version: 1.0 
   !! date: 02/26/24
   !!
-  !! test program to look at triple junctions and fundamental zones
+  !! test program 
 
 use mod_kinds
 use mod_global
@@ -62,10 +62,12 @@ program EMplay
 type(e_T)                    :: eu 
 
 real(kind=wp), allocatable   :: refpat(:,:), defpat(:,:)
+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)
+                                minx, miny, xi1max, xi2max, normdp, oldnorm, oldW(3,3), horiginal(8), CIC, sol(8), &
+                                homographies(8,1000)
 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(kind=irg)            :: nx, ny, nxy, nbx, nby, i, ii, j, NSR, cnt, nxSR, nySR, jj  
 real(wp)                     :: tol
 
 
@@ -88,8 +90,8 @@ program EMplay
 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')
 
@@ -109,11 +111,30 @@ program EMplay
 
 call DIC%getbsplines(refp=.TRUE., verify=.TRUE.)
 
+! read 1000 random homographies from a file and run the fitting, then 
+! write result to another file for comparison with an IDL routine.
+open(unit=dataunit,file='/Volumes/Drive2/playarea/DIC/test/homographies.data',status='unknown',form='unformatted')
+read(dataunit) homographies
+close(dataunit,status='keep')
+
+open(unit=dataunit,file='/Volumes/Drive2/playarea/DIC/test/homography-fits2.data',status='unknown',form='unformatted')
+
+do jj=1,1000
 ! 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 /)
+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 /)
 call DIC%applyHomography(horiginal, 0.5_wp, 0.5_wp)
 
+! store both patterns in a binary file 
+! allocate( tmppat( nx, ny ) )
+! open(unit=dataunit,file='patterns.data',status='unknown',form='unformatted')
+! tmppat = DIC%getpattern('r', nx, ny)
+! write(dataunit) tmppat
+! tmppat = DIC%getpattern('d', nx, ny)
+! write(dataunit) tmppat
+! close(dataunit,status='keep')
+! deallocate(tmppat)
+
 Woriginal = DIC%getShapeFunction(horiginal)
 
 ! get the derivatives of the reference pattern to determine the Hessian
@@ -155,12 +176,12 @@ 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
+    if (normdp.lt.1.1*oldnorm) then
         oldnorm = normdp
         oldW = W 
     else
@@ -170,16 +191,25 @@ program EMplay
 end do 
 
 W = matrixInvert_wp( W )
-W = W / W(3,3)
+!W = W / W(3,3)
 hg = DIC%getHomography(W)
 write (*,*) '' 
 write (*,*) ' Final homography : '
 write (*,*) DIC%getHomography(W)
 write (*,*) '' 
 write (*,*) ' Target homography : '
 write (*,*) horiginal
-
-write (*,*) ' differences : '
+write (*,*) '' 
+write (*,*) ' Differences : '
 write (*,*) horiginal-hg
 
+! write results to data file (single precision because IDL has a bug for double precision)
+write (dataunit) real(hg) 
+write (dataunit) real(normdp)
+write (dataunit) ii
+
+end do ! loop over jj 
+
+close(dataunit,status='keep')
+
 end program EMplay

resources/templatecodes.txt

@@ -79,6 +79,7 @@
 101EMgrainviz
 102EMCBED2DQC
 103EMSRCBED
+104EMppEBSD
 105EMDIRAM
 106EMEBSDpc
 107EMEBSDPCA
@@ -107,6 +108,7 @@
 282EMEBSDdefect
 287EMHREBSD
 288EMHREBSDpreview
+289EMHREBSDDIC
 290EMgetIPF
 291EMorav
 300EMShapeAmplitude