title: Using inpaint function in OpenCV via Python to interpolate broken river-data in a watershed
body: <h2>Background</h2>
  
  <p>A raster file collected via LIDAR records the topography of a watershed. To properly model the watershed, the river must appear continuous without any breaks or interruptions. The roads in the raster file appear like dams that interrupt the river as seen in the picture below</p>
  
  <p><a href="http://i.stack.imgur.com/6Whor.png" rel="nofollow">Specific Area Under Consideration in the Watershed</a></p>
  
  <h2>Objective</h2>
  
  <p>These river breaks are the main problem and I am trying but failing to remove them.</p>
  
  <h2>Approach</h2>
  
  <p>Via Python, I used the various tools and prebuilt functions in the OpenCV library. The primary function I used in this approach is the cv2.inpaint function. This function takes in an <strong>image file</strong> and a <strong>mask file</strong> and interpolates the original image wherever the mask file pixels are nonzero.</p>
  
  <p>The main step here is determining the mask file which I did by detecting the corners at the break in the river. The mask file will guide the inpaint function to fill in the pixels according to the patterns in the surrounding pixels.</p>
  
  <h2>Problem</h2>
  
  <p>My issue is that this happens from all directions whereas I require it to only extrapolate pixel data from the river itself. The image below shows the flawed result: inpaint works but it considers data from outside the river too.</p>
  
  <p><a href="http://i.stack.imgur.com/itTIX.png" rel="nofollow">Inpainted Result</a></p>
  
  <p>Here is my code if you are so kind as to help:</p>
  
  <pre><code>import scipy.io as sio
  import numpy as np
  from matplotlib import pyplot as plt
  import cv2
  
  matfile = sio.loadmat('box.mat') ## box.mat original raster file linked below
  ztopo = matfile['box']
  
  #Crop smaller section for analysis
  ztopo2 = ztopo[200:500, 1400:1700]
  
  
  ## Step 1) Isolate river
  river = ztopo2.copy()
  river[ztopo2&lt;217.5] = 0
  #This will become problematic for entire watershed w/o proper slicing
  
  
  ## Step 2) Detect Corners
  dst = cv2.cornerHarris(river,3,7,0.04)
  # cornerHarris arguments adjust qualities of corner markers
  
  # Dilate Image (unnecessary)
  #dst = cv2.dilate(dst,None)  
  
  # Threshold for an optimal value, it may vary depending on the image.
  # This adjusts what defines a corner
  river2 = river.copy()
  river2[dst&gt;0.01*dst.max()]=[150]
  
  
  ## Step 3) Remove river and keep corners
  
  #Initiate loop to isolate detected corners
  i=0
  j=0
  rows,columns = river2.shape
  for i in np.arange(rows):
      for j in np.arange(columns):
          if river2[i,j] != 150:
              river2[i,j] = 0
          j = j + 1
      i = i + 1
  
  # Save corners as new image for import during next step.
  # Import must be via cv2 as thresholding and contour detection can only work on BGR files. Sio import in line 6 (matfile = sio.loadmat('box.mat')) imports 1 dimensional image rather than BGR.
  cv2.imwrite('corners.png', river2)
  
  
  ## Step 4) Create mask image by defining and filling a contour around the previously detected corners
  
  #Step 4 code retrieved from http://dsp.stackexchange.com/questions/2564/opencv-c-connect-nearby-contours-based-on-distance-between-them
  #Article: OpenCV/C++ connect nearby contours based on distance between them
  
  #Initiate function to specify features of contour connections
  def find_if_close(cnt1,cnt2):
      row1,row2 = cnt1.shape[0],cnt2.shape[0]
      for i in xrange(row1):
          for j in xrange(row2):
              dist = np.linalg.norm(cnt1[i]-cnt2[j])
              if abs(dist) &lt; 50 :
                  return True
              elif i==row1-1 and j==row2-1:
                  return False
  
  #import image of corners created in step 3 so thresholding can function properly
  img = cv2.imread('corners.png')
  
  #Thresholding and Finding contours only works on grayscale image
  gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
  ret,thresh = cv2.threshold(gray,127,255,cv2.THRESH_BINARY)
  contours,hier = cv2.findContours(thresh,cv2.RETR_EXTERNAL,2)
  
  LENGTH = len(contours)
  status = np.zeros((LENGTH,1))
  
  for i,cnt1 in enumerate(contours):
      x = i    
      if i != LENGTH-1:
          for j,cnt2 in enumerate(contours[i+1:]):
              x = x+1
              dist = find_if_close(cnt1,cnt2)
              if dist == True:
                  val = min(status[i],status[x])
                  status[x] = status[i] = val
              else:
                  if status[x]==status[i]:
                      status[x] = i+1
  
  unified = []
  maximum = int(status.max())+1
  for i in xrange(maximum):
      pos = np.where(status==i)[0]
      if pos.size != 0:
          cont = np.vstack(contours[i] for i in pos)
          hull = cv2.convexHull(cont) # I don't know why/how this is used
          unified.append(hull)
  
  
  cv2.drawContours(img,unified,-1,(0,255,0),1) #Last argument specifies contour width
  cv2.drawContours(thresh,unified,-1,255,-1)
  
  # Thresh is the filled contour while img is the contour itself
  # The contour surrounds the corners
  
  #cv2.imshow('pic', thresh) #Produces black and white image
  
  
  ## Step 5) Merge via inpaint
  river = np.uint8(river)
  ztopo2 = np.uint8(ztopo2)
  
  thresh[thresh&gt;0] = 1
  #new = river.copy()
  merged = cv2.inpaint(river,thresh,12,cv2.INPAINT_TELEA)
  
  plt.imshow(merged)
  plt.colorbar()
  plt.show()
  </code></pre>