Initial attempt at including diffusion

Signed-off by: David Rowenhorst <[email protected]>

pyebsdindex/nlpar.py

@@ -60,4 +60,8 @@
 
 __all__ = [
     "NLPAR",
-]
+]
+
+class DIFF_NLPAR(NLPAR):
+  def __init__(self, **kwargs):
+    pass

pyebsdindex/nlpar_cpu.py

@@ -44,10 +44,12 @@
 
 
 class NLPAR:
-  def __init__(self, filename=None,  lam=0.7, searchradius=3,dthresh=0.0, nrows = None, ncols = None, **kwargs):
+  def __init__(self, filename=None,  lam=0.7, searchradius=3,dthresh=0.0, diff_offset=1.0,
+               nrows = None, ncols = None, **kwargs):
     self.lam = lam
     self.searchradius = searchradius
     self.dthresh = dthresh
+    self.diff_offset = diff_offset,
     self.filepath = None
     self.hdfdatapath = None
     self.filepathout = None
@@ -304,7 +306,7 @@ def d2norm(d2, n2, dij, sigma):
     return np.mean(lamopt_values, axis = 0).flatten()
 
   def calcnlpar(self, chunksize=0, searchradius=None, lam = None, dthresh = None, saturation_protect=True, automask=True,
-               filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,verbose=2,
+               filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,verbose=2, diff_offset=None,
                 **kwargs):
 
     if lam is not None:
