Learning Convolutions by Inventing Computer Vision

The color image shows three channels (RGB) producing full color, while the grayscale image has only one channel representing brightness. They show the same scene, but the grayscale version loses all color information.

# Display the color image
show_images([color_img], titles=['Color Image'])

# Display the grayscale image
show_images([gray_img], titles=['Grayscale Image'], cmap_list=['gray'])

# Display both side by side
show_images([color_img, gray_img],
            titles=['Color Image', 'Grayscale Image'],
            cmap_list=[None, 'gray'])

A color image has shape (H, W, 3). The third dimension is 3 because each pixel stores three values: Red, Green, and Blue. A grayscale image has shape (H, W) with no third dimension because each pixel is a single brightness value. Values range from 0 to 255 (uint8). color_arr[0, 0] gives the RGB triplet of the top-left pixel; color_arr[0, 0, 0] gives just its red intensity.

color_arr = np.array(color_img)
gray_arr  = np.array(gray_img)

print("--- Color Array ---")
print("Shape:", color_arr.shape)   # (H, W, 3)
print("Dtype:", color_arr.dtype)   # uint8
print("Min:",   color_arr.min(), "  Max:", color_arr.max())
print("Pixel at [0,0]:", color_arr[0, 0])         # [R, G, B] of top-left pixel
print("Red value at [0,0]:", color_arr[0, 0, 0])  # Red channel only

print("\n--- Grayscale Array ---")
print("Shape:", gray_arr.shape)    # (H, W)
print("Dtype:", gray_arr.dtype)    # uint8

print("\nFirst 5x5 block of grayscale:")
print(gray_arr[:5, :5])

height = color_arr.shape[0]
width  = color_arr.shape[1]
print(f"Image size: {height} rows x {width} columns")

# Center pixel
center_row = height // 2
center_col = width // 2
print("Center pixel:", color_arr[center_row, center_col])

# Bottom-right corner pixel
print("Bottom-right pixel:", color_arr[height - 1, width - 1])

# First row of grayscale
print("First row of grayscale:", gray_arr[0, :])
print("Number of values:", len(gray_arr[0, :]))

# Average brightness using ONLY Python loops
total = 0
count = 0
h, w = gray_arr.shape
for r in range(h):
    for c in range(w):
        total += int(gray_arr[r, c])
        count += 1

avg_brightness = total / count
print(f"Average brightness (manual loop): {avg_brightness:.2f}")
print(f"Verify with numpy:                {gray_arr.mean():.2f}")

# Top-left 20x20 patch of grayscale
gray_patch = gray_arr[:20, :20]
show_images([gray_patch], titles=['Grayscale top-left 20x20'], cmap_list=['gray'])

# Top-left 20x20 patch of color
color_patch = color_arr[:20, :20]
show_images([color_patch], titles=['Color top-left 20x20'])

# Middle 20x20 patch
mid_r = gray_arr.shape[0] // 2
mid_c = gray_arr.shape[1] // 2
gray_mid  = gray_arr[mid_r-10:mid_r+10, mid_c-10:mid_c+10]
color_mid = color_arr[mid_r-10:mid_r+10, mid_c-10:mid_c+10]
show_images([gray_mid, color_mid],
            titles=['Gray middle 20x20', 'Color middle 20x20'],
            cmap_list=['gray', None])

test[1, 0, 1] is the Green channel (index 1) of the pixel at row 1, column 0. Since row 1 is all green [0, 255, 0], the value is 255.

# 1. Top-left pixel
print("Top-left pixel R,G,B:", color_arr[0, 0])

# 2. Build the 4x4 test array
test = np.zeros((4, 4, 3), dtype=np.uint8)
test[0, :] = [255, 0, 0]      # Row 0: red
test[1, :] = [0, 255, 0]      # Row 1: green
test[2, :] = [0, 0, 255]      # Row 2: blue
test[3, :] = [255, 255, 255]  # Row 3: white
show_images([test], titles=['Test 4x4 color array'])

# 3. Channel access
print("test[0,0,0] (Red of first pixel):   ", test[0, 0, 0])   # 255
print("test[1,0,1] (Green of second row):  ", test[1, 0, 1])   # 255

Bright areas in the red channel mean those pixels have a high red component. Dark areas have little red. Each channel is independent of the others.

height = color_arr.shape[0]
width  = color_arr.shape[1]

