adds a normalized cross correlation as an alternative similarity metric

NamelistTemplates/EMEBSDDI.template

@@ -9,6 +9,8 @@
 !
 ! 'dynamic' for on the fly indexing or 'static' for pre calculated dictionary
  indexingmode = 'dynamic',
+! use the normalized dot product 'ndp' or the normalized cross correlation 'ncc' as the similarity metric
+ similaritymetric = 'ndp',
 !
 !###################################################################
 ! DICTIONARY PARAMETERS: COMMON TO 'STATIC' AND 'DYNAMIC'

Source/EMOpenCLLib/Indexingmod.f90

@@ -85,6 +85,7 @@ module Indexingmod
 !> @date 02/07/20 MDG 4.0 integration with new wrapper routine; adds optional arguments 
 !> @date 02/22/20 MDG 4.1 split of call back routines
 !> @date 06/01/20 MDG 5.0 change in handling of pattern binning (as of EMsoft version 5.0.3)
+!> @date 10/06/20 MDG 5.1 added normalized cross correlation as optional similarity metric
 !--------------------------------------------------------------------------
 subroutine EBSDDIdriver(Cnmldeffile, Cprogname, cproc, ctimeproc, cerrorproc, objAddress, cancel) &
            bind(c, name='EBSDDIdriver') 
@@ -266,8 +267,8 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
 real(kind=dbl)                                      :: prefactor
 character(fnlen)                                    :: xtalname
 integer(kind=irg)                                   :: binx,biny,TID,nthreads,Emin,Emax, iiistart, iiiend, jjend
-real(kind=sgl)                                      :: sx,dx,dxm,dy,dym,rhos,x,projweight, dp, mvres, nel, emult
-real(kind=sgl)                                      :: dc(3),quat(4),ixy(2),bindx
+real(kind=sgl)                                      :: sx,dx,dxm,dy,dym,rhos,x,projweight, dp, mvres, nel, emult, mean, sdev
+real(kind=sgl)                                      :: dc(3),quat(4),ixy(2),bindx, Nval, Nval2
 integer(kind=irg)                                   :: nix,niy,nixp,niyp
 real(kind=sgl)                                      :: euler(3)
 integer(kind=irg)                                   :: indx
@@ -302,6 +303,13 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
 ! deal with the namelist stuff
 call GetEBSDIndexingNameList(nmldeffile,dinl)
 
+if (dinl%similaritymetric.eq.'ndp') then 
+  call Message('--> Using normalized dot product as similarity metric')
+end if 
+if (dinl%similaritymetric.eq.'ncc') then 
+  call Message('--> Using normalized cross correlation as similarity metric')
+end if 
+
 ! binned pattern array size for experimental patterns
 binx = dinl%exptnumsx/dinl%binning
 biny = dinl%exptnumsy/dinl%binning
@@ -953,6 +961,10 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
   cancelled = .FALSE.
 end if 
 
+Nval = 1.0/float(binx*biny)
+Nval2 = 1.0/float(binx*biny-1)
+
+
 dictionaryloop: do ii = 1,cratio+1
     results = 0.0
 
@@ -976,7 +988,8 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
     end if
 
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID,iii,jj,ll,mm,pp,ierr,io_int, vlen, tock, ttime, dicttranspose, EBSDdictpatflt) &
-!$OMP& PRIVATE(binned, ma, mi, EBSDpatternintd, EBSDpatterninteger, EBSDpatternad, quat, imagedictflt, imagedictfltflip)
+!$OMP& PRIVATE(binned, ma, mi, EBSDpatternintd, EBSDpatterninteger, EBSDpatternad, quat, imagedictflt, imagedictfltflip) &
+!$OMP& PRIVATE(mean, sdev)
 
 ! allocate the local arrays that are used by each thread
     allocate(EBSDpatterninteger(binx,biny))
@@ -1032,8 +1045,11 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
 
         call InnerProdGPU(cl_expt,cl_dict,Ne,Nd,correctsize,results,numd,dinl%devid,kernel,context,command_queue)
 
+        if (dinl%similaritymetric.eq.'ncc') results = results * Nval
+
         dp =  maxval(results)
         if (dp.gt.mvres) mvres = dp
+        write (*,*) jj,pp,'; max/min ncc = ',dp, minval(results) 
 
 ! this might be simplified later for the remainder of the patterns
         do qq = 1,ppendE(jj)
@@ -1139,23 +1155,33 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
          call CalcEBSDPatternSingleFull(jpar,quat,accum_e_MC,mLPNH,mLPSH,EBSDdetector%rgx,&
                                         EBSDdetector%rgy,EBSDdetector%rgz,binned,Emin,Emax,mask,prefactor)
 