@@ -313,12 +315,16 @@ def calcnlpar(self, chunksize=0, searchradius=None, lam = None, dthresh = None,
     if dthresh is not None:
       self.dthresh = dthresh
 
+    if diff_offset is not None:
+      self.diff_offset = diff_offset
+
     if searchradius is not None:
       self.searchradius = searchradius
 
     lam = np.float32(self.lam)
     dthresh = np.float32(self.dthresh)
     sr = np.int64(self.searchradius)
+    diff_offset = np.float32(self.diff_offset)
 
     if filename is not None:
       self.setfile(filepath=filename)
@@ -429,7 +435,7 @@ def calcnlpar(self, chunksize=0, searchradius=None, lam = None, dthresh = None,
       #dataout = data
 
       dataout = self.nlpar_nb(data,lam, sr, dthresh, sigchunk,
-                              rowcountread,ncols,indices,saturation_protect)
+                              rowcountread,ncols,indices,saturation_protect, diff_offset=diff_offset)
 
       dataout = dataout.reshape(rowcountread, ncols, phw)
       dataout = dataout[j-rowstartread:, :, : ]
@@ -646,7 +652,7 @@ def sigma_numba(data, nn, nrows, ncols, rowstartcount, colstartcount, indices, s
 
   @staticmethod
   @numba.jit(nopython=True,cache=True,fastmath=False,parallel=True)
-  def nlpar_nb(data,lam, sr, dthresh, sigma, nrows,ncols,indices,saturation_protect=True):
+  def nlpar_nb(data,lam, sr, dthresh, sigma, nrows,ncols,indices,saturation_protect=True, diff_offset = np.float32(1.0)):
     def getpairid(idx0, idx1):
       idx0_t = int(idx0)
       idx1_t = int(idx1)
@@ -661,7 +667,7 @@ def getpairid(idx0, idx1):
     shpdata = data.shape
     shpind = indices.shape
     winsz = np.int32((2*sr+1)**2)
-
+    diff_step =  np.ones((winsz), dtype=np.float32)
 
     mxval = np.max(data)
     if saturation_protect == False:
@@ -690,7 +696,9 @@ def getpairid(idx0, idx1):
 
             if indx_nn == indx_0:
               weights[counter] = np.float32(-1.0e6)
+              diff_step[counter] = 1.0 / diff_offset
             else:
+              diff_step[counter] =  diff_offset
               pairid = getpairid(indx_0, indx_nn)
               if pairid in pairdict:
                 weights[counter] = pairdict[pairid]
@@ -720,6 +728,7 @@ def getpairid(idx0, idx1):
 
           weights[i_nn] = np.maximum(weights[i_nn]-dthresh, numba.float32(0.0))
           weights[i_nn] = np.exp(-1.0 * weights[i_nn] * lam2)
+          weights[i_nn] *= diff_step[i_nn]
           sum += weights[i_nn]
 
         for i_nn in range(winsz):

pyebsdindex/opencl/clnlpar.cl

@@ -268,7 +268,8 @@ __kernel void calcnlpar(
       const long npatpoint, 
       const float maxlim, 
       const float lam2, 
-      const float dthresh){
+      const float dthresh,
+      const float diff_offset){
   //IDs of work-item represent x and y coordinates in image
   //const long4 calclim =  crlimits[0];
   const long x = get_global_id(0)+crlimits[0];
@@ -293,11 +294,13 @@ __kernel void calcnlpar(
 
   float d[512]; // taking a risk here that noone will want a SR > 10
   float n[512]; 
+  float diff_step[512];
 
 
   for(j=0; j < nnn; ++j){
       d[j] = 0.0;
       n[j] = 1.0e-6; 
+      diff_step[j] = 1.0 ;
   }
 
 
@@ -340,7 +343,14 @@ __kernel void calcnlpar(
               d1 *= mask1; 
 
               dd = sum16(&d1);
-              dd = (indx_ij == indx0) ? -1.0 : dd; // mark the center point 
+              if (indx_ij == indx0) {
+                dd = -1.0;
+                diff_step[count] = 1.0 / diff_offset;
+              } else{
+                dd = dd;
+                diff_step[count] = diff_offset;
+              }
+              //dd = (indx_ij == indx0) ? -1.0 : dd; // mark the center point 
 
               d[count] += dd;
               n[count] += sum16(&mask1);
@@ -388,10 +398,12 @@ __kernel void calcnlpar(
                 nn = 1.0;
               }
 
-              dd -= dthresh; 
+              dd -= dthresh;
+
               dd = dd >= 0.0 ? dd : 0.0; 
 
               dd = exp(-1.0*dd*lam2); 
+              dd *= diff_step[count]; 
               sum += dd; 
               d[count] = dd;
 

pyebsdindex/opencl/nlpar_cl.py

@@ -304,8 +304,8 @@ def calcsigma_cl(self,nn=1,saturation_protect=True,automask=True, normalize_d=Fa
 
 
   def calcnlpar_cl(self, searchradius=None, lam = None, dthresh = None, saturation_protect=True, automask=True,
-                   filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,
-                   gpu_id = None, verbose=2, **kwargs):
+                   filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,diff_offset= None,
+                   gpu_id = None, verbose=2,   **kwargs):
 
     class OpenCLClalcError(Exception):
       pass
@@ -315,9 +315,13 @@ class OpenCLClalcError(Exception):
 
     if dthresh is not None:
       self.dthresh = dthresh
+
     if self.dthresh is None:
       self.dthresh = 0.0
 
+    if diff_offset is not None:
+      self.diff_offset = diff_offset
+
     if searchradius is not None:
       self.searchradius = searchradius
 
@@ -330,6 +334,7 @@ class OpenCLClalcError(Exception):
     lam = np.float32(self.lam)
     dthresh = np.float32(self.dthresh)
     sr = np.int64(self.searchradius)
+    diff_offset = np.float32(self.diff_offset)
 
     if filename is not None:
       self.setfile(filepath=filename)
@@ -561,7 +566,8 @@ class OpenCLClalcError(Exception):
                                  np.int64(npat_point),
                                  np.float32(mxval),
                                  np.float32(1.0/lam**2),
-                                 np.float32(dthresh) )
+                                 np.float32(dthresh),
+                                 np.float32(diff_offset))
 
           data = data.astype(np.float32) # prepare to receive data back from GPU
           data.reshape(-1)[:] = 0.0

pyebsdindex/opencl/nlpar_clray.py

@@ -370,8 +370,8 @@ def _sigmachunkcalc_cl(self, data, calclim, clparams=None, saturation_protect=Tr
 
 
 
-  def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, saturation_protect=True, automask=True,
-                filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,
+  def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, saturation_protect=True, automask = True,
+                filename=None, fileout=None, reset_sigma=False, backsub = False, rescale = False,diff_offset = None,
                 verbose = 2, gpu_id = None, **kwargs):
 
     if lam is not None:
@@ -384,6 +384,9 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, saturat
     if self.dthresh is None:
       self.dthresh = 0.0
 
+    if diff_offset is not None:
+      self.diff_offset = diff_offset
+
     if searchradius is not None:
       self.searchradius = searchradius
 
@@ -396,6 +399,7 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, saturat
     lam = np.float32(self.lam)
     dthresh = np.float32(self.dthresh)
     sr = np.int64(self.searchradius)
+    diff_offset = np.float32(self.diff_offset)
 
     if filename is not None:
       self.setfile(filepath=filename)
@@ -477,7 +481,8 @@ def calcnlpar_clray(self, searchradius=None, lam = None, dthresh = None, saturat
                                          reset_sigma=reset_sigma,
                                          backsub = backsub,
                                          rescale = rescale,
-                                         gpu_id= gpu_id)
+                                         gpu_id= gpu_id,
+                                         diff_offset=diff_offset)
 
     target_mem = clparams.gpu[gpu_id].max_mem_alloc_size//6
     max_mem = clparams.gpu[gpu_id].global_mem_size*0.4
@@ -601,6 +606,7 @@ def _nlparchunkcalc_cl(self, data, calclim, clparams=None):
     sr = np.int64(self.searchradius)
     nnn = int((2 * sr + 1) ** 2)
     dthresh = np.float32(self.dthresh)
+    diff_offset = np.float32(self.diff_offset)
     #print(chunks[2], chunks[3])
     #print(lam, sr, dthresh)
 
@@ -662,7 +668,8 @@ def _nlparchunkcalc_cl(self, data, calclim, clparams=None):
                   np.int64(npat_point),
                   np.float32(mxval),
                   np.float32(1.0 / lam ** 2),
-                  np.float32(dthresh))
+                  np.float32(dthresh), 
+                  np.float32(diff_offset))
 
     data = data.astype(np.float32)  # prepare to receive data back from GPU
     data.reshape(-1)[:] = 0.0