# Extract Red channel
red_channel = np.zeros((height, width), dtype=np.uint8)
for r in range(height):
    for c in range(width):
        red_channel[r, c] = color_arr[r, c, 0]

show_images([red_channel], titles=['Red Channel (as grayscale)'], cmap_list=['gray'])

# Extract Green and Blue channels the same way
green_channel = np.zeros((height, width), dtype=np.uint8)
blue_channel  = np.zeros((height, width), dtype=np.uint8)
for r in range(height):
    for c in range(width):
        green_channel[r, c] = color_arr[r, c, 1]
        blue_channel[r, c]  = color_arr[r, c, 2]

# Display all three channels side by side
show_images(
    [red_channel, green_channel, blue_channel],
    titles=['Red', 'Green', 'Blue'],
    cmap_list=['Reds', 'Greens', 'Blues']
)

The recombined image looks identical to the original because each pixel's RGB values are reconstructed exactly (R + 0 + 0 = R, etc.). Using int() and min(255, ...) prevents uint8 overflow.

# Build red_image: (H, W, 3) with only red channel filled
red_image = np.zeros((height, width, 3), dtype=np.uint8)
for r in range(height):
    for c in range(width):
        red_image[r, c, 0] = red_channel[r, c]

show_images([red_image], titles=['Red channel as color'])

# Build green_image and blue_image similarly
green_image = np.zeros((height, width, 3), dtype=np.uint8)
blue_image  = np.zeros((height, width, 3), dtype=np.uint8)
for r in range(height):
    for c in range(width):
        green_image[r, c, 1] = green_channel[r, c]
        blue_image[r, c, 2]  = blue_channel[r, c]

show_images([red_image, green_image, blue_image],
            titles=['Red', 'Green', 'Blue'])

# Add them together using a loop (handle overflow)
combined = np.zeros((height, width, 3), dtype=np.uint8)
for r in range(height):
    for c in range(width):
        for ch in range(3):
            val = int(red_image[r, c, ch]) + int(green_image[r, c, ch]) + int(blue_image[r, c, ch])
            combined[r, c, ch] = min(255, val)

show_images([color_arr, combined],
            titles=['Original', 'Recombined R+G+B'])

Three common ideas for combining R, G, B into a single brightness value:

1. Average: gray = (R + G + B) / 3. Treats all three channels equally.
2. Lightness: gray = (max(R, G, B) + min(R, G, B)) / 2. Uses only the extremes.
3. Luminosity (weighted): gray = 0.299R + 0.587G + 0.114B. Matches human perception by weighting green most heavily and blue least.

The simple average assumes R, G, and B contribute equally to perceived brightness, which is not true for human vision. Our eyes are most sensitive to green and least to blue.

def rgb_to_gray_average(arr):
    h, w = arr.shape[0], arr.shape[1]
    result = np.zeros((h, w), dtype=np.uint8)
    for r in range(h):
        for c in range(w):
            R, G, B = int(arr[r, c, 0]), int(arr[r, c, 1]), int(arr[r, c, 2])
            gray = int((R + G + B) / 3)
            result[r, c] = max(0, min(255, gray))
    return result

def rgb_to_gray_lightness(arr):
    h, w = arr.shape[0], arr.shape[1]
    result = np.zeros((h, w), dtype=np.uint8)
    for r in range(h):
        for c in range(w):
            R, G, B = int(arr[r, c, 0]), int(arr[r, c, 1]), int(arr[r, c, 2])
            gray = int((max(R, G, B) + min(R, G, B)) / 2)
            result[r, c] = max(0, min(255, gray))
    return result

def rgb_to_gray_luminosity(arr):
    h, w = arr.shape[0], arr.shape[1]
    result = np.zeros((h, w), dtype=np.uint8)
    for r in range(h):
        for c in range(w):
            R, G, B = int(arr[r, c, 0]), int(arr[r, c, 1]), int(arr[r, c, 2])
            gray = int(0.299 * R + 0.587 * G + 0.114 * B)
            result[r, c] = max(0, min(255, gray))
    return result

def rgb_to_gray_mymethod(arr):
    """Desaturation: use the median of R, G, B."""
    h, w = arr.shape[0], arr.shape[1]
    result = np.zeros((h, w), dtype=np.uint8)
    for r in range(h):
        for c in range(w):
            vals = [int(arr[r, c, 0]), int(arr[r, c, 1]), int(arr[r, c, 2])]
            vals.sort()
            result[r, c] = vals[1]  # median
    return result

