Primer commit con todos los laboratorios

1. First Lab/canny.py

@@ -0,0 +1,120 @@
+import numpy as np
+import scipy.ndimage as ndi
+from scipy.ndimage import (gaussian_filter, convolve,
+                           generate_binary_structure, binary_erosion, label)
+
+
+def smooth_with_function_and_mask(image, function, mask):
+
+    not_mask = np.logical_not(mask)
+    bleed_over = function(mask.astype(float))
+    masked_image = np.zeros(image.shape, image.dtype)
+    masked_image[mask] = image[mask]
+    smoothed_image = function(masked_image)
+    output_image = smoothed_image / (bleed_over + np.finfo(float).eps)
+    return output_image
+
+
+def canny(image, sigma, low_threshold, high_threshold, mask=None):
+
+    if mask is None:
+        mask = np.ones(image.shape, dtype=bool)
+    fsmooth = lambda x: gaussian_filter(x, sigma, mode='constant')
+    smoothed = smooth_with_function_and_mask(image, fsmooth, mask)
+    jsobel = ndi.sobel(smoothed, axis=1)
+    isobel = ndi.sobel(smoothed, axis=0)
+    abs_isobel = np.abs(isobel)
+    abs_jsobel = np.abs(jsobel)
+    magnitude = np.hypot(isobel, jsobel)
+
+    s = generate_binary_structure(2, 2)
+    eroded_mask = binary_erosion(mask, s, border_value=0)
+    eroded_mask = eroded_mask & (magnitude > 0)
+
+    local_maxima = np.zeros(image.shape,bool)
+    #----- 0 to 45 degrees ------
+    pts_plus = (isobel >= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel)
+    pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel)
+    pts = pts_plus | pts_minus
+    pts = eroded_mask & pts
+    # Get the magnitudes shifted left to make a matrix of the points to the
+    # right of pts. Similarly, shift left and down to get the points to the
+    # top right of pts.
+    c1 = magnitude[1:, :][pts[:-1, :]]
+    c2 = magnitude[1:, 1:][pts[:-1, :-1]]
+    m = magnitude[pts]
+    w = abs_jsobel[pts] / abs_isobel[pts]
+    c_plus = c2 * w + c1 * (1 - w) <= m
+    c1 = magnitude[:-1, :][pts[1:, :]]
+    c2 = magnitude[:-1, :-1][pts[1:, 1:]]
+    c_minus = c2 * w + c1 * (1 - w) <= m
+    local_maxima[pts] = c_plus & c_minus
+    #----- 45 to 90 degrees ------
+    # Mix diagonal and vertical
+    #
+    pts_plus = (isobel >= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel)
+    pts_minus = (isobel <= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel)
+    pts = pts_plus | pts_minus
+    pts = eroded_mask & pts
+    c1 = magnitude[:, 1:][pts[:, :-1]]
+    c2 = magnitude[1:, 1:][pts[:-1, :-1]]
+    m = magnitude[pts]
+    w = abs_isobel[pts] / abs_jsobel[pts]
+    c_plus = c2 * w + c1 * (1 - w) <= m
+    c1 = magnitude[:, :-1][pts[:, 1:]]
+    c2 = magnitude[:-1, :-1][pts[1:, 1:]]
+    c_minus = c2 * w + c1 * (1 - w) <= m
+    local_maxima[pts] = c_plus & c_minus
+    #----- 90 to 135 degrees ------
+    # Mix anti-diagonal and vertical
+    #
+    pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel <= abs_jsobel)
+    pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel <= abs_jsobel)
+    pts = pts_plus | pts_minus
+    pts = eroded_mask & pts
+    c1a = magnitude[:, 1:][pts[:, :-1]]
+    c2a = magnitude[:-1, 1:][pts[1:, :-1]]
+    m = magnitude[pts]
+    w = abs_isobel[pts] / abs_jsobel[pts]
+    c_plus = c2a * w + c1a * (1.0 - w) <= m
+    c1 = magnitude[:, :-1][pts[:, 1:]]
+    c2 = magnitude[1:, :-1][pts[:-1, 1:]]
+    c_minus = c2 * w + c1 * (1.0 - w) <= m
+    cc = np.logical_and(c_plus,c_minus)
+    local_maxima[pts] = c_plus & c_minus
+    #----- 135 to 180 degrees ------
+    # Mix anti-diagonal and anti-horizontal
+    #
+    pts_plus = (isobel <= 0) & (jsobel >= 0) & (abs_isobel >= abs_jsobel)
+    pts_minus = (isobel >= 0) & (jsobel <= 0) & (abs_isobel >= abs_jsobel)
+    pts = pts_plus | pts_minus
+    pts = eroded_mask & pts
+    c1 = magnitude[:-1, :][pts[1:, :]]
+    c2 = magnitude[:-1, 1:][pts[1:, :-1]]
+    m = magnitude[pts]
+    w = abs_jsobel[pts] / abs_isobel[pts]
+    c_plus = c2 * w + c1 * (1 - w) <= m
+    c1 = magnitude[1:, :][pts[:-1, :]]
+    c2 = magnitude[1:,:-1][pts[:-1,1:]]
+    c_minus = c2 * w + c1 * (1 - w) <= m
+    local_maxima[pts] = c_plus & c_minus
+    #
+    #---- Create two masks at the two thresholds.
+    #
+    high_mask = local_maxima & (magnitude >= high_threshold)
+    low_mask = local_maxima & (magnitude >= low_threshold)
+    #
+    # Segment the low-mask, then only keep low-segments that have
+    # some high_mask component in them
+    #
+    labels,count = label(low_mask, np.ndarray((3, 3),bool))
+    if count == 0:
+        return low_mask
+
+    sums = (np.array(ndi.sum(high_mask,labels,
+                             np.arange(count,dtype=np.int32) + 1),
+                     copy=False, ndmin=1))
+    good_label = np.zeros((count + 1,),bool)
+    good_label[1:] = sums > 0
+    output_mask = good_label[labels]
+    return output_mask