+
          if (dinl%scalingmode .eq. 'gam') then
            binned = binned**dinl%gammavalue
          end if
 
-! hi pass filtering
+! hi pass filtering (disabled)
 !      rdata = dble(binned)
 !      fdata = HiPassFilter(rdata,dims,w)
 !      binned = sngl(fdata)
 
+
+         if (dinl%similaritymetric.eq.'ndp') then 
 ! adaptive histogram equalization
-         ma = maxval(binned)
-         mi = minval(binned)
-
-         EBSDpatternintd = ((binned - mi)/ (ma-mi))
-         EBSDpatterninteger = nint(EBSDpatternintd*255.0)
-         EBSDpatternad =  adhisteq(dinl%nregions,binx,biny,EBSDpatterninteger)
-         binned = float(EBSDpatternad)
+           ma = maxval(binned)
+           mi = minval(binned)
+
+           EBSDpatternintd = ((binned - mi)/ (ma-mi))
+           EBSDpatterninteger = nint(EBSDpatternintd*255.0)
+           EBSDpatternad =  adhisteq(dinl%nregions,binx,biny,EBSDpatterninteger)
+           binned = float(EBSDpatternad)
+         else  ! use normalized cross correlation 
+           binned = binned * mask
+           mean = sum(binned) * Nval
+           binned = binned - mean 
+           sdev = sqrt(Nval2 * sum( binned*binned ))
+           binned = binned / sdev 
+         end if 
 
          imagedictflt = 0.0
          imagedictfltflip = 0.0
@@ -1166,13 +1192,15 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
          end do
 
 ! normalize and apply circular mask 
-         imagedictflt(1:L) = imagedictflt(1:L) * masklin(1:L)
-         vlen = vecnorm(imagedictflt(1:correctsize))
-         if (vlen.ne.0.0) then
-           imagedictflt(1:correctsize) = imagedictflt(1:correctsize)/vlen
-         else
-           imagedictflt(1:correctsize) = 0.0
-         end if
+         if (dinl%similaritymetric.eq.'ndp') then 
+           imagedictflt(1:L) = imagedictflt(1:L) * masklin(1:L)
+           vlen = vecnorm(imagedictflt(1:correctsize))
+           if (vlen.ne.0.0) then
+             imagedictflt(1:correctsize) = imagedictflt(1:correctsize)/vlen
+           else
+             imagedictflt(1:correctsize) = 0.0
+           end if
+         end if 
 
          dict((pp-1)*correctsize+1:pp*correctsize) = imagedictflt(1:correctsize)
 

Source/EMsoftHDFLib/patternmod.f90

@@ -1470,6 +1470,7 @@ end subroutine patternmod_errormessage
 !> @date 04/01/18 MDG 1.0 original
 !> @date 04/04/18 MDG 1.1 added optional exptIQ parameter
 !> @date 06/02/20 MDG 2.0 changed handling of pattern binning (version 5.0.3)
+!> @date 10/06/20 MDG 2.1 add option for normalized cross correlation
 !--------------------------------------------------------------------------
 recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, masklin, correctsize, totnumexpt, &
                                         epatterns, exptIQ)
@@ -1515,11 +1516,11 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
 integer(HSIZE_T)                            :: dims3(3), offset3(3)
 integer(kind=irg),parameter                 :: iunitexpt = 41, itmpexpt = 42, itmpexpt2 = 43
 integer(kind=irg)                           :: tickstart, tstop 
-real(kind=sgl)                              :: vlen, tmp, ma, mi, io_real(1)
+real(kind=sgl)                              :: vlen, tmp, ma, mi, io_real(1), Nval, Nval2, mean, sdev
 real(kind=dbl)                              :: w, Jres, sclfct
 integer(kind=irg),allocatable               :: EBSDpint(:,:)
 real(kind=sgl),allocatable                  :: tmpimageexpt(:), EBSDPattern(:,:), exppatarray(:), EBSDpat(:,:), pat1D(:), &
-                                               Pepatterns(:,:)
+                                               Pepatterns(:,:), mask(:,:)
 real(kind=dbl),allocatable                  :: rrdata(:,:), ffdata(:,:), ksqarray(:,:), inpat(:,:), outpat(:,:)
 complex(kind=dbl),allocatable               :: hpmask(:,:)
 complex(C_DOUBLE_COMPLEX),pointer           :: inp(:,:), outp(:,:)
@@ -1529,6 +1530,11 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
 
 call Message('Preprocessing experimental patterns')
 
+if (ebsdnl%similaritymetric.eq.'ncc') then
+  allocate(mask(ebsdnl%exptnumsx,ebsdnl%exptnumsy))
+  mask = reshape(masklin,(/ebsdnl%exptnumsx,ebsdnl%exptnumsy/))
+end if 
+
 !===================================================================================
 ! define a bunch of mostly integer parameters
 