The luminosity method produces the most perceptually accurate grayscale because it mirrors how human eyes perceive brightness. Green gets the most weight (0.587) because our eyes have more green-sensitive cone cells. Blue gets the least (0.114) because blue light appears dimmer to us. The biggest differences appear in areas with strong color saturation, especially reds and blues.

# Apply all methods
gray_average   = rgb_to_gray_average(color_arr)
gray_lightness = rgb_to_gray_lightness(color_arr)
gray_luminosity = rgb_to_gray_luminosity(color_arr)
gray_mymethod  = rgb_to_gray_mymethod(color_arr)

# Display all side by side
show_images(
    [color_arr, gray_average, gray_lightness, gray_luminosity, gray_mymethod],
    titles=['Original', 'Average', 'Lightness', 'Luminosity', 'Median'],
    cmap_list=[None, 'gray', 'gray', 'gray', 'gray']
)

# Compute pixel-wise difference between average and luminosity
h, w = gray_average.shape
diff = np.zeros((h, w), dtype=np.uint8)
for r in range(h):
    for c in range(w):
        diff[r, c] = abs(int(gray_average[r, c]) - int(gray_luminosity[r, c]))

show_images([diff], titles=['Difference: Average vs Luminosity'], cmap_list=['hot'])

1. Output pixel (r, c) corresponds to the 2x2 block of input pixels starting at (2*r, 2*c): namely (2r, 2c), (2r, 2c+1), (2r+1, 2c), (2r+1, 2c+1).

2. Three possible values to assign: (a) just take one of them, e.g. the top-left value (nearest neighbor / subsampling); (b) average all four values (average pooling); (c) take the maximum of the four (max pooling).

3. Yes, there is information loss. Four pixel values are reduced to one, so you cannot perfectly reconstruct the original. Different 2x2 blocks can produce the same output value.

def downsample_nearest(img):
    """Downsamples (H, W) to (H//2, W//2) by taking every other pixel."""
    h, w = img.shape
    out_h, out_w = h // 2, w // 2
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            result[r, c] = img[2 * r, 2 * c]
    return result

def downsample_average(img):
    """Downsamples by averaging each 2x2 block."""
    h, w = img.shape
    out_h, out_w = h // 2, w // 2
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            val = (int(img[2*r, 2*c]) + int(img[2*r, 2*c+1])
                 + int(img[2*r+1, 2*c]) + int(img[2*r+1, 2*c+1])) / 4
            result[r, c] = int(val)
    return result

def downsample_max(img):
    """Downsamples by taking the maximum of each 2x2 block."""
    h, w = img.shape
    out_h, out_w = h // 2, w // 2
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            result[r, c] = max(int(img[2*r, 2*c]),   int(img[2*r, 2*c+1]),
                               int(img[2*r+1, 2*c]), int(img[2*r+1, 2*c+1]))
    return result

def downsample_min(img):
    """Downsamples by taking the minimum of each 2x2 block."""
    h, w = img.shape
    out_h, out_w = h // 2, w // 2
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            result[r, c] = min(int(img[2*r, 2*c]),   int(img[2*r, 2*c+1]),
                               int(img[2*r+1, 2*c]), int(img[2*r+1, 2*c+1]))
    return result

# Test
small = downsample_nearest(img_100)
print("Output shape:", small.shape)  # (50, 50)

Nearest neighbor is the sharpest but can miss detail from skipped pixels. Average pooling is smoothest because it blends all four values. Max pooling preserves the brightest activations in each block, making it preferred in neural networks because strong feature responses (e.g., an edge) are kept even if neighboring pixels are weaker.

small_nearest = downsample_nearest(img_100)
small_average = downsample_average(img_100)
small_max     = downsample_max(img_100)
small_min     = downsample_min(img_100)

show_images(
    [img_100, small_nearest, small_average, small_max, small_min],
    titles=['Original 100x100', 'Nearest', 'Average', 'Max', 'Min'],
    cmap_list=['gray'] * 5
)

# Difference map: nearest vs average
h, w = small_nearest.shape
diff_down = np.zeros((h, w), dtype=np.uint8)
for r in range(h):
    for c in range(w):
        diff_down[r, c] = abs(int(small_nearest[r, c]) - int(small_average[r, c]))

