diff --git a/NamelistTemplates/EMHREBSDDIC.template b/NamelistTemplates/EMHREBSDDIC.template new file mode 100644 index 0000000..39ee511 --- /dev/null +++ b/NamelistTemplates/EMHREBSDDIC.template @@ -0,0 +1,38 @@ + &HREBSDDICdata +!################################################################### +! EBSD Pattern Parameters will be read from an EMDI output file +!################################################################### +! name of the dot product file (.h5) path relative to EMdatapathname + DIfile = 'undefined', +! output HDF5 file; path relative to EMdatapathname + datafile = 'undefined', +! relative path to the namelist file for EMppEBSD program + ppEBSDnml = 'undefined' +! coordinates of the reference pattern to be used (0-based coordinates) + patx = 0, + paty = 0, +! width of the exclusion zone surrounding the sub-region + nbx = 0, + nby = 0, +! number of threads (default = 1) + nthreads = 1, +! +! name of experimental pattern file; path relative to EMdatapathname +! for this program, this should be a preprocessed data file (EMppEBSD output) + patternfile = 'undefined', + +!################################################################### +! Stiffness tensor parameters (unit:GPa) +!################################################################### +! type of crystal (cubic: 'cub' or hcp: 'hex') +crystal= 'cub', + +! cubic +C11 = 276.0, +C12 = 159.0, +C44 = 132.0, + +! additional terms for HCP +C13 = 0.0 +C33 = 0.0 + / diff --git a/NamelistTemplates/EMppEBSD.template b/NamelistTemplates/EMppEBSD.template new file mode 100644 index 0000000..3372efe --- /dev/null +++ b/NamelistTemplates/EMppEBSD.template @@ -0,0 +1,51 @@ + &ppEBSDdata +! The line above must not be changed +! +! The values below are the default values for this program +! +! width of data set in pattern input file + ipf_wd = 100, +! height of data set in pattern input file + ipf_ht = 100, +! define the region of interest as x0 y0 w h; leave all at 0 for full field of view +! region of interest has the point (x0,y0) as its upper left corner and is w x h patterns + ROI = 0 0 0 0, +! to use a custom mask, enter the mask filename here; leave undefined for standard mask option + maskfile = 'undefined', +! mask or not + maskpattern = 'n', +! mask radius (in pixels, AFTER application of the binning operation) + maskradius = 240, +! hi pass filter w parameter; 0.05 is a reasonable value + hipassw = 0.05, +! number of regions for adaptive histogram equalization + nregions = 10, +! pattern dimensions (actual) + numsx = 640, + numsy = 480, +! +!################################################################### +! INPUT FILE PARAMETERS: COMMON TO 'STATIC' AND 'DYNAMIC' +!################################################################### +! +! name of datafile where the patterns are stored; path relative to EMdatapathname + exptfile = 'undefined', +! input file type parameter: Binary, EMEBSD, EMEBSD32i, EMEBSD32f, TSLHDF, TSLup2, +! OxfordHDF, OxfordBinary, BrukerHDF, NORDIF + inputtype = 'Binary', +! here we enter the HDF group names and data set names as individual strings (up to 10) +! enter the full path of a data set in individual strings for each group, in the correct order, +! and with the data set name as the last name; leave the remaining strings empty (they should all +! be empty for the Binary and TSLup2 formats) + HDFstrings = '' '' '' '' '' '' '' '' '' '', +! +!################################################################### +! OTHER FILE PARAMETERS: COMMON TO 'STATIC' AND 'DYNAMIC' +!################################################################### +! +! temporary data storage file name ; path will be relative to EMdatapathname + tmpfile = 'EMEBSD_pp.data', +! +! number of threads for parallel pre-processing (0 to use all available) + nthreads = 1, + / diff --git a/Source/EMsoftOOLib/mod_DIC.f90 b/Source/EMsoftOOLib/mod_DIC.f90 index 4fae037..1575759 100644 --- a/Source/EMsoftOOLib/mod_DIC.f90 +++ b/Source/EMsoftOOLib/mod_DIC.f90 @@ -106,6 +106,7 @@ module mod_DIC procedure, pass(self) :: setverbose_ procedure, pass(self) :: setpattern_ procedure, pass(self) :: getpattern_ + procedure, pass(self) :: DIC_cleanup_ final :: DIC_destructor @@ -123,6 +124,7 @@ module mod_DIC generic, public :: setverbose => setverbose_ generic, public :: setpattern => setpattern_ generic, public :: getpattern => getpattern_ + generic, public :: cleanup => DIC_cleanup_ end type DIC_T @@ -134,7 +136,7 @@ module mod_DIC contains !-------------------------------------------------------------------------- -type(DIC_T) function DIC_constructor( nx, ny ) result(DIC) +recursive type(DIC_T) function DIC_constructor( nx, ny ) result(DIC) !DEC$ ATTRIBUTES DLLEXPORT :: DIC_constructor !! author: MDG !! version: 1.0 @@ -209,6 +211,28 @@ subroutine DIC_destructor( self ) end subroutine DIC_destructor +!-------------------------------------------------------------------------- +subroutine DIC_cleanup_( self ) +!DEC$ ATTRIBUTES DLLEXPORT :: DIC_cleanup_ + !! author: MDG + !! version: 1.0 + !! date: 12/11/24 + !! + !! deallocate some arrays + +IMPLICIT NONE + +class(DIC_T), INTENT(INOUT) :: self + +call self%sdef%destroy() + +if (allocated(self%defpat)) deallocate(self%defpat) +if (allocated(self%wtarget)) deallocate(self%wtarget) +if (allocated(self%tarzmn)) deallocate(self%tarzmn) +if (allocated(self%residuals)) deallocate(self%residuals) + +end subroutine DIC_cleanup_ + !-------------------------------------------------------------------------- recursive subroutine setverbose_(self, v) !DEC$ ATTRIBUTES DLLEXPORT :: setverbose_ @@ -254,18 +278,21 @@ recursive subroutine setpattern_(self, rp, pattern) character(1),INTENT(IN) :: rp real(wp),INTENT(IN) :: pattern(:,:) +type(IO_T) :: Message integer(kind=irg) :: sz(2) sz = shape(pattern) if (rp.eq.'r') then - allocate( self%refpat(0:sz(1)-1,0:sz(2)-1) ) + if (.not.allocated(self%defpat)) allocate( self%refpat(0:sz(1)-1,0:sz(2)-1) ) self%refpat = pattern + if (self%verbose) call Message%printMessage(' setpattern_::reference pattern set ') end if if (rp.eq.'d') then - allocate( self%defpat(0:sz(1)-1,0:sz(2)-1) ) + if (.not.allocated(self%defpat)) allocate( self%defpat(0:sz(1)-1,0:sz(2)-1) ) self%defpat = pattern + if (self%verbose) call Message%printMessage(' setpattern_::target pattern set ') end if end subroutine setpattern_ @@ -289,12 +316,16 @@ recursive function getpattern_(self, rp, nx, ny) result(pattern) integer(kind=irg),INTENT(IN) :: ny real(kind=sgl) :: pattern(nx, ny) +type(IO_T) :: Message + if (rp.eq.'r') then pattern = real(self%refpat) + if (self%verbose) call Message%printMessage(' getpattern_::reference pattern retrieved ') end if if (rp.eq.'d') then pattern = real(self%defpat) + if (self%verbose) call Message%printMessage(' getpattern_::target pattern retrieved ') end if end function getpattern_ @@ -332,7 +363,7 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads) 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') + if (self%verbose) call Message%printMessage(' getbsplines_::b-splines for reference pattern completed') end if end if @@ -342,7 +373,7 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads) 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') + if (self%verbose) call Message%printMessage(' getbsplines_::b-splines for target pattern completed') end if end if @@ -358,7 +389,7 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads) 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) + call Message%WriteValue(' getbsplines_::max b-spline reconstruction error of reference pattern:', io_real, 1) if (errmax >= self%tol) then call Message%printMessage(' ** reference reconstruction test failed ** ') else @@ -379,7 +410,11 @@ recursive subroutine getbsplines_(self, verify, refp, defp, grads) 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') + if (self%verbose) then + call Message%printMessage(' getbsplines_::gradient of reference pattern completed') + write (*,*) ' getbsplines_::gradx(1:2,1:2) : ', self%gradx(1:2,1:2) + write (*,*) ' getbsplines_::grady(1:2,1:2) : ', self%grady(1:2,1:2) + end if 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 @@ -419,22 +454,25 @@ recursive subroutine defineSR_(self, nbx, nby, PCx, PCy) 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 /) ) +if (self%verbose) call Message%printMessage(' defineSR_::referenceSR array allocated') ! 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 +if (self%verbose) call Message%printMessage(' defineSR_::xiX, xiY arrays allocated') ! allocate array for the product of the gradient and the Jacobian allocate( self%GradJac(0:7,0:self%NSR-1) ) +if (self%verbose) call Message%printMessage(' defineSR_::GradJac array allocated') ! allocate the residuals array allocate( self%residuals(0:self%NSR-1) ) +if (self%verbose) call Message%printMessage(' defineSR_::residuals array allocated') ! and finally the targetdef array allocate( self%targetSR(0:self%NSR-1) ) - -if (self%verbose) call Message%printMessage(' sub-region arrays allocated') +if (self%verbose) call Message%printMessage(' defineSR_::targetSR array allocated') end subroutine defineSR_ @@ -464,7 +502,8 @@ recursive subroutine applyZMN_(self, doreference, dotarget) 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) + if (self%verbose) call Message%WriteValue(' applyZMN_::reference mean/stdev : ', io_real, 2) + ! write (*,*) ' applyZMN_::refzmn range ', minval(self%refzmn), maxval(self%refzmn) end if end if @@ -474,7 +513,8 @@ recursive subroutine applyZMN_(self, doreference, dotarget) 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) + if (self%verbose) call Message%WriteValue(' applyZMN_::target mean/stdev : ', io_real, 2) + ! write (*,*) ' applyZMN_::tarzmn range ', minval(self%tarzmn), maxval(self%tarzmn) end if end if @@ -535,6 +575,9 @@ recursive subroutine applyHomography_(self, h, PCx, PCy, dotarget) cnt = cnt + 1 end do end do +if (self%verbose) then + write (*,*) ' applyHomography_::XiPrime(0:1,0:1) : ', self%XiPrime(0:1,0:1) +end if ! and interpolate/extrapolate the deformed pattern as needed if (.not.allocated(self%defpat)) allocate( self%defpat( 0:lnx-1, 0:lny-1 )) @@ -553,10 +596,15 @@ recursive subroutine applyHomography_(self, h, PCx, PCy, dotarget) ! reduce to interval [0,1] self%defpat = self%defpat - minval(self%defpat) self%defpat = self%defpat/maxval(self%defpat) +if (self%verbose) then + call Message%printMessage(' applyHomography_::defpat interpolated') + write (*,*) ' applyHomography_::range(defpat) : ',minval(self%defpat), maxval(self%defpat) +end if ! 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 @@ -585,7 +633,7 @@ recursive subroutine getresiduals_(self, CIC) real(wp),INTENT(OUT),OPTIONAL :: CIC type(IO_T) :: Message -real(kind=dbl) :: io_real(1) +real(kind=dbl) :: io_real(2) self%residuals = self%refzmn-self%tarzmn @@ -594,7 +642,10 @@ recursive subroutine getresiduals_(self, CIC) CIC = self%CIC if (self%verbose) then io_real(1) = dble(CIC) - call Message%WriteValue(' CIC value : ', io_real, 1) + call Message%WriteValue(' getresiduals_::CIC value : ', io_real, 1) + io_real(1) = minval(self%residuals) + io_real(2) = maxval(self%residuals) + call Message%WriteValue(' getresiduals_::residuals range : ', io_real, 2) end if end if @@ -637,7 +688,10 @@ recursive subroutine getHessian_(self, Hessian) if (present(Hessian)) Hessian = self%Hessian -if (self%verbose) call Message%printMessage(' reference Hessian computed') +if (self%verbose) then + call Message%printMessage(' getHessian_::reference Hessian computed') + write (*,*) Hessian +end if end subroutine getHessian_ @@ -658,8 +712,11 @@ recursive subroutine solveHessian_(self, SOL, normdp) real(wp),INTENT(INOUT) :: SOL(8) real(wp),INTENT(INOUT) :: normdp +type(IO_T) :: Message + real(wp) :: Hessiancopy(8,8) ! local copy real(wp) :: xi1max, xi2max +real(kind=dbl) :: io_real(2) ! LAPACK variables character :: UPLO ='U' @@ -673,18 +730,26 @@ recursive subroutine solveHessian_(self, SOL, normdp) ! compute gradCIC (right-hand side of the equation) self%gradCIC = matmul(self%GradJac, self%residuals) self%gradCIC = self%gradCIC * 2.0_wp/self%refstdev +if (self%verbose) then + call Message%printMessage(' solveHessian_::gradCIC computed') + io_real(1) = minval(self%gradCIC) + io_real(2) = maxval(self%gradCIC) + call Message%WriteValue(' solveHessian_::gradCIC range : ',io_real,2) +end if ! 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) +if (self%verbose) write (*,*) ' solveHessian_::SOL : ', SOL ! 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) +if (self%verbose) write (*,*) ' solveHessian_::normdp : ', normdp end subroutine solveHessian_ diff --git a/Source/EMsoftOOLib/program_mods/mod_HREBSDDIC.f90 b/Source/EMsoftOOLib/program_mods/mod_HREBSDDIC.f90 new file mode 100644 index 0000000..35bbaf8 --- /dev/null +++ b/Source/EMsoftOOLib/program_mods/mod_HREBSDDIC.f90 @@ -0,0 +1,732 @@ +! ################################################################### +! Copyright (c) 2013-2023, 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_HREBSDDIC + !! author: MDG (OO version) + !! version: 1.0 + !! date: 11/07/24 + !! + !! class definition for the EMHREBSDDIC program + !! uses the inverse compositional Gauss–Newton algorithm from ATEX + +use mod_kinds +use mod_global +use mod_ppEBSD, only : ppEBSDNameListType + +IMPLICIT NONE + +! namelist for the EMHREBSDDIC program +type, public :: HREBSDDICNameListType + integer(kind=irg) :: patx + integer(kind=irg) :: paty + integer(kind=irg) :: nbx + integer(kind=irg) :: nby + integer(kind=irg) :: nthreads + real(kind=sgl) :: C11 + real(kind=sgl) :: C12 + real(kind=sgl) :: C44 + real(kind=sgl) :: C13 + real(kind=sgl) :: C33 + character(fnlen) :: patternfile + character(fnlen) :: DIfile + character(fnlen) :: datafile + character(fnlen) :: ppEBSDnml + character(3) :: crystal +end type HREBSDDICNameListType + +! class definition +type, public :: HREBSD_DIC_T +private + character(fnlen) :: nmldeffile = 'EMHREBSDDIC.nml' + type(HREBSDDICNameListType) :: nml + type(ppEBSDNameListType) :: ppEBSDnml + +contains +private + procedure, pass(self) :: readNameList_ + procedure, pass(self) :: writeHDFNameList_ + procedure, pass(self) :: getNameList_ + procedure, pass(self) :: HREBSD_DIC_ + + generic, public :: getNameList => getNameList_ + generic, public :: writeHDFNameList => writeHDFNameList_ + generic, public :: readNameList => readNameList_ + generic, public :: HREBSD_DIC => HREBSD_DIC_ + +end type HREBSD_DIC_T + +! the constructor routine for this class +interface HREBSD_DIC_T + module procedure HREBSDDIC_constructor +end interface HREBSD_DIC_T + +contains + +!-------------------------------------------------------------------------- +type(HREBSD_DIC_T) function HREBSDDIC_constructor( nmlfile ) result(HREBSDDIC) +!DEC$ ATTRIBUTES DLLEXPORT :: HREBSDDIC_constructor +!! author: MDG +!! version: 1.0 +!! date: 11/07/24 +!! +!! constructor for the HREBSD_DIC_T Class; + +use mod_ppEBSD, only : ppEBSD_T + +IMPLICIT NONE + +character(fnlen), INTENT(IN) :: nmlfile +type(ppEBSD_T) :: ppEBSD + +call HREBSDDIC%readNameList(nmlfile) + +! we also need to read in the namelist file from the EMppEBSD program +! this file defines where the original patterns are located as well +! as information about their size etc. +ppEBSD = ppEBSD_T( HREBSDDIC%nml%ppEBSDnml ) +HREBSDDIC%ppEBSDnml = ppEBSD%get_nml() + +end function HREBSDDIC_constructor + +!-------------------------------------------------------------------------- +subroutine HREBSDDIC_destructor(self) +!DEC$ ATTRIBUTES DLLEXPORT :: HREBSDDIC_destructor +!! author: MDG +!! version: 1.0 +!! date: 11/07/24 +!! +!! destructor for the HREBSD_DIC_T Class + +IMPLICIT NONE + +type(HREBSD_DIC_T), INTENT(INOUT) :: self + +call reportDestructor('HREBSD_DIC_T') + +end subroutine HREBSDDIC_destructor + +!-------------------------------------------------------------------------- +subroutine readNameList_(self, nmlfile, initonly) +!DEC$ ATTRIBUTES DLLEXPORT :: readNameList_ +!! author: MDG +!! version: 1.0 +!! date: 11/07/24 +!! +!! read the namelist from an nml file for the HREBSD_DIC_T Class + +use mod_io +use mod_EMsoft + +IMPLICIT NONE + +class(HREBSD_DIC_T), INTENT(INOUT) :: self +character(fnlen),INTENT(IN) :: nmlfile + !! full path to namelist file +logical,OPTIONAL,INTENT(IN) :: initonly + !! fill in the default values only; do not read the file + +type(EMsoft_T) :: EMsoft +type(IO_T) :: Message +logical :: skipread = .FALSE. + +integer(kind=irg) :: patx +integer(kind=irg) :: paty +integer(kind=irg) :: nbx +integer(kind=irg) :: nby +integer(kind=irg) :: nthreads +real(kind=sgl) :: C11 +real(kind=sgl) :: C12 +real(kind=sgl) :: C44 +real(kind=sgl) :: C13 +real(kind=sgl) :: C33 +character(fnlen) :: patternfile +character(fnlen) :: DIfile +character(fnlen) :: datafile +character(fnlen) :: ppEBSDnml +character(3) :: crystal + +namelist / HREBSDDICdata / patx, paty, nthreads, C11, C12, C44, C13, C33, DIfile, & + datafile, patternfile, DIfile, ppEBSDnml, crystal, nbx, nby + +patx = 0 +paty = 0 +nbx = 0 +nby = 0 +nthreads = 1 +C11 = 276.0 +C12 = 159.0 +C44 = 132.0 +C13 = 0.0 +C33 = 0.0 +patternfile = 'undefined' +datafile = 'undefined' +DIfile = 'undefined' +ppEBSDnml = 'undefined' +crystal = 'cub' + +if (present(initonly)) then + if (initonly) skipread = .TRUE. +end if + +if (.not.skipread) then +! read the namelist file + open(UNIT=dataunit,FILE=trim(nmlfile),DELIM='apostrophe',STATUS='old') + read(UNIT=dataunit,NML=HREBSDDICdata) + close(UNIT=dataunit,STATUS='keep') + +! check for required entries + if (trim(patternfile).eq.'undefined') then + call Message%printError('readNameList:',' patternfile name is undefined in '//nmlfile) + end if + + if (trim(DIfile).eq.'undefined') then + call Message%printError('readNameList:',' DIfile file name is undefined in '//nmlfile) + end if + + if (trim(ppEBSDnml).eq.'undefined') then + call Message%printError('readNameList:',' ppEBSDnml file name is undefined in '//nmlfile) + end if + + if (trim(datafile).eq.'undefined') then + call Message%printError('readNameList:',' datafile file name is undefined in '//nmlfile) + end if + +end if + +self%nml%patx = patx +self%nml%paty = paty +self%nml%nbx = nbx +self%nml%nby = nby +self%nml%nthreads = nthreads +self%nml%C11 = C11 +self%nml%C12 = C12 +self%nml%C44 = C44 +self%nml%C13 = C13 +self%nml%C33 = C33 +self%nml%patternfile = patternfile +self%nml%DIfile = DIfile +self%nml%datafile = datafile +self%nml%ppEBSDnml = ppEBSDnml +self%nml%crystal = crystal + +end subroutine readNameList_ + +!-------------------------------------------------------------------------- +function getNameList_(self) result(nml) +!DEC$ ATTRIBUTES DLLEXPORT :: getNameList_ +!! author: MDG +!! version: 1.0 +!! date: 11/07/24 +!! +!! pass the namelist for the HREBSD_DIC_T Class to the calling program + +IMPLICIT NONE + +class(HREBSD_DIC_T), INTENT(INOUT) :: self +type(HREBSDDICNameListType) :: nml + +nml = self%nml + +end function getNameList_ + +!-------------------------------------------------------------------------- +recursive subroutine writeHDFNameList_(self, HDF, HDFnames) +!DEC$ ATTRIBUTES DLLEXPORT :: writeHDFNameList_ +!! author: MDG +!! version: 1.0 +!! date: 11/07/24 +!! +!! write namelist to HDF file + +use mod_HDFsupport +use mod_HDFnames +use stringconstants + +use ISO_C_BINDING + +IMPLICIT NONE + +class(HREBSD_DIC_T), INTENT(INOUT) :: self +type(HDF_T), INTENT(INOUT) :: HDF +type(HDFnames_T), INTENT(INOUT) :: HDFnames + +integer(kind=irg),parameter :: n_int = 5, n_real = 5 +integer(kind=irg) :: hdferr, io_int(n_int) +real(kind=sgl) :: io_real(n_real) +character(20) :: intlist(n_int), reallist(n_real) +character(fnlen) :: dataset, sval(1),groupname +character(fnlen,kind=c_char) :: line2(1),line10(10) + +associate( enl => self%nml ) + +! create the group for this namelist +hdferr = HDF%createGroup(HDFnames%get_NMLlist()) + +! write all the single integers +io_int = (/ enl%nthreads, enl%patx, enl%paty, enl%nbx, enl%nby /) +intlist(1) = 'nthreads' +intlist(2) = 'patx' +intlist(3) = 'paty' +intlist(4) = 'nbx' +intlist(5) = 'nby' +call HDF%writeNMLintegers(io_int, intlist, n_int) + +! write all the single reals +io_real = (/ enl%C11, enl%C12, enl%C44, enl%C13, enl%C33 /) +reallist(1) = 'C11' +reallist(2) = 'C12' +reallist(3) = 'C44' +reallist(4) = 'C13' +reallist(5) = 'C33' +call HDF%writeNMLreals(io_real, reallist, n_real) + +! write all the strings +dataset = 'DIfile' ! SC_DIfile +line2(1) = trim(enl%DIfile) +hdferr = HDF%writeDatasetStringArray(dataset, line2, 1) +if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create DIfile dataset', hdferr) + +dataset = SC_datafile +line2(1) = trim(enl%datafile) +hdferr = HDF%writeDatasetStringArray(dataset, line2, 1) +if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create datafile dataset', hdferr) + +dataset = 'patternfile' ! SC_patternfile +line2(1) = trim(enl%patternfile) +hdferr = HDF%writeDatasetStringArray(dataset, line2, 1) +if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create patternfile dataset', hdferr) + +dataset = 'ppEBSDnml' +line2(1) = trim(enl%ppEBSDnml) +hdferr = HDF%writeDatasetStringArray(dataset, line2, 1) +if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create ppEBSDnml dataset', hdferr) + +dataset = 'crystal' +line2(1) = trim(enl%crystal) +hdferr = HDF%writeDatasetStringArray(dataset, line2, 1) +if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create crystal dataset', hdferr) + +! pop this group off the stack +call HDF%pop() +call HDF%pop() + +end associate + +end subroutine writeHDFNameList_ + +!-------------------------------------------------------------------------- +subroutine HREBSD_DIC_(self, EMsoft, progname, HDFnames) +!DEC$ ATTRIBUTES DLLEXPORT :: HREBSD_DIC_ +!! author: MDG +!! version: 1.0 +!! date: 11/21/24 +!! +!! perform the computations + +use mod_EMsoft +use HDF5 +use mod_HDFsupport +use mod_HDFnames +use mod_OMPsupport +use mod_quaternions +use mod_rotations +use mod_io +use mod_HDFsupport +use mod_math +use HDF5 +use mod_ppEBSD +use mod_DIC +use mod_DIfiles +use mod_timing +use mod_memory +use mod_fftw3 +use stringconstants +use mod_platformsupport +use bspline_module +use bspline_kinds_module, only: wp, ip +use ISO_C_BINDING + +IMPLICIT NONE + +class(HREBSD_DIC_T), INTENT(INOUT) :: self +type(EMsoft_T), INTENT(INOUT) :: EMsoft +character(fnlen), INTENT(INOUT) :: progname +type(HDFnames_T), INTENT(INOUT) :: HDFnames + +type(DIfile_T) :: DIFT +type(QuaternionArray_T) :: qAR +type(IO_T) :: Message +type(Timing_T) :: timer +type(memory_T) :: mem, memth +type(HDF_T) :: HDF +type(e_T) :: eu +type(o_T) :: om +type(q_T) :: qu +type(a_T) :: ax +type(Quaternion_T) :: quat +type(DIC_T) :: DIC +type(HDFnames_T) :: HDFnames2 + +character(fnlen) :: inpfile, HDFstring +real(kind=dbl) :: stepsizes(3), fpar(3), sig, totaltilt, ave, cosang, sinang +real(kind=sgl) :: io_real(6), ss, tstop +real(kind=dbl),allocatable :: PC(:,:), homographies(:,:), normdp(:) +real(kind=sgl),allocatable :: exptpattern(:,:), targetpattern(:,:) +integer(kind=irg),allocatable :: nit(:) +integer(kind=irg) :: hdferr, binx, biny, L, recordsize, patsz, ROI_size, sz(2), i, j, kk, numangles, & + istat, interp_grid, interp_size, info, offset, nbx, nby, ii, jj, ierr, io_int(2), & + 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 +character(11) :: dstr +character(15) :: tstrb +character(15) :: tstre +! integer(HSIZE_T) :: dims2(2), dims3(3), offset3(3) +real(kind=dbl) :: C(6,6) +character(fnlen) :: DIfile, groupname, dataset, datagroupname, outname, fname, hostname +logical :: f_exists, g_exists, overwrite = .TRUE. +integer(kind=4) :: hnstat + +! is this program executed from home or work ? [code will be removed in Release version] +hnstat = system_hostnm(hostname) +write (*,*) trim(hostname) + +! this program reads a dot product EMDI file, including all the experimental parameters +! as well as a pre-processed pattern file generated by EMppEBSD. + +associate( enl => self%nml, dinl=>DIFT%nml ) + +call setRotationPrecision('d') +call openFortranHDFInterface() + +! if (1.eq.0) then +! read all parameters from the EMDI output file +DIfile = EMsoft%generateFilePath('EMdatapathname',trim(enl%DIfile)) + +HDFnames2 = HDFnames_T() +call HDFnames2%set_NMLfiles(SC_NMLfiles) +call HDFnames2%set_NMLfilename(SC_DictionaryIndexingNML) +call HDFnames2%set_NMLparameters(SC_NMLparameters) +call HDFnames2%set_NMLlist(SC_DictionaryIndexingNameListType) + +call DIFT%readDotProductFile(EMsoft, HDF, HDFnames2, DIfile, hdferr, & + getPhi1=.TRUE., & + getPhi=.TRUE., & + getPhi2=.TRUE.) + +fpar = (/ real(dinl%exptnumsx), real(dinl%exptnumsy), real(dinl%delta) /) + +! end if + +! use the memory class for most array allocations +mem = Memory_T() + +call Message%printMessage(' Completed reading DI data from file '//trim(DIfile)) + +! make sure that the patternfile actually exists +f_exists = .FALSE. +fname = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(enl%patternfile) +inquire(file=trim(fname), exist=f_exists) + +if (.not.f_exists) then + call Message%printError('HREBSD_DIC_',' pattern file not found ') +end if + +! define some pattern dimension parameters +binx = dinl%exptnumsx +biny = dinl%exptnumsy +L = binx * biny +recordsize = 4 * L +patsz = L + +write (*,*) ' pattern dimensions : ', binx, biny, L, self%ppEBSDnml%ROI + +! initialize the timing routines +timer = Timing_T() +tstrb = timer%getTimeString() + +! file exists to the first thing we do is extract the reference pattern; this is a simple +! binary file generated by the EMppEBSD program +call mem%alloc(exptpattern, (/ binx,biny /), 'exptpattern', 0.0 ) +open(unit=dataunit,file=trim(fname),& + status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr) +if (sum(self%ppEBSDnml%ROI).ne.0) then + offset = enl%paty * self%ppEBSDnml%ROI(3) + enl%patx + 1 + numpats = (self%ppEBSDnml%ROI(3)-self%ppEBSDnml%ROI(1)+1) * (self%ppEBSDnml%ROI(4)-self%ppEBSDnml%ROI(2)+1) +else + offset = enl%paty * dinl%ipf_wd + enl%patx + 1 + numpats = dinl%ipf_wd * dinl%ipf_ht +end if +write (*,*) ' reference pattern offset = ', offset, recordsize +! offset = 1 +! numpats = 1000 +read(dataunit,rec=offset) exptpattern +close(dataunit,status='keep') + +write (*,*) ' reference entries : ', exptpattern(1:2,1:2) +write (*,*) 'reference pattern : ', minval(exptpattern), maxval(exptpattern) + +! allocate global arrays to store the output of the DIC routine. +call mem%alloc(homographies, (/ 8, numpats /), 'homographies', initval=0.D0 ) +call mem%alloc(normdp, (/ numpats /), 'normdp', initval=0.D0 ) +call mem%alloc(nit, (/ numpats /), 'nit', initval=0 ) + +! open the pattern file so that each thread can read from it +open(unit=dataunit,file=trim(fname),& + status='unknown',form='unformatted',access='direct',recl=recordsize,iostat=ierr) + +call OMP_setNThreads(enl%nthreads) + +memth = memory_T( nt = enl%nthreads ) + +! 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) + +! set the thread ID +TID = OMP_GET_THREAD_NUM() + +!copy the exptpattern into the thread's DIC class +DIC = DIC_T( binx, biny ) +call DIC%setverbose( .FALSE. ) +! call DIC%setverbose( .TRUE. ) +call DIC%setpattern( 'r', dble(exptpattern) ) + +! 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) + +! generate the b-splines for the reference pattern and verify the accuracy +! also get the derivatives of the reference pattern to determine the Hessian +call DIC%getbsplines(refp=.TRUE., verify=.TRUE., grads=.TRUE.) + +! zero-mean and normalize the referenceSR array +call DIC%applyZMN(doreference=.TRUE.) + +! compute the Hessian matrix +call DIC%getHessian(Hessian) + +! 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 + + call DIC%setpattern( 'd', dble(targetpattern) ) + call DIC%getbsplines(defp=.TRUE.) + +! 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,500 + ! write (*,*) '-----------------' + ! write (*,*) ' iteration # ',jj + call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.) + call DIC%applyZMN(dotarget=.TRUE.) + call DIC%getresiduals(CIC) + call DIC%solveHessian(SOL, ndp) + +! convert to a shape function and take the inverse + Wnew = DIC%getShapeFunction(reshape(SOL, (/8/))) + Winv = matrixInvert_wp( Wnew ) + W = matmul( W, Winv ) + ! W = W / W(3,3) + if (ndp.lt.oldnorm) then + oldnorm = ndp + oldW = W + hg = DIC%getHomography(W) + else + W = oldW + ndp = oldnorm + EXIT + end if + end do + + if (mod(ii,25).eq.0) then + io_int(1) = ii + 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) + nit(ii) = jj + call DIC%cleanup() + ! call Message%printMessage(' -------------------------- ') +end do +!$OMP END DO +call memth%dealloc(targetpattern, 'targetpattern', TID=TID) +!$OMP BARRIER +!$OMP END PARALLEL + +close(dataunit,status='keep') + + + +! ! stiffness matrix in the crystal frame +! C = 0.D0 +! C(1,1)=dble(enl%C11) +! C(2,2)=dble(enl%C11) +! C(1,2)=dble(enl%C12) +! C(2,1)=dble(enl%C12) +! C(4,4)=dble(enl%C44) +! C(5,5)=dble(enl%C44) + +! ! add equivalent entries for cubic crystal +! if (enl%crystal.eq.'cub') then +! C(3,3)=dble(enl%C11) +! C(1,3)=dble(enl%C12) +! C(2,3)=dble(enl%C12) +! C(3,1)=dble(enl%C12) +! C(3,2)=dble(enl%C12) +! C(6,6)=dble(enl%C44) +! ! or for hexagonal close packed crystal +! else if (enl%crystal.eq.'hex') then +! C(3,3)=dble(enl%C33) +! C(1,3)=dble(enl%C13) +! C(2,3)=dble(enl%C13) +! C(3,1)=dble(enl%C13) +! C(3,2)=dble(enl%C13) +! C(6,6)=dble(enl%C11-enl%C12)/2.D0 +! else +! call Message%printError('HREBSD_','Undefined crystal type '//trim(enl%crystal)) +! end if + +! call Message%printMessage(' Stiffness tensor (unit:GPa) = ') + +! do i = 1, 6 +! io_real(1:6) = sngl(C(i,1:6)) +! call Message%WriteValue('',io_real,6) +! end do + + + + + + + +! prepare the HDF5 output file +call timer%makeTimeStamp() +dstr = timer%getDateString() +tstre = timer%getTimeString() + +! Create a new file using the default properties. +fname = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(enl%datafile) + +hdferr = HDF%createFile(fname) + +! write the EMheader to the file + datagroupname = trim(HDFnames%get_ProgramData()) + call HDF%writeEMheader(EMsoft,dstr, tstrb, tstre, progname, datagroupname) + +! add the Duration field to the EMheader group + hdferr = HDF%openGroup(HDFnames%get_EMheader()) + hdferr = HDF%openGroup(HDFnames%get_ProgramData()) + +dataset = SC_Duration + call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr) + tstop = 0 + if (g_exists) then + hdferr = HDF%writeDatasetFloat(dataset, tstop, overwrite) + else + hdferr = HDF%writeDatasetFloat(dataset, tstop) + end if + call HDF%pop() + call HDF%pop() + + ! open or create a namelist group to write all the namelist files into +groupname = SC_NMLfiles + hdferr = HDF%createGroup(HDFnames%get_NMLfiles()) + +! read the text file and write the array to the file +dataset = trim(HDFnames%get_NMLfilename()) + hdferr = HDF%writeDatasetTextFile(dataset, EMsoft%nmldeffile) + +! leave this group + call HDF%pop() + +! create a namelist group to write all the namelist files into + hdferr = HDF%createGroup(HDFnames%get_NMLparameters()) + + call self%writeHDFNameList(HDF, HDFnames) + +! then the remainder of the data in a EMData group + hdferr = HDF%createGroup(HDFnames%get_EMData()) + hdferr = HDF%createGroup(HDFnames%get_ProgramData()) + +dataset = 'homographies' + call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr) + if (g_exists) then + hdferr = HDF%writeDatasetDoubleArray(dataset, homographies, 8, numpats, overwrite) + else + hdferr = HDF%writeDatasetDoubleArray(dataset, homographies, 8, numpats) + end if + +dataset = 'ndp' + call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr) + if (g_exists) then + hdferr = HDF%writeDatasetDoubleArray(dataset, normdp, numpats, overwrite) + else + hdferr = HDF%writeDatasetDoubleArray(dataset, normdp, numpats) + end if + +dataset = 'nit' + call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr) + if (g_exists) then + hdferr = HDF%writeDatasetIntegerArray(dataset, nit, numpats, overwrite) + else + hdferr = HDF%writeDatasetIntegerArray(dataset, nit, numpats) + end if + +!call HDF%pop() +!call HDF%popall() + +call closeFortranHDFInterface() + +end associate + +! call mem%allocated_memory_use(' after dealloc') + +end subroutine HREBSD_DIC_ + + +end module mod_HREBSDDIC diff --git a/Source/EMsoftOOLib/program_mods/mod_ppEBSD.f90 b/Source/EMsoftOOLib/program_mods/mod_ppEBSD.f90 new file mode 100644 index 0000000..06df0aa --- /dev/null +++ b/Source/EMsoftOOLib/program_mods/mod_ppEBSD.f90 @@ -0,0 +1,1037 @@ +! ################################################################### +! Copyright (c) 2013-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_ppEBSD + !! author: MDG + !! version: 1.0 + !! date: 12/06/24 + !! + !! class definition for the EMppEBSD program + +use mod_kinds +use mod_global + +IMPLICIT NONE + +! namelist for the EMgetppEBSD program +type, public :: ppEBSDNameListType + integer(kind=irg) :: ipf_ht + integer(kind=irg) :: ipf_wd + integer(kind=irg) :: maskradius + integer(kind=irg) :: numsx + integer(kind=irg) :: numsy + integer(kind=irg) :: nthreads + integer(kind=irg) :: nregions + integer(kind=irg) :: ROI(4) + real(kind=dbl) :: hipassw + character(1) :: maskpattern + character(1) :: filterpattern + character(fnlen) :: exptfile + character(fnlen) :: tmpfile + character(fnlen) :: tiffname + character(fnlen) :: maskfile + character(fnlen) :: inputtype + character(fnlen) :: HDFstrings(10) +end type ppEBSDNameListType + +! class definition +type, public :: ppEBSD_T +private + character(fnlen) :: nmldeffile = 'EMppEBSD.nml' + type(ppEBSDNameListType) :: nml + +contains +private + procedure, pass(self) :: readNameList_ + procedure, pass(self) :: getNameList_ + procedure, pass(self) :: ppEBSD_ + procedure, pass(self) :: get_nml_ + procedure, pass(self) :: get_ipf_ht_ + procedure, pass(self) :: get_ipf_wd_ + procedure, pass(self) :: get_maskradius_ + procedure, pass(self) :: get_numsx_ + procedure, pass(self) :: get_numsy_ + procedure, pass(self) :: get_nthreads_ + procedure, pass(self) :: get_nregions_ + procedure, pass(self) :: get_ROI_ + procedure, pass(self) :: get_hipassw_ + procedure, pass(self) :: get_maskpattern_ + procedure, pass(self) :: get_exptfile_ + procedure, pass(self) :: get_tmpfile_ + procedure, pass(self) :: get_maskfile_ + procedure, pass(self) :: get_inputtype_ + procedure, pass(self) :: get_HDFstrings_ + procedure, pass(self) :: set_ipf_ht_ + procedure, pass(self) :: set_ipf_wd_ + procedure, pass(self) :: set_maskradius_ + procedure, pass(self) :: set_numsx_ + procedure, pass(self) :: set_numsy_ + procedure, pass(self) :: set_nthreads_ + procedure, pass(self) :: set_nregions_ + procedure, pass(self) :: set_ROI_ + procedure, pass(self) :: set_hipassw_ + procedure, pass(self) :: set_maskpattern_ + procedure, pass(self) :: set_exptfile_ + procedure, pass(self) :: set_tmpfile_ + procedure, pass(self) :: set_maskfile_ + procedure, pass(self) :: set_inputtype_ + procedure, pass(self) :: set_HDFstrings_ + + generic, public :: getNameList => getNameList_ + generic, public :: readNameList => readNameList_ + generic, public :: ppEBSD => ppEBSD_ + generic, public :: get_nml => get_nml_ + generic, public :: get_ipf_ht => get_ipf_ht_ + generic, public :: get_ipf_wd => get_ipf_wd_ + generic, public :: get_maskradius => get_maskradius_ + generic, public :: get_numsx => get_numsx_ + generic, public :: get_numsy => get_numsy_ + generic, public :: get_nthreads => get_nthreads_ + generic, public :: get_nregions => get_nregions_ + generic, public :: get_ROI => get_ROI_ + generic, public :: get_hipassw => get_hipassw_ + generic, public :: get_maskpattern => get_maskpattern_ + generic, public :: get_exptfile => get_exptfile_ + generic, public :: get_tmpfile => get_tmpfile_ + generic, public :: get_maskfile => get_maskfile_ + generic, public :: get_inputtype => get_inputtype_ + generic, public :: get_HDFstrings => get_HDFstrings_ + generic, public :: set_ipf_ht => set_ipf_ht_ + generic, public :: set_ipf_wd => set_ipf_wd_ + generic, public :: set_maskradius => set_maskradius_ + generic, public :: set_numsx => set_numsx_ + generic, public :: set_numsy => set_numsy_ + generic, public :: set_nthreads => set_nthreads_ + generic, public :: set_nregions => set_nregions_ + generic, public :: set_ROI => set_ROI_ + generic, public :: set_hipassw => set_hipassw_ + generic, public :: set_maskpattern => set_maskpattern_ + generic, public :: set_exptfile => set_exptfile_ + generic, public :: set_tmpfile => set_tmpfile_ + generic, public :: set_maskfile => set_maskfile_ + generic, public :: set_inputtype => set_inputtype_ + generic, public :: set_HDFstrings => set_HDFstrings_ +end type ppEBSD_T + +! the constructor routine for this class +interface ppEBSD_T + module procedure ppEBSD_constructor +end interface ppEBSD_T + +contains + +!-------------------------------------------------------------------------- +type(ppEBSD_T) function ppEBSD_constructor( nmlfile ) result(ppEBSD) +!DEC$ ATTRIBUTES DLLEXPORT :: ppEBSD_constructor +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! constructor for the ppEBSD_T Class; reads the name list + +IMPLICIT NONE + +character(fnlen), OPTIONAL :: nmlfile + +call ppEBSD%readNameList(nmlfile) + +end function ppEBSD_constructor + +!-------------------------------------------------------------------------- +subroutine ppEBSD_destructor(self) +!DEC$ ATTRIBUTES DLLEXPORT :: ppEBSD_destructor +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! destructor for the ppEBSD_T Class + +IMPLICIT NONE + +type(ppEBSD_T), INTENT(INOUT) :: self + +call reportDestructor('ppEBSD_T') + +end subroutine ppEBSD_destructor + +!-------------------------------------------------------------------------- +function get_nml_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_nml_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get nml from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +type(ppEBSDNameListType) :: out + +out = self%nml + +end function get_nml_ + +!-------------------------------------------------------------------------- +function get_ipf_ht_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_ipf_ht_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get ipf_ht from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%ipf_ht + +end function get_ipf_ht_ + +!-------------------------------------------------------------------------- +subroutine set_ipf_ht_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_ipf_ht_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set ipf_ht in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%ipf_ht = inp + +end subroutine set_ipf_ht_ + +!-------------------------------------------------------------------------- +function get_ipf_wd_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_ipf_wd_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get ipf_wd from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%ipf_wd + +end function get_ipf_wd_ + +!-------------------------------------------------------------------------- +subroutine set_ipf_wd_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_ipf_wd_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set ipf_wd in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%ipf_wd = inp + +end subroutine set_ipf_wd_ + +!-------------------------------------------------------------------------- +function get_maskradius_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_maskradius_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get maskradius from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%maskradius + +end function get_maskradius_ + +!-------------------------------------------------------------------------- +subroutine set_maskradius_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_maskradius_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set maskradius in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%maskradius = inp + +end subroutine set_maskradius_ + +!-------------------------------------------------------------------------- +function get_numsx_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_numsx_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get numsx from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%numsx + +end function get_numsx_ + +!-------------------------------------------------------------------------- +subroutine set_numsx_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_numsx_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set numsx in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%numsx = inp + +end subroutine set_numsx_ + +!-------------------------------------------------------------------------- +function get_numsy_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_numsy_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get numsy from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%numsy + +end function get_numsy_ + +!-------------------------------------------------------------------------- +subroutine set_numsy_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_numsy_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set numsy in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%numsy = inp + +end subroutine set_numsy_ + +!-------------------------------------------------------------------------- +function get_nthreads_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_nthreads_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get nthreads from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%nthreads + +end function get_nthreads_ + +!-------------------------------------------------------------------------- +subroutine set_nthreads_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_nthreads_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set nthreads in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%nthreads = inp + +end subroutine set_nthreads_ + +!-------------------------------------------------------------------------- +function get_nregions_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_nregions_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get nregions from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out + +out = self%nml%nregions + +end function get_nregions_ + +!-------------------------------------------------------------------------- +subroutine set_nregions_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_nregions_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set nregions in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp + +self%nml%nregions = inp + +end subroutine set_nregions_ + +!-------------------------------------------------------------------------- +function get_ROI_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_ROI_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get ROI from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg) :: out(4) + +out = self%nml%ROI + +end function get_ROI_ + +!-------------------------------------------------------------------------- +subroutine set_ROI_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_ROI_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set ROI in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +integer(kind=irg), INTENT(IN) :: inp(4) + +self%nml%ROI = inp + +end subroutine set_ROI_ + +!-------------------------------------------------------------------------- +function get_hipassw_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_hipassw_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get hipassw from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +real(kind=dbl) :: out + +out = self%nml%hipassw + +end function get_hipassw_ + +!-------------------------------------------------------------------------- +subroutine set_hipassw_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_hipassw_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set hipassw in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +real(kind=dbl), INTENT(IN) :: inp + +self%nml%hipassw = inp + +end subroutine set_hipassw_ + +!-------------------------------------------------------------------------- +function get_maskpattern_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_maskpattern_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get maskpattern from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(1) :: out + +out = self%nml%maskpattern + +end function get_maskpattern_ + +!-------------------------------------------------------------------------- +subroutine set_maskpattern_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_maskpattern_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set maskpattern in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(1), INTENT(IN) :: inp + +self%nml%maskpattern = inp + +end subroutine set_maskpattern_ + +!-------------------------------------------------------------------------- +function get_exptfile_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_exptfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get exptfile from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = self%nml%exptfile + +end function get_exptfile_ + +!-------------------------------------------------------------------------- +subroutine set_exptfile_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_exptfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set exptfile in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%exptfile = inp + +end subroutine set_exptfile_ + +!-------------------------------------------------------------------------- +function get_tmpfile_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_tmpfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get tmpfile from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = self%nml%tmpfile + +end function get_tmpfile_ + +!-------------------------------------------------------------------------- +subroutine set_tmpfile_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_tmpfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set tmpfile in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%tmpfile = inp + +end subroutine set_tmpfile_ + +!-------------------------------------------------------------------------- +function get_maskfile_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_maskfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get maskfile from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = self%nml%maskfile + +end function get_maskfile_ + +!-------------------------------------------------------------------------- +subroutine set_maskfile_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_maskfile_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set maskfile in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%maskfile = inp + +end subroutine set_maskfile_ + +!-------------------------------------------------------------------------- +function get_inputtype_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_inputtype_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get inputtype from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen) :: out + +out = self%nml%inputtype + +end function get_inputtype_ + +!-------------------------------------------------------------------------- +subroutine set_inputtype_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_inputtype_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set inputtype in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp + +self%nml%inputtype = inp + +end subroutine set_inputtype_ + +!-------------------------------------------------------------------------- +function get_HDFstrings_(self) result(out) +!DEC$ ATTRIBUTES DLLEXPORT :: get_HDFstrings_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! get HDFstrings from the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen) :: out(10) + +out = self%nml%HDFstrings + +end function get_HDFstrings_ + +!-------------------------------------------------------------------------- +subroutine set_HDFstrings_(self,inp) +!DEC$ ATTRIBUTES DLLEXPORT :: set_HDFstrings_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! set HDFstrings in the ppEBSD_T class + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen), INTENT(IN) :: inp(10) + +self%nml%HDFstrings = inp + +end subroutine set_HDFstrings_ + +!-------------------------------------------------------------------------- +subroutine readNameList_(self, nmlfile, initonly) +!DEC$ ATTRIBUTES DLLEXPORT :: readNameList_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! read the namelist from an nml file for the ppEBSD_T Class + +use mod_io + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +character(fnlen),INTENT(IN) :: nmlfile + !! full path to namelist file +logical,OPTIONAL,INTENT(IN) :: initonly + !! fill in the default values only; do not read the file + +type(IO_T) :: Message +logical :: skipread = .FALSE. + +integer(kind=irg) :: ipf_ht +integer(kind=irg) :: ipf_wd +integer(kind=irg) :: maskradius +integer(kind=irg) :: numsx +integer(kind=irg) :: numsy +integer(kind=irg) :: nthreads +integer(kind=irg) :: nregions +integer(kind=irg) :: ROI(4) +real(kind=dbl) :: hipassw +character(1) :: maskpattern +character(fnlen) :: exptfile +character(fnlen) :: tmpfile +character(fnlen) :: maskfile +character(fnlen) :: inputtype +character(fnlen) :: HDFstrings(10) + +! define the IO namelist to facilitate passing variables to the program. +namelist / ppEBSDdata / numsx, numsy, nregions, maskpattern, nthreads, ipf_ht, ipf_wd, exptfile, maskradius, inputtype, & + tmpfile, maskfile, HDFstrings, hipassw, ROI + +! set the input parameters to default values + ipf_ht = 100 + ipf_wd = 100 + maskfile = 'undefined' + maskpattern = 'n' + maskradius = 240 + hipassw = 0.05 + nregions = 10 + numsx = 0 + numsy = 0 + ROI = (/ 0, 0, 0, 0 /) + exptfile = 'undefined' + inputtype = 'Binary' + HDFstrings = '' + tmpfile = 'EMEBSDDict_tmp.data' + nthreads = 1 + +if (present(initonly)) then + if (initonly) skipread = .TRUE. +end if + +if (.not.skipread) then +! read the namelist file + open(UNIT=dataunit,FILE=trim(nmlfile),DELIM='apostrophe',STATUS='old') + read(UNIT=dataunit,NML=ppEBSDdata) + close(UNIT=dataunit,STATUS='keep') + +! check for required entries + if (trim(exptfile).eq.'undefined') then + call Message%printError('readNameList:',' experimental file name is undefined in '//nmlfile) + end if + + if (numsx.eq.0) then + call Message%printError('readNameList:',' patterns size numsx is zero in '//nmlfile) + end if + + if (numsy.eq.0) then + call Message%printError('readNameList:',' patterns size numsy is zero in '//nmlfile) + end if + end if + +! if we get here, then all appears to be ok, and we need to fill in the enl fields +self%nml%ipf_ht = ipf_ht +self%nml%ipf_wd = ipf_wd +self%nml%maskradius = maskradius +self%nml%numsx = numsx +self%nml%numsy = numsy +self%nml%nthreads = nthreads +self%nml%nregions = nregions +self%nml%ROI = ROI +self%nml%hipassw = hipassw +self%nml%maskpattern = maskpattern +self%nml%exptfile = exptfile +self%nml%tmpfile = tmpfile +self%nml%maskfile = maskfile +self%nml%inputtype = inputtype +self%nml%HDFstrings = HDFstrings + +end subroutine readNameList_ + +!-------------------------------------------------------------------------- +function getNameList_(self) result(nml) +!DEC$ ATTRIBUTES DLLEXPORT :: getNameList_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! pass the namelist for the ppEBSD_T Class to the calling program + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +type(ppEBSDNameListType) :: nml + +nml = self%nml + +end function getNameList_ + +!-------------------------------------------------------------------------- +subroutine ppEBSD_(self, EMsoft, progname) +!DEC$ ATTRIBUTES DLLEXPORT :: ppEBSD_ +!! author: MDG +!! version: 1.0 +!! date: 12/06/24 +!! +!! perform the computations + +use mod_EMsoft +use mod_patterns +use mod_filters +use mod_io +use HDF5 +use mod_HDFsupport +use mod_DIfiles +use mod_memory + +use, intrinsic :: iso_fortran_env + +IMPLICIT NONE + +class(ppEBSD_T), INTENT(INOUT) :: self +type(EMsoft_T), INTENT(INOUT) :: EMsoft +character(fnlen), INTENT(INOUT) :: progname + +type(IO_T) :: Message +type(HDF_T) :: HDF +type(memory_T) :: mem + +integer(kind=irg) :: num,ierr,irec,istat +integer(kind=irg) :: L,totnumexpt,imght,imgwd,nnk, recordsize, iii, hdferr,& + recordsize_correct, patsz +real(kind=sgl),allocatable :: mask(:,:),masklin(:) + +integer(kind=irg) :: i,j,ii,jj,kk,ll,mm,pp,qq, io_int(3), iiistart, iiiend, jjend +real(kind=sgl) :: mi, ma, io_real(2), tstart, tmp, vlen, tstop +integer(kind=irg) :: binx,biny,TID,nthreads +integer(kind=irg) :: correctsize +logical :: f_exists, ROIselected +character(1000) :: charline +type(EBSDDINameListType) :: dinl +character(fnlen) :: fname + +call openFortranHDFInterface() + +associate(ppnl=>self%nml) + +! memory class +mem = memory_T() + +! copy various constants from the namelist +L = ppnl%numsx*ppnl%numsy +totnumexpt = ppnl%ipf_wd*ppnl%ipf_ht +imght = ppnl%numsx +imgwd = ppnl%numsy +recordsize = L*4 +binx = ppnl%numsx +biny = ppnl%numsy + +! make sure that correctsize is a multiple of 16; if not, make it so +! (probably not necessary in this program but kep in place for consistency) +if (mod(L,16) .ne. 0) then + correctsize = 16*ceiling(float(L)/16.0) +else + correctsize = L +end if + +! determine the experimental and dictionary sizes in bytes +recordsize_correct = correctsize*4 +patsz = correctsize + +ROIselected = .FALSE. +if (sum(ppnl%ROI).ne.0) ROIselected = .TRUE. + +!========================================= +! ALLOCATION AND INITIALIZATION OF ARRAYS +!========================================= +call mem%alloc(mask, (/ binx,biny /), 'mask', 1.0) +call mem%alloc(masklin, (/ L /), 'masklin', 0.0) + +!===================================================== +! define the circular mask if necessary and convert to 1D vector +!===================================================== + +if (trim(ppnl%maskfile).ne.'undefined') then +! read the mask from file; the mask can be defined by a 2D array of 0 and 1 values +! that is stored in row form as strings, e.g. +! 0000001110000000 +! 0000011111000000 +! ... etc +! + f_exists = .FALSE. + fname = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(ppnl%maskfile) + inquire(file=trim(fname), exist=f_exists) + if (f_exists.eqv..TRUE.) then + mask = 0.0 + open(unit=dataunit,file=trim(fname),status='old',form='formatted') + do jj=biny,1,-1 + read(dataunit,"(A)") charline + do ii=1,binx + if (charline(ii:ii).eq.'1') mask(ii,jj) = 1.0 + end do + end do + close(unit=dataunit,status='keep') + else + call Message%printError('ppEBSD','maskfile '//trim(fname)//' does not exist') + end if +else + if (ppnl%maskpattern.eq.'y') then + do ii = 1,biny + do jj = 1,binx + if((ii-biny/2)**2 + (jj-binx/2)**2 .ge. ppnl%maskradius**2) then + mask(jj,ii) = 0.0 + end if + end do + end do + end if +end if + +! convert the mask to a linear (1D) array +do ii = 1,biny + do jj = 1,binx + masklin((ii-1)*binx+jj) = mask(jj,ii) + end do +end do + + +!===================================================== +! Preprocess all the experimental patterns and store +! them in a temporary file as vectors +!===================================================== + +! first, make sure that this file does not already exist +f_exists = .FALSE. +fname = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(ppnl%tmpfile) +inquire(file=trim(fname), exist=f_exists) + +! if the file exists, delete and recreate it +if (f_exists) then + open(unit=dataunit,file=trim(fname),& + status='unknown',form='unformatted',access='direct',recl=recordsize_correct,iostat=ierr) + close(unit=dataunit,status='delete') +end if + +! copy the relevant ppnl parameters into the dinl structure +dinl%nthreads = ppnl%nthreads +dinl%hipassw = ppnl%hipassw +dinl%ROI = ppnl%ROI +dinl%ipf_wd = ppnl%ipf_wd +dinl%ipf_ht = ppnl%ipf_ht +dinl%tmpfile = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(ppnl%tmpfile) +dinl%exptfile = trim(EMsoft%generateFilePath('EMdatapathname'))//trim(ppnl%exptfile) +dinl%inputtype = ppnl%inputtype +dinl%HDFstrings = ppnl%HDFstrings +dinl%nregions = ppnl%nregions +dinl%DIModality = 'EBSD' + +call PreProcessPatterns(EMsoft, HDF, .FALSE., dinl, binx, biny, masklin, correctsize, totnumexpt) + +!===================================================== +call Message%printMessage(' -> created pre-processed EBSD pattern file '//trim(fname)) + +! close the fortran HDF5 interface +call closeFortranHDFInterface() + +end associate + +call mem%dealloc(mask, 'mask') +call mem%dealloc(masklin, 'masklin') + +! call mem%allocated_memory_use(' from end of ppEBSD_ subroutine ') + +end subroutine ppEBSD_ + + + +end module mod_ppEBSD diff --git a/Source/EMsoftOOLib/stringconstants.in b/Source/EMsoftOOLib/stringconstants.in index cc505dc..2b675b4 100644 --- a/Source/EMsoftOOLib/stringconstants.in +++ b/Source/EMsoftOOLib/stringconstants.in @@ -315,6 +315,9 @@ OrderParameter;"CPLP" ppEBSD;"ppEBSD" ppEBSDNameList;"ppEBSDNameList" ppEBSDNML;"ppEBSDNML" +HREBSDDIC;"HREBSDDIC" +HREBSDDICNameList;"HREBSDDICNameList" +HREBSDDICNML;"HREBSDDICNML" PEDZA;"PEDZA" PEDZANML;"PEDZANML" PEDZANameList;"PEDZANameList" diff --git a/Source/SEM/CMakeLists.txt b/Source/SEM/CMakeLists.txt index 6b7b487..75e18cd 100644 --- a/Source/SEM/CMakeLists.txt +++ b/Source/SEM/CMakeLists.txt @@ -287,6 +287,15 @@ if(EMsoftOO_ENABLE_HDF5_SUPPORT) INCLUDE_DIRS ${EMsoftOOLib_BINARY_DIR} ) + Add_EMsoftOO_Executable(TARGET EMHREBSDDIC + SOURCES ${APP_DIR}/EMHREBSDDIC.f90 + LINK_LIBRARIES ${EXE_LINK_LIBRARIES} + TEMPLATE ${TMPLT_DIR}/EMHREBSDDIC.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} diff --git a/Source/SEM/EMHREBSDDIC.f90 b/Source/SEM/EMHREBSDDIC.f90 new file mode 100644 index 0000000..cdb3907 --- /dev/null +++ b/Source/SEM/EMHREBSDDIC.f90 @@ -0,0 +1,67 @@ +! ################################################################### +! 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 EMHREBSDDIC + !! author: MDG + !! version: 1.0 + !! date: 06/28/23 + !! + !! f90 version of the DIC implementation by Ernould et al. + +use mod_kinds +use mod_global +use mod_EMsoft +use mod_HREBSDDIC +use mod_HDFnames +use stringconstants + +IMPLICIT NONE + +character(fnlen) :: progname = 'EMHREBSDDIC.f90' +character(fnlen) :: progdesc = 'Implementation of the HR-EBSD-DIC algorithm' + +type(EMsoft_T) :: EMsoft +type(HREBSD_DIC_T) :: HREBSDDIC +type(HDFnames_T) :: HDFnames + +! print the EMsoft header and handle any command line arguments +EMsoft = EMsoft_T( progname, progdesc, tpl = (/ 289 /) ) + +! deal with the namelist stuff +HREBSDDIC = HREBSD_DIC_T(EMsoft%nmldeffile) + +HDFnames = HDFnames_T() +call HDFnames%set_ProgramData(SC_HREBSDDIC) +call HDFnames%set_NMLlist(SC_HREBSDDICNameList) +call HDFnames%set_NMLfilename(SC_HREBSDDICNML) +call HDFnames%set_Variable(SC_HREBSDDIC) + +! perform the computations +call HREBSDDIC%HREBSD_DIC(EMsoft, progname, HDFnames) + +end program EMHREBSDDIC diff --git a/Source/TestPrograms/play.f90 b/Source/TestPrograms/play.f90 index d72e3bf..fc06537 100644 --- a/Source/TestPrograms/play.f90 +++ b/Source/TestPrograms/play.f90 @@ -47,6 +47,7 @@ program EMplay use mod_DIC use HDF5 use mod_HDFsupport +use mod_platformsupport use bspline_module use bspline_kinds_module, only: wp, ip @@ -60,6 +61,7 @@ program EMplay type(q_T) :: qu type(o_T) :: om type(e_T) :: eu +type(IO_T) :: Message real(kind=wp), allocatable :: refpat(:,:), defpat(:,:) real(kind=sgl),allocatable :: tmppat(:,:) @@ -67,9 +69,12 @@ program EMplay 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, jj +integer(kind=irg) :: nx, ny, nxy, nbx, nby, i, ii, j, NSR, cnt, nxSR, nySR, jj, recordsize, ierr real(wp) :: tol +integer(kind=4) :: hnstat +character(fnlen) :: fname, gname, hostname +hnstat = system_hostnm(hostname) EMsoft = EMsoft_T( progname, progdesc ) @@ -86,14 +91,43 @@ program EMplay ! instantiate the DIC class ! this also initializes the x and y coordinate arrays DIC = DIC_T( nx, ny ) +call DIC%setverbose(.FALSE.) +! call DIC%setverbose(.TRUE.) nxy = nx * ny +recordsize = nxy * 4 -allocate(refpat(0:nx-1,0:ny-1), defpat(0:nx-1,0:ny-1)) +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/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') 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') +! 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/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) + + ! read(20) refpat +! read(dataunit, rec=1) 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) + call DIC%setpattern('r', refpat) @@ -108,47 +142,61 @@ program EMplay nby = 30 call DIC%defineSR( nbx, nby, 0.5_wp, 0.5_wp) +call DIC%getbsplines(refp=.TRUE., verify=.TRUE., grads=.TRUE.) -call DIC%getbsplines(refp=.TRUE., verify=.TRUE.) +!zero-mean and normalize the referenceSR and targetSR arrays +call DIC%applyZMN(doreference=.TRUE.) + +! compute the Hessian matrix +call DIC%getHessian( Hessian ) ! 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') +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') -open(unit=dataunit,file='/Volumes/Drive2/playarea/DIC/test/homography-fits2.data',status='unknown',form='unformatted') +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 +do jj=1, 1000 + 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 /) -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) + 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) + +write (dataunit,rec=jj) DIC%getpattern('d',nx,ny) Woriginal = DIC%getShapeFunction(horiginal) ! get the derivatives of the reference pattern to determine the Hessian -call DIC%getbsplines(defp=.TRUE., grads=.TRUE.) +call DIC%getbsplines(defp=.TRUE.) -! zero-mean and normalize the referenceSR and targetSR arrays -call DIC%applyZMN(doreference=.TRUE., dotarget=.TRUE.) +! ! zero-mean and normalize the referenceSR and targetSR arrays +! call DIC%applyZMN(dotarget=.TRUE.) -! initial CIC function (to be minimized ...) -call DIC%getresiduals( CIC ) -write (*,*) ' initial CIC value : ', CIC - -! compute the Hessian matrix -call DIC%getHessian( Hessian ) +! ! initial CIC function (to be minimized ...) +! call DIC%getresiduals( CIC ) +! stop ! initialize the homography coefficients to zero ! in a more general implementation we might initialize them to the @@ -159,7 +207,7 @@ program EMplay oldnorm = 100.0_wp ! and here we start the loop -do ii=1,100 +do ii=1,500 write (*,*) ' iteration # ',ii call DIC%applyHomography(hg, 0.5_wp, 0.5_wp, dotarget=.TRUE.) @@ -178,38 +226,41 @@ program EMplay W = matmul( W, Winv ) ! W = W / W(3,3) hg = DIC%getHomography(W) - write (*,*) hg - write (*,*) ' norm(deltap) = ', normdp - write (*,*) '------------' - if (normdp.lt.1.1*oldnorm) then + ! write (*,*) hg + ! write (*,*) ' norm(deltap) = ', normdp + ! write (*,*) '------------' + if (normdp.lt.oldnorm) then oldnorm = normdp oldW = W else W = oldW EXIT end if + ! call Message%printMessage(' ---------------------- ') 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 (*,*) '' -write (*,*) ' Differences : ' -write (*,*) horiginal-hg +! write (*,*) '' +! write (*,*) ' Final homography : ' +! write (*,*) DIC%getHomography(W) +! write (*,*) '' +! write (*,*) ' Target homography : ' +! write (*,*) horiginal +! 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 +write (unit=28,FMT='(9(F12.8,","),I4)') real(hg), real(normdp), ii + +! if (jj.eq.3) stop + call DIC%cleanup() end do ! loop over jj +close(28,status='keep') close(dataunit,status='keep') end program EMplay