In Lesson 1 we saw that perception starts with raw sensor signals and must produce usable information. Real robot cameras generate large arrays \( I(x,y) \) of intensities (or depth, color, etc.). Feature extraction compresses these arrays into sparse, informative primitives—points or curves— that stay stable under moderate changes of viewpoint, illumination, and noise.

Formally, let an image be a function \( I:\Omega\subset\mathbb{R}^2 \rightarrow \mathbb{R} \). A feature detector is a map \( \mathcal{F}(I) = \{f_1,\dots,f_m\} \) producing features \( f_i \) such as edge pixels, corner keypoints, or blob centers.

flowchart TD
  A["Raw image I(x,y)"] --> B["Preprocess (denoise, normalize)"]
  B --> C["Compute local measures (grad, curvature, scale responses)"]
  C --> D["Detect primitives"]
  D --> E["Describe / store for later matching"]
  E --> F["Used by higher perception modules"]
        

We now build the main three local measures for primitives.

An edge is a curve along which intensity changes sharply. The central tool is the gradient. Assuming \( I \) is differentiable, define:

\[ \nabla I(x,y) = \begin{bmatrix} I_x(x,y) \\ I_y(x,y) \end{bmatrix}, \quad I_x=\frac{\partial I}{\partial x}, \quad I_y=\frac{\partial I}{\partial y}. \]

The directional derivative of \( I \) in a unit direction \( \mathbf{u}=[\cos\theta,\sin\theta]^T \) is:

\[ D_{\mathbf{u}} I = \nabla I \cdot \mathbf{u}. \]

Proposition (maximum change direction). The direction maximizing \( D_{\mathbf{u}}I \) is \( \mathbf{u}^\star = \nabla I/\|\nabla I\| \), and the maximum rate is \( \|\nabla I\| \).

Proof. By Cauchy–Schwarz:

\[ D_{\mathbf{u}} I = \nabla I \cdot \mathbf{u} \le \|\nabla I\|\;\|\mathbf{u}\| = \|\nabla I\|. \]

Equality holds iff \( \mathbf{u} \) is parallel to \( \nabla I \), hence \( \mathbf{u}^\star=\nabla I/\|\nabla I\| \). ■

Therefore, edge strength is naturally measured by the gradient magnitude: \( G(x,y)=\|\nabla I(x,y)\| \). Edge pixels are those where \( G \) is large and locally maximal along the gradient direction.

In digital images, derivatives are approximated by finite differences or filters. A simple centered difference is:

\[ I_x(x,y) \approx \frac{I(x+h,y)-I(x-h,y)}{2h}, \qquad I_y(x,y) \approx \frac{I(x,y+h)-I(x,y-h)}{2h}. \]

With pixel spacing \( h=1 \), this is a convolution with kernel \( [-\tfrac12,0,\tfrac12] \). More robust gradient estimators (Sobel, Scharr) combine differentiation with smoothing.

The celebrated Canny detector further (i) smooths with Gaussian, (ii) computes gradients, (iii) performs non-maximum suppression, (iv) uses hysteresis thresholds. We keep these steps conceptual here.

Corners occur where intensity changes strongly in two directions. The classical derivation starts from the shift–difference energy:

\[ E(\Delta x,\Delta y) = \sum_{(x,y)\in W} \big(I(x+\Delta x,y+\Delta y)-I(x,y)\big)^2, \]

where \( W \) is a small window around \( (x,y) \). For small shifts, apply a first-order Taylor expansion:

\[ I(x+\Delta x,y+\Delta y) \approx I(x,y) + I_x\Delta x + I_y\Delta y. \]

Substitute into \( E \):

\[ E(\Delta x,\Delta y) \approx \sum_{W} (I_x\Delta x + I_y\Delta y)^2 = \begin{bmatrix}\Delta x&\Delta y\end{bmatrix} \mathbf{M} \begin{bmatrix}\Delta x\\\Delta y\end{bmatrix}, \]

where the structure tensor (second-moment matrix) is