show_images([diff_down], titles=['Difference: Nearest vs Average'], cmap_list=['hot'])

Strategy A — Nearest neighbor replication: Each output pixel (r, c) copies from the closest input pixel (r//2, c//2). Every input pixel is "stretched" into a 2x2 block. Simple but blocky.

Strategy B — Interpolation: Output pixel (r, c) maps to input position (r/2, c/2), which may be fractional. Blend between the surrounding input pixels using their distances as weights. Smoother but slightly blurry.

def upsample_nearest(img, scale=2):
    """Upsamples (H, W) to (H*scale, W*scale) by pixel replication."""
    h, w = img.shape
    out_h, out_w = h * scale, w * scale
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            result[r, c] = img[r // scale, c // scale]
    return result

def lerp(a, b, t):
    """Linear interpolation: t=0 gives a, t=1 gives b."""
    return a + t * (b - a)

def upsample_bilinear(img, scale=2):
    """Upsamples using bilinear interpolation."""
    h, w = img.shape
    out_h, out_w = h * scale, w * scale
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            in_r = r / scale
            in_c = c / scale
            r0 = int(in_r)
            c0 = int(in_c)
            r1 = min(r0 + 1, h - 1)
            c1 = min(c0 + 1, w - 1)
            dr = in_r - r0
            dc = in_c - c0
            top    = lerp(float(img[r0, c0]), float(img[r0, c1]), dc)
            bottom = lerp(float(img[r1, c0]), float(img[r1, c1]), dc)
            value  = lerp(top, bottom, dr)
            result[r, c] = max(0, min(255, int(value)))
    return result

def upsample_repeat_rows_cols(img, scale=2):
    """Upsamples by repeating each row and column."""
    h, w = img.shape
    out_h, out_w = h * scale, w * scale
    result = np.zeros((out_h, out_w), dtype=np.uint8)
    for r in range(out_h):
        for c in range(out_w):
            result[r, c] = img[r // scale, c // scale]
    return result

big_nearest = upsample_nearest(img_50)
print("Output shape:", big_nearest.shape)  # (100, 100)

Nearest neighbor is blocky because each input pixel occupies a 2x2 block with no transition. Bilinear is smoother because it blends between neighbors, but the trade-off is slight blurriness. Bilinear typically achieves a lower MAE because averaging reduces extreme errors.

big_nearest  = upsample_nearest(img_50)
big_bilinear = upsample_bilinear(img_50)

show_images(
    [img_100, big_nearest, big_bilinear],
    titles=['Original 100x100', 'Nearest (100x100)', 'Bilinear (100x100)'],
    cmap_list=['gray'] * 3
)

# Mean absolute error (using loops)
def mean_absolute_error(img_a, img_b):
    h, w = img_a.shape
    total = 0
    for r in range(h):
        for c in range(w):
            total += abs(int(img_a[r, c]) - int(img_b[r, c]))
    return total / (h * w)

mae_nearest  = mean_absolute_error(big_nearest, img_100)
mae_bilinear = mean_absolute_error(big_bilinear, img_100)
print(f"MAE (Nearest):  {mae_nearest:.4f}")
print(f"MAE (Bilinear): {mae_bilinear:.4f}")

def resize_nearest(img, target_h, target_w):
    """Resizes img to (target_h, target_w) using nearest neighbor."""
    h, w = img.shape
    result = np.zeros((target_h, target_w), dtype=np.uint8)
    for r in range(target_h):
        for c in range(target_w):
            in_r = int(r * h / target_h)
            in_c = int(c * w / target_w)
            in_r = min(in_r, h - 1)
            in_c = min(in_c, w - 1)
            result[r, c] = img[in_r, in_c]
    return result

def resize_bilinear(img, target_h, target_w):
    """Resizes img to (target_h, target_w) using bilinear interpolation."""
    h, w = img.shape
    result = np.zeros((target_h, target_w), dtype=np.uint8)
    for r in range(target_h):
        for c in range(target_w):
            in_r = r * (h - 1) / max(1, target_h - 1)
            in_c = c * (w - 1) / max(1, target_w - 1)
            r0 = int(in_r)
            c0 = int(in_c)
            r1 = min(r0 + 1, h - 1)
            c1 = min(c0 + 1, w - 1)
            dr = in_r - r0
            dc = in_c - c0
            top    = lerp(float(img[r0, c0]), float(img[r0, c1]), dc)
            bottom = lerp(float(img[r1, c0]), float(img[r1, c1]), dc)
            value  = lerp(top, bottom, dr)
            result[r, c] = max(0, min(255, int(value)))
    return result

# Resize to 150x150
img_150_nearest  = resize_nearest(img_50, 150, 150)
img_150_bilinear = resize_bilinear(img_50, 150, 150)
show_images([img_50, img_150_nearest, img_150_bilinear],
            titles=['Original 50x50', 'Nearest 150x150', 'Bilinear 150x150'],
            cmap_list=['gray'] * 3)

# Resize to 75x75
img_75_nearest  = resize_nearest(img_50, 75, 75)
img_75_bilinear = resize_bilinear(img_50, 75, 75)
show_images([img_50, img_75_nearest, img_75_bilinear],
            titles=['Original 50x50', 'Nearest 75x75', 'Bilinear 75x75'],
            cmap_list=['gray'] * 3)

Nearest neighbor is very blocky at 20x magnification because each original pixel becomes a 20x20 block of identical values. Bilinear is blurry because it linearly interpolates, creating smooth gradients that wash out details. To preserve sharp edges while upsampling, you could use edge-aware interpolation: detect edges first, then only interpolate along the edge (not across it). Bicubic interpolation also helps by using a wider neighborhood for smoother curves.

patch = img_50[:10, :10]
print("Patch shape:", patch.shape)

patch_200_nearest  = resize_nearest(patch, 200, 200)
patch_200_bilinear = resize_bilinear(patch, 200, 200)

show_images([patch, patch_200_nearest, patch_200_bilinear],
            titles=['10x10 Patch', 'Nearest (200x200)', 'Bilinear (200x200)'],
            cmap_list=['gray', 'gray', 'gray'],
            figsize=(15, 5))

1. Edge problem: At the edges, the 3x3 neighborhood extends outside the image boundaries. One solution is to only compute output where the kernel fits entirely ("valid" mode), making the output smaller. Another is to pad the image with zeros.

2. Output: The output is still a 2D grid of numbers (an image). For a 3x3 kernel on an HxW image in "valid" mode, the output is (H-2) x (W-2).

3. All weights = 1/9: This computes the average of the 3x3 neighborhood at each pixel — a box blur that smooths the image.

4. Center weight = 1, rest = 0: This copies the center pixel unchanged — the identity operation. The output is the same image (just slightly cropped in valid mode).

def conv(image, weights):
    """
    Applies a convolution filter to a grayscale image.
    image  : 2D numpy array of shape (H, W)
    weights: 2D numpy array of shape (K, K)
    Returns: 2D numpy array of shape (H-K+1, W-K+1)
    """
    H, W = image.shape
    K    = weights.shape[0]
    out_H = H - K + 1
    out_W = W - K + 1
    output = np.zeros((out_H, out_W), dtype=np.float64)

    for r in range(out_H):
        for c in range(out_W):
            total = 0.0
            for i in range(K):
                for j in range(K):
                    total += image[r + i, c + j] * weights[i, j]
            output[r, c] = total

    return output

# Sanity check: identity filter
identity_kernel = np.array([
    [0, 0, 0],
    [0, 1, 0],
    [0, 0, 0]
], dtype=np.float64)

result = conv(img_100.astype(np.float64), identity_kernel)
print("Output shape:", result.shape)        # (98, 98)
print(f"result[0,0] = {result[0,0]:.1f}")
print(f"img_100[1,1] = {img_100[1,1]}")
print("Match:", abs(result[0,0] - float(img_100[1,1])) < 1e-9)

Now that you've seen what each kernel does:

• Kernel A is a blur filter — it averages the neighborhood, smoothing the image.
• Kernel B is a sharpening filter — it enhances differences between a pixel and its neighbors.
• Kernel C is a horizontal edge detector (Sobel) — it responds to vertical changes in brightness.
• Kernel D is a vertical edge detector (Sobel) — it responds to horizontal changes in brightness.

The 5x5 blur is blurrier than the 3x3 blur because it averages over 25 pixels instead of 9 — a larger neighborhood smooths more detail. The Gaussian-like kernel produces a more natural-looking blur because it gives higher weight to nearby pixels and less to corners. Repeated blurring progressively erases fine details until the image becomes uniformly smooth.

# Task 1: Custom emboss filter
my_kernel = np.array([
    [-2, -1, 0],
    [-1,  1, 1],
    [ 0,  1, 2]
], dtype=np.float64)
result_mine = conv(img_100.astype(np.float64), my_kernel)
show_conv_result(img_100, result_mine, title='Emboss Filter')

# Task 2: 5x5 blur
kernel_blur_5x5 = np.full((5, 5), 1.0 / 25.0)
result_blur_5 = conv(img_100.astype(np.float64), kernel_blur_5x5)
show_conv_result(img_100, result_blur_5, title='5x5 Blur')

# Task 3: Gaussian-like blur
gaussian_approx = np.array([
    [1, 2, 1],
    [2, 4, 2],
    [1, 2, 1]
], dtype=np.float64)
gaussian_approx = gaussian_approx / gaussian_approx.sum()  # normalize to sum=1
result_gaussian = conv(img_100.astype(np.float64), gaussian_approx)
show_conv_result(img_100, result_gaussian, title='Gaussian-like Blur')

# Task 4: Apply blur multiple times
img_float = img_100.astype(np.float64)
kernel_A = np.full((3, 3), 1.0 / 9.0)
current = img_float
for _ in range(5):
    current = conv(current, kernel_A)

show_conv_result(img_100, current, title='5x Blur')

def pad_image(img, pad_size):
    """Pads a (H, W) image with pad_size zeros on all sides."""
    H, W = img.shape
    new_H = H + 2 * pad_size
    new_W = W + 2 * pad_size
    result = np.zeros((new_H, new_W), dtype=img.dtype)
    for r in range(H):
        for c in range(W):
            result[r + pad_size, c + pad_size] = img[r, c]
    return result

def conv_same(image, weights):
    """Convolution with 'same' padding: output is the same size as input."""
    K = weights.shape[0]
    pad_size = (K - 1) // 2
    padded = pad_image(image, pad_size)
    return conv(padded, weights)

# Test
padded = pad_image(img_100, pad_size=1)
print("Padded shape:", padded.shape)  # (102, 102)

kernel_A = np.full((3, 3), 1.0 / 9.0)
result_same = conv_same(img_100.astype(np.float64), kernel_A)
print("conv_same output shape:", result_same.shape)  # (100, 100)

You might want different kernels per channel if, for example, you want to blur the blue channel more to reduce noise (blue sensors are noisier in cameras) while keeping the green channel sharp.

def conv_color(image_rgb, weights):
    """Applies conv_same independently to each color channel."""
    H, W = image_rgb.shape[0], image_rgb.shape[1]
    result = np.zeros((H, W, 3), dtype=np.float64)
    for ch in range(3):
        channel = image_rgb[:, :, ch].astype(np.float64)
        result[:, :, ch] = conv_same(channel, weights)
    return result

# Apply blur to color image
kernel_A = np.full((3, 3), 1.0 / 9.0)
color_blurred = conv_color(color_arr, kernel_A)

# Display helper for color convolution output
def show_color_conv(original, output, title='Color Convolution'):
    lo = output.min()
    hi = output.max()
    if hi > lo:
        display = ((output - lo) / (hi - lo) * 255).astype(np.uint8)
    else:
        display = np.zeros_like(output, dtype=np.uint8)
    show_images([original, display], titles=['Original', title])

show_color_conv(color_arr, color_blurred, title='Color Blur (Kernel A)')

# Apply Sobel horizontal to color image
kernel_C = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=np.float64)
color_edges = conv_color(color_arr, kernel_C)
show_color_conv(color_arr, color_edges, title='Color Sobel Horizontal')

import math

# Step 1: Convert to grayscale (luminosity)
step1_gray = rgb_to_gray_luminosity(color_arr)

# Step 2: Resize to 50x50
step2_small = resize_bilinear(step1_gray, 50, 50)

# Step 3: Blur (Kernel A)
kernel_A = np.full((3, 3), 1.0 / 9.0)
step3_blur = conv_same(step2_small.astype(np.float64), kernel_A)

# Step 4: Edge detection (Sobel H and V, then combine)
kernel_C = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=np.float64)
kernel_D = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=np.float64)
step4_edges_h = conv_same(step3_blur, kernel_C)
step4_edges_v = conv_same(step3_blur, kernel_D)