@@ -1685,13 +1691,16 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
   call WriteValue(' binned size  : ',io_int,2,"(I4,' x ',I4)")
 end if 
 
+Nval = 1.0/float(Px*Py)
+Nval2 = 1.0/float(Px*Py-1)
+
 ! ===================================================================================
 ! we do one row at a time
 prepexperimentalloop: do iii = iiistart,iiiend
 
 ! start the OpenMP portion
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, jj, kk, mi, ma, istat, ip, op) &
-!$OMP& PRIVATE(tmpimageexpt, EBSDPat, rrdata, ffdata, EBSDpint, vlen, tmp, inp, outp)
+!$OMP& PRIVATE(tmpimageexpt, EBSDPat, rrdata, ffdata, EBSDpint, vlen, tmp, inp, outp, mean, sdev)
 
 ! set the thread ID
     TID = OMP_GET_THREAD_NUM()
@@ -1736,6 +1745,8 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
 ! then loop in parallel over all patterns to perform the preprocessing steps
 !$OMP DO SCHEDULE(DYNAMIC)
     do jj=1,jjend
+
+
 ! convert imageexpt to 2D EBSD Pattern array
         do kk=1,Py
           EBSDPat(1:Px,kk) = exppatarray((jj-1)*Ppatsz+(kk-1)*Px+1:(jj-1)*Ppatsz+kk*Px)
@@ -1751,25 +1762,37 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
         ffdata = applyHiPassFilter(rrdata, (/ Px, Py /), w, hpmask, inp, outp, HPplanf, HPplanb)
         EBSDPat = sngl(ffdata)
 
+! in ncc mode, we multiply by the mask before scaling the intensities
+
+        if (ebsdnl%similaritymetric.eq.'ndp') then 
 ! adaptive histogram equalization
-        ma = maxval(EBSDPat)
-        mi = minval(EBSDPat)
+          ma = maxval(EBSDPat)
+          mi = minval(EBSDPat)
 
-        EBSDpint = nint(((EBSDPat - mi) / (ma-mi))*255.0)
-        EBSDPat = float(adhisteq(ebsdnl%nregions,Px,Py,EBSDpint))
+          EBSDpint = nint(((EBSDPat - mi) / (ma-mi))*255.0)
+          EBSDPat = float(adhisteq(ebsdnl%nregions,Px,Py,EBSDpint))
+        else  ! use normalized cross correlation so subtract mean and divided by standard deviation
+          EBSDPat = EBSDPat * mask
+          mean = sum(EBSDPat) * Nval
+          EBSDPat = EBSDPat - mean 
+          sdev = sqrt(Nval2 * sum( EBSDPat*EBSDPat ))
+          EBSDPat = EBSDPat / sdev
+        end if
 
 ! convert back to 1D vector
         do kk=1,Py
           exppatarray((jj-1)*Ppatsz+(kk-1)*Px+1:(jj-1)*Ppatsz+kk*Px) = EBSDPat(1:Px,kk)
         end do
 
 ! apply circular mask and normalize for the dot product computation
-        exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) * masklin(1:PL)
-        vlen = vecnorm(exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL))
-        if (vlen.ne.0.0) then
-          exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL)/vlen
-        else
-          exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = 0.0
+        if (ebsdnl%similaritymetric.eq.'ndp') then 
+          exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) * masklin(1:PL)
+          vlen = vecnorm(exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL))
+          if (vlen.ne.0.0) then
+            exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL)/vlen
+          else
+            exppatarray((jj-1)*Ppatsz+1:(jj-1)*Ppatsz+PL) = 0.0
+          end if
         end if
     end do
 !$OMP END DO