\[ \mathbf{M} = \sum_{W} \begin{bmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{bmatrix}. \]

Let eigenvalues of \( \mathbf{M} \) be \( \lambda_1,\lambda_2 \ge 0 \). Then:

• flat region: \( \lambda_1,\lambda_2 \) small
• edge: one large, one small
• corner: both large

Harris–Stephens cornerness avoids explicit eigendecomposition by using invariants:

\[ R = \det(\mathbf{M}) - k\;(\operatorname{tr}\mathbf{M})^2 = \lambda_1\lambda_2 - k(\lambda_1+\lambda_2)^2, \quad 0.04 \le k \le 0.06. \]

Corners are points where \( R \) is large and positive. This is computationally efficient and widely used in robotics.

Blobs represent compact regions brighter or darker than their surroundings. A robust way to detect them is to search for extrema in scale space.

Define a Gaussian-smoothed image at scale \( \sigma \):

\[ L(x,y;\sigma) = (G(\cdot,\cdot;\sigma) * I)(x,y), \quad G(x,y;\sigma)=\frac{1}{2\pi\sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}. \]

The Laplacian of Gaussian (LoG) response is:

\[ \nabla^2 L = \frac{\partial^2 L}{\partial x^2} + \frac{\partial^2 L}{\partial y^2}. \]

To compare responses across scales, use the scale-normalized LoG:

\[ \mathcal{B}(x,y;\sigma) = \sigma^2 \nabla^2 L(x,y;\sigma). \]

Scale-normalization idea. If an image is resized by factor \( a \), the Laplacian scales as \( \nabla^2 L_a(x,y;\sigma) = a^{-2}\nabla^2 L(ax,ay;a\sigma) \). Multiplying by \( \sigma^2 \) compensates this:

\[ \sigma^2 \nabla^2 L_a(x,y;\sigma) = (a\sigma)^2 \nabla^2 L(ax,ay;a\sigma). \]

Thus extrema of \( \mathcal{B} \) are (approximately) invariant to scale. In practice, a Difference-of-Gaussians (DoG) approximates LoG:

\[ \text{DoG}(x,y;\sigma) = L(x,y;k\sigma) - L(x,y;\sigma) \approx (k-1)\sigma^2 \nabla^2 L(x,y;\sigma), \quad k \approx 1.6. \]

flowchart LR
  I0["I(x,y)"] --> GE["Edges: large |grad I|"]
  I0 --> HC["Corners: large det(M), small trace(M)"]
  I0 --> BL["Blobs: extrema of sigma^2 * Laplacian(Gauss * I)"]
  GE --> Use1["Boundaries / outlines"]
  HC --> Use2["Repeatable keypoints"]
  BL --> Use3["Region-like landmarks"]
        

Below are simple examples using standard robotics/vision libraries. We also include a short from-scratch edge detector in Python to connect to the math.

5.1 Python (OpenCV + NumPy)

import cv2
import numpy as np

img = cv2.imread("scene.png", cv2.IMREAD_GRAYSCALE)

# --- Edges (Sobel magnitude) ---
Ix = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
Iy = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)
G  = np.sqrt(Ix**2 + Iy**2)
edges = (G > 50).astype(np.uint8) * 255

# --- Corners (Harris) ---
harris = cv2.cornerHarris(np.float32(img), blockSize=2, ksize=3, k=0.04)
corners = np.argwhere(harris > 0.01 * harris.max())

# --- Blobs (DoG / SIFT-like detector, but only detection) ---
blob_detector = cv2.SimpleBlobDetector_create()
keypoints = blob_detector.detect(img)

print("corners:", len(corners), "blobs:", len(keypoints))
      

5.2 Python from scratch (finite-difference edges)

import numpy as np

def edge_map_finite_difference(I, thresh=30.0):
    # I: 2D float image
    Kx = np.array([[-0.5, 0.0, 0.5]])
    Ky = Kx.T
    # convolution (naive)
    Ix = np.zeros_like(I)
    Iy = np.zeros_like(I)
    for y in range(1, I.shape[0]-1):
        for x in range(1, I.shape[1]-1):
            patch = I[y-1:y+2, x-1:x+2]
            Ix[y,x] = np.sum(patch * Kx)
            Iy[y,x] = np.sum(patch * Ky)
    G = np.sqrt(Ix*Ix + Iy*Iy)
    return (G > thresh).astype(np.uint8)

# example:
# edges = edge_map_finite_difference(img.astype(float))
      

5.3 C++ (OpenCV)

#include <opencv2/opencv.hpp>
#include <iostream>
int main(){
    cv::Mat img = cv::imread("scene.png", cv::IMREAD_GRAYSCALE);

    // Edges
    cv::Mat Ix, Iy, G;
    cv::Sobel(img, Ix, CV_64F, 1, 0, 3);
    cv::Sobel(img, Iy, CV_64F, 0, 1, 3);
    cv::magnitude(Ix, Iy, G);
    cv::Mat edges = G > 50;

    // Corners
    cv::Mat harris;
    cv::cornerHarris(img, harris, 2, 3, 0.04);
    std::vector<cv::Point> corners;
    for(int y=0;y<harris.rows;y++)
      for(int x=0;x<harris.cols;x++)
        if(harris.at<float>(y,x) > 0.01f*harris.at<float>(0,0))
          corners.push_back(cv::Point(x,y));

    // Blobs
    auto detector = cv::SimpleBlobDetector::create();
    std::vector<cv::KeyPoint> keypoints;
    detector->detect(img, keypoints);

    std::cout << "corners=" << corners.size()
              << " blobs=" << keypoints.size() << std::endl;
    return 0;
}
      

5.4 Java (OpenCV)

import org.opencv.core.*;
import org.opencv.imgcodecs.Imgcodecs;
import org.opencv.imgproc.Imgproc;
import org.opencv.features2d.*;

public class FeaturesDemo {
  static { System.loadLibrary(Core.NATIVE_LIBRARY_NAME); }