h, w = step4_edges_h.shape
step4_edges = np.zeros((h, w), dtype=np.float64)
for r in range(h):
    for c in range(w):
        step4_edges[r, c] = math.sqrt(
            step4_edges_h[r, c] ** 2 + step4_edges_v[r, c] ** 2
        )

# Normalize edge magnitude to [0, 255]
lo, hi = step4_edges.min(), step4_edges.max()
if hi > lo:
    edges_display = ((step4_edges - lo) / (hi - lo) * 255).astype(np.uint8)
else:
    edges_display = np.zeros_like(step4_edges, dtype=np.uint8)

# Display all steps
show_images(
    [color_arr, step1_gray, step2_small,
     np.clip(step3_blur, 0, 255).astype(np.uint8), edges_display],
    titles=['1. Color', '2. Grayscale', '3. 50x50', '4. Blurred', '5. Edges'],
    cmap_list=[None, 'gray', 'gray', 'gray', 'gray']
)

# Kernel 1: Diagonal edge detector
# Positive on one diagonal, negative on the other
my_kernel_1 = np.array([
    [-1,  0,  1],
    [ 0,  0,  0],
    [ 1,  0, -1]
], dtype=np.float64)
result_1 = conv_same(img_100.astype(np.float64), my_kernel_1)
show_conv_result(img_100, result_1, title='Diagonal Edge Detector')
# Observation: highlights diagonal edges (top-left to bottom-right vs
# top-right to bottom-left)