@@ -1857,7 +1880,7 @@ recursive subroutine PreProcessPatterns(nthreads, inRAM, ebsdnl, binx, biny, mas
 ! convert outpat to a 1D array for writing to file or epatterns array
     pat1D = sngl(reshape(outpat, (/ BL /) ))
 ! normalize the pattern 
-    pat1D = pat1D/vecnorm(pat1D)
+    if (ebsdnl%similaritymetric.eq.'ndp') pat1D = pat1D/vecnorm(pat1D)
     if (inRAM.eqv..TRUE.) then
       epatterns(1:Bpatsz, jj) = pat1D(1:BL)
     else 

Source/EMsoftLib/NameListHandlers.f90

@@ -4107,10 +4107,11 @@ recursive subroutine GetBSENameList(nmlfile, enl, initonly)
 character(fnlen)        :: useangles
 character(fnlen)        :: imagefile
 character(fnlen)        :: masterfile
+character(fnlen)        :: Kosselmasterfile
 character(fnlen)        :: datafile
 
 ! define the IO namelist to facilitate passing variables to the program.
-namelist  / BSEdata / energymin, energymax, beamcurrent, dwelltime, gammavalue, workingdistance, BSEdistance, &
+namelist  / BSEdata / energymin, energymax, beamcurrent, dwelltime, gammavalue, workingdistance, BSEdistance, Kosselmasterfile, &
                       rin, rout, NsqL, nthreads, scalingmode, useangles, imagefile, masterfile, datafile, incidence
 
 
@@ -4131,6 +4132,7 @@ recursive subroutine GetBSENameList(nmlfile, enl, initonly)
 useangles = 'original'
 imagefile = 'undefined'
 masterfile = 'undefined'
+Kosselmasterfile = 'undefined'
 datafile = 'undefined'
 
 if (present(initonly)) then
@@ -4159,6 +4161,10 @@ recursive subroutine GetBSENameList(nmlfile, enl, initonly)
   call FatalError('GetBSENameList:',' master pattern file name is undefined in '//nmlfile)
  end if
 
+ if (trim(Kosselmasterfile).eq.'undefined') then
+  call FatalError('GetBSENameList:',' Kossel master pattern file name is undefined in '//nmlfile)
+ end if
+
  if (trim(datafile).eq.'undefined') then
   call FatalError('GetBSENameList:',' DI input file name is undefined in '//nmlfile)
  end if
@@ -4182,6 +4188,7 @@ recursive subroutine GetBSENameList(nmlfile, enl, initonly)
 enl%useangles = useangles
 enl%imagefile = trim(imagefile)
 enl%masterfile = trim(masterfile)
+enl%Kosselmasterfile = trim(Kosselmasterfile)
 enl%datafile = trim(datafile)
 
 end subroutine GetBSENameList
@@ -6460,6 +6467,7 @@ recursive subroutine GetEBSDIndexingNameList(nmlfile, enl, initonly)
 character(1)                                      :: maskpattern
 character(1)                                      :: keeptmpfile
 character(3)                                      :: scalingmode
+character(3)                                      :: similaritymetric
 character(3)                                      :: Notify
 character(fnlen)                                  :: dotproductfile
 character(fnlen)                                  :: masterfile
@@ -6497,7 +6505,7 @@ recursive subroutine GetEBSDIndexingNameList(nmlfile, enl, initonly)
 scalingmode, maskpattern, energyaverage, L, omega, nthreads, energymax, datafile, angfile, ctffile, &
 ncubochoric, numexptsingle, numdictsingle, ipf_ht, ipf_wd, nnk, nnav, exptfile, maskradius, inputtype, &
 dictfile, indexingmode, hipassw, stepX, stepY, tmpfile, avctffile, nosm, eulerfile, Notify, maskfile, &
-section, HDFstrings, ROI, keeptmpfile, multidevid, usenumd, nism, isangle, refinementNMLfile
+section, HDFstrings, ROI, keeptmpfile, multidevid, usenumd, nism, isangle, refinementNMLfile, similaritymetric
 
 ! set the input parameters to default values (except for xtalname, which must be present)
 ncubochoric     = 50
@@ -6539,6 +6547,7 @@ recursive subroutine GetEBSDIndexingNameList(nmlfile, enl, initonly)
 maskpattern     = 'n'           ! 'y' or 'n' to include a circular mask
 Notify          = 'Off'
 scalingmode     = 'not'         ! intensity selector ('lin', 'gam', or 'not')
+similaritymetric = 'ndp'
 masterfile      = 'undefined'   ! filename
 dotproductfile  = 'undefined'
 energymin       = 10.0
@@ -6664,6 +6673,7 @@ recursive subroutine GetEBSDIndexingNameList(nmlfile, enl, initonly)
 enl%beamcurrent   = beamcurrent
 enl%dwelltime     = dwelltime
 enl%scalingmode   = scalingmode
+enl%similaritymetric = similaritymetric
 enl%ncubochoric   = ncubochoric
 enl%omega         = omega
 enl%energymin     = energymin

Source/EMsoftLib/NameListTypedefs.f90

@@ -630,6 +630,7 @@ module NameListTypedefs
         character(fnlen)        :: useangles
         character(fnlen)        :: imagefile
         character(fnlen)        :: masterfile
+        character(fnlen)        :: Kosselmasterfile
         character(fnlen)        :: datafile
 end type BSENameListType
 
@@ -1349,6 +1350,7 @@ module NameListTypedefs
         character(1)            :: maskpattern
         character(3)            :: scalingmode
         character(3)            :: Notify
+        character(3)            :: similaritymetric
         !character(3)            :: eulerconvention
         !character(3)            :: outputformat
         character(1)            :: keeptmpfile