  public static void main(String[] args){
    Mat img = Imgcodecs.imread("scene.png", Imgcodecs.IMREAD_GRAYSCALE);

    // Edges
    Mat Ix = new Mat(), Iy = new Mat(), G = new Mat();
    Imgproc.Sobel(img, Ix, CvType.CV_64F, 1, 0);
    Imgproc.Sobel(img, Iy, CvType.CV_64F, 0, 1);
    Core.magnitude(Ix, Iy, G);
    Mat edges = new Mat();
    Imgproc.threshold(G, edges, 50, 255, Imgproc.THRESH_BINARY);

    // Corners (Harris)
    Mat harris = new Mat();
    Imgproc.cornerHarris(img, harris, 2, 3, 0.04);

    // Blobs
    SimpleBlobDetector detector = SimpleBlobDetector.create();
    MatOfKeyPoint keypoints = new MatOfKeyPoint();
    detector.detect(img, keypoints);

    System.out.println("blobs=" + keypoints.toArray().length);
  }
}
      

5.5 MATLAB / Simulink

I = im2double(imread('scene.png'));

% Edges (Sobel)
[Gx, Gy] = imgradientxy(I,'sobel');
G = hypot(Gx, Gy);
edges = G > 0.2;

% Corners (Harris)
corners = detectHarrisFeatures(I);

% Blobs (LoG)
blobs = detectSURFFeatures(I); % uses scale-space extrema idea

fprintf('corners=%d blobs=%d\n', corners.Count, blobs.Count);
      

Simulink note: a typical block chain is Image From File → Gaussian Filter → Sobel Filter / Harris Corner Detector / Blob Analysis → Display. You can assemble these using the Computer Vision Toolbox blocks.

Problem 1 (Directional derivative maximum): Show that for any unit vector \( \mathbf{u} \), \( D_{\mathbf{u}}I \le \|\nabla I\| \), and find the maximizing direction.

Solution: Using Cauchy–Schwarz,

\[ D_{\mathbf{u}}I = \nabla I \cdot \mathbf{u} \le \|\nabla I\|\,\|\mathbf{u}\|=\|\nabla I\|. \]

Equality holds iff \( \mathbf{u} \parallel \nabla I \). Thus \( \mathbf{u}^\star=\nabla I/\|\nabla I\| \).

Problem 2 (Harris corner derivation): Starting from \( E(\Delta x,\Delta y)=\sum_W (I(x+\Delta x,y+\Delta y)-I(x,y))^2 \), derive the quadratic form with matrix \( \mathbf{M} \).

Solution: First-order Taylor:

\[ I(x+\Delta x,y+\Delta y)-I(x,y) \approx I_x\Delta x + I_y\Delta y. \]

Then

\[ E(\Delta x,\Delta y) \approx \sum_W (I_x\Delta x + I_y\Delta y)^2 = \begin{bmatrix}\Delta x&\Delta y\end{bmatrix} \left(\sum_W \begin{bmatrix} I_x^2 & I_xI_y\\ I_xI_y & I_y^2 \end{bmatrix}\right) \begin{bmatrix}\Delta x\\\Delta y\end{bmatrix}. \]

Hence \( \mathbf{M}=\sum_W \begin{bmatrix}I_x^2&I_xI_y\\I_xI_y&I_y^2\end{bmatrix} \).

Problem 3 (Edge vs. corner eigenvalues): Let \( \mathbf{M} \) have eigenvalues \( \lambda_1\ge\lambda_2 \ge 0 \). Explain why edges correspond to \( \lambda_1\gg\lambda_2 \) and corners to \( \lambda_1\approx\lambda_2\gg 0 \).

Solution: The energy for a shift is \( E(\Delta)=\Delta^T\mathbf{M}\Delta \). In the eigenbasis, this is \( E = \lambda_1\alpha_1^2 + \lambda_2\alpha_2^2 \). If only one eigenvalue is large, shifts along the small-eigenvalue direction cause little change, matching a 1D edge. If both are large, shifts in any direction increase energy, implying a 2D corner-like junction.

Problem 4 (DoG approximates LoG): Use a first-order expansion in scale to show \( L(x,y;k\sigma)-L(x,y;\sigma) \approx (k-1)\sigma^2\nabla^2 L(x,y;\sigma) \).

Solution: Consider \( L(\sigma) \) as a function of \( \sigma \). For small \( \delta\sigma=(k-1)\sigma \),

\[ L(k\sigma) \approx L(\sigma) + \delta\sigma\frac{\partial L}{\partial \sigma}. \]

In scale-space theory, \( \frac{\partial L}{\partial \sigma} = \sigma \nabla^2 L \). Therefore,

\[ L(k\sigma)-L(\sigma) \approx (k-1)\sigma\cdot (\sigma\nabla^2 L) = (k-1)\sigma^2\nabla^2 L. \]

So DoG is proportional to the normalized LoG response.

We introduced feature extraction as a key step turning raw images into sparse, meaningful primitives for robotics. Edges arise from large gradient magnitudes, corners from large two-directional intensity variation captured by the structure tensor, and blobs from extrema of scale-normalized Laplacian responses. These families form the backbone of many robot vision systems, even when later combined with more advanced pipelines.