# Kernel 2: Laplacian (all-direction edge detector)
# Negative ring around a large positive center; sums to zero
my_kernel_2 = np.array([
    [-1, -1, -1],
    [-1,  8, -1],
    [-1, -1, -1]
], dtype=np.float64)
result_2 = conv_same(img_100.astype(np.float64), my_kernel_2)
show_conv_result(img_100, result_2, title='Laplacian (All Edges)')
# Observation: detects edges in ALL directions, not just horizontal or vertical

# Kernel 3: Corner-only detector
# Only the corners have non-zero weights
my_kernel_3 = np.array([
    [ 1,  0, -1],
    [ 0,  0,  0],
    [-1,  0,  1]
], dtype=np.float64)
result_3 = conv_same(img_100.astype(np.float64), my_kernel_3)
show_conv_result(img_100, result_3, title='Corner Detector')
# Observation: responds to diagonal texture changes; highlights corners
# and diagonal patterns

Repeated sharpening amplifies noise and creates extreme values. Repeated blurring progressively erases detail until the image becomes nearly uniform. Sobel followed by blur produces smooth, thick edge lines. In a CNN: Layer 1 detects simple features (edges, gradients). Layer 2 combines edges into textures and corners. Deeper layers detect complex structures like shapes, object parts, and eventually entire objects.