1. First Lab/canny_skimage.py

@@ -0,0 +1,89 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import ndimage as ndi
+from skimage.util import random_noise
+from skimage import feature
+from skimage import filters
+
+
+# Generate noisy image of a square
+image = np.zeros((128, 128), dtype=float)
+image[:, :-64] = 1
+
+#image1 = ndi.rotate(image, 15, mode='constant')
+#image_noise = ndi.gaussian_filter(image, 4)
+
+lvar = filters.gaussian(image, sigma=4) + 1e-10
+image_random = random_noise(np.zeros((128, 128), dtype=float), mode='localvar', local_vars=lvar*0.5)
+image_random1 = random_noise(image, mode='localvar', local_vars=lvar*0.5)
+image_canny= image_random + image
+
+image_canny= image_random + image
+#image_random = random_noise(image, mode='localvar', mean=0.1)
+
+# Compute the Canny filter for two values of sigma
+edges1 = feature.canny(image_canny)
+edges2 = feature.canny(image, sigma=3)
+
+def fom(img, img_gold_std, alpha = DEFAULT_ALPHA):
+
+
+    # To avoid oversmoothing, we apply canny edge detection with very low
+    # standard deviation of the Gaussian kernel (sigma = 0.1).
+    edges_img = canny(img, 0.1, 20, 50)
+    edges_gold = canny(img_gold_std, 0.1, 20, 50)
+
+    # Compute the distance transform for the gold standard image.
+    dist = distance_transform_edt(np.invert(edges_gold))
+
+    fom = 1.0 / np.maximum(
+        np.count_nonzero(edges_img),
+        np.count_nonzero(edges_gold))
+
+    N, M = img.shape
+
+    for i in range(0, N):
+        for j in range(0, M):
+            if edges_img[i, j]:
+                fom += 1.0 / ( 1.0 + dist[i, j] * dist[i, j] * alpha)
+
+    fom /= np.maximum(
+        np.count_nonzero(edges_img),
+        np.count_nonzero(edges_gold))    
+
+    return fom
+
+
+
+
+
+
+
+
+# display results
+fig, ax = plt.subplots(nrows=1, ncols=5, figsize=(8, 3))
+
+ax[0].imshow(image, cmap='gray')
+ax[0].set_title('image', fontsize=20)
+
+ax[1].imshow(image_random, cmap='gray')
+ax[1].set_title('ruido', fontsize=20)
+
+ax[1].imshow(image_canny, cmap='gray')
+ax[1].set_title('suma de ruido', fontsize=20)
+
+ax[2].imshow(edges1, cmap='gray')
+ax[2].set_title('canny', fontsize=20)
+
+#ax[3].imshow(edges1, cmap='gray')
+#ax[3].set_title('Canny filter', fontsize=20)
+
+#ax[4].imshow(image_, cmap='gray')
+#ax[4].set_title('noise', fontsize=20)
+
+
+for a in ax:
+    a.axis('off')
+
+fig.tight_layout()
+plt.show()

1. First Lab/kuan.py

@@ -0,0 +1,74 @@
+
+import numpy as np
+import scipy.stats as stats
+
+DEFAULT_WINDOW_SIZE = 7
+DEFAULT_CU = 0.25
+
+def variation(window):
+    """
+    Calculates the coefficient of variation, that is, the ratio of the standard
+    deviation to the mean.
+    """
+    return stats.variation(window, None)
+
+def weight(window, cu):
+    """
+    Calculates the weight parameter of the Kuan filter. The weight defines
+    weather the Kuan filter will act as an identity filter or as a mean filter.
+    """
+    ci = variation(window)
+    ci_2 = ci * ci
+    cu_2 = cu * cu
+
+    if cu_2 > ci_2:
+        w = 0.0 # do not allow negative values
+    else:
+        w = (1.0 - (cu_2 / ci_2)) / (1.0 + cu_2)
+
+    return w
+
+def filter(img, window_size = DEFAULT_WINDOW_SIZE, cu = DEFAULT_CU):
+    """
+    Filter the given image, using the passed in window size and coefficient of
+    variation of a smooth area.
+    Input arguments:
+        img: a numpy array representing the image to be filtered.
+        window_size: the size.
+        cu: the coefficient of variation of a smooth region of img.
+    Output:
+        the filtered image.
+    """
+
+    img = np.float64(img)
+    img_filtered = np.zeros_like(img)
+
+    N, M = img.shape
+    win_offset = window_size / 2
+
+    for i in range(0, N):
+        xleft = i - win_offset
+        xright = i + win_offset
+
+        # border conditions in x axis
+        if xleft < 0:
+            xleft = 0
+        if xright >= N:
+            xright = N
+
+        for j in range(0, M):
+            yup = j - win_offset
+            ydown = j + win_offset
+
+            # border conditions in y axis
+            if yup < 0:
+                yup = 0
+            if ydown >= M:
+                ydown = M
+
+            # Calculate new value for pixel [i,j]
+            window = img[xleft:xright, yup:ydown]
+            w = weight(window, cu)
+            img_filtered[i, j] = round((img[i, j] * w) + (window.mean() * (1.0 - w)))
+
+    return img_filtered