# Define kernels
kernel_A = np.full((3, 3), 1.0 / 9.0)
kernel_B = np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]], dtype=np.float64)
kernel_C = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=np.float64)

# Task 1: Apply Kernel B (sharpen) 10 times
current = img_100.astype(np.float64)
snapshots_b = []
for i in range(10):
    current = conv_same(current, kernel_B)
    if i + 1 in [1, 3, 5, 10]:
        snap = current.copy()
        lo, hi = snap.min(), snap.max()
        if hi > lo:
            snap = ((snap - lo) / (hi - lo) * 255).astype(np.uint8)
        else:
            snap = np.zeros_like(snap, dtype=np.uint8)
        snapshots_b.append((i + 1, snap))

show_images([s for _, s in snapshots_b],
            titles=[f'{n}x Kernel B' for n, _ in snapshots_b],
            cmap_list=['gray'] * len(snapshots_b))

# Task 2: Apply Kernel A (blur) 20 times
current = img_100.astype(np.float64)
for _ in range(20):
    current = conv_same(current, kernel_A)
blur_20 = np.clip(current, 0, 255).astype(np.uint8)
show_images([img_100, blur_20],
            titles=['Original', '20x Blur'],
            cmap_list=['gray', 'gray'])

# Task 3: Kernel C then Kernel A (edge detect then blur)
edges = conv_same(img_100.astype(np.float64), kernel_C)
smoothed_edges = conv_same(edges, kernel_A)
show_conv_result(img_100, smoothed_edges, title='Sobel then Blur')