Computer Vision Fundamentals for Drones
Computer vision transforms raw camera feeds into spatial understanding—enabling drones to perceive, interpret, and navigate their environment. This tutorial covers the essential CV techniques for drone applications.
The Computer Vision Pipeline
Sponsored
Every drone vision system follows a structured pipeline from pixel capture to scene understanding.
Architecture Diagram
**Real-world analogy:** Computer vision for drones is like giving eyes and a brain to an aircraft. The camera captures what's there, but the CV pipeline interprets what it means—distinguishing a road from a river, a building from a tree.
Image Preprocessing
Preprocessing normalizes raw camera data for consistent analysis across varying conditions.
Real-world analogy: Like adjusting your glasses and squinting in bright sunlight—preprocessing helps the vision system "see" clearly regardless of conditions.
python
from scipy.ndimage import gaussian_filter, median_filter
class ImagePreprocessor:
"""Complete image preprocessing pipeline for drone imagery."""
def __init__(self, target_size=(640, 480)):
self.target_size = target_size
def remove_lens_distortion(self, image, k1=-0.2, k2=0.05, p1=0.01, p2=-0.01):
"""
Correct lens barrel/pincushion distortion.
In production, use calibrated camera matrices.
"""
h, w = image.shape[:2]
cx, cy = w / 2, h / 2
# Create coordinate grids
y, x = np.mgrid[0:h, 0:w].astype(np.float32)
x_norm = (x - cx) / cx
y_norm = (y - cy) / cy
# Radial distortion
r2 = x_norm**2 + y_norm**2
radial = 1 + k1 * r2 + k2 * r2**2
# Apply distortion
x_distorted = x_norm * radial + 2 * p1 * x_norm * y_norm + p2 * (r2 + 2 * x_norm**2)
y_distorted = y_norm * radial + p1 * (r2 + 2 * y_norm**2) + 2 * p2 * x_norm * y_norm
x_map = (x_distorted * cx + cx).astype(np.float32)
y_map = (y_distorted * cy + cy).astype(np.float32)
# Simple nearest-neighbor remapping (for demonstration)
result = np.zeros_like(image)
for i in range(h):
for j in range(w):
src_x = int(np.clip(x_map[i, j], 0, w-1))
src_y = int(np.clip(y_map[i, j], 0, h-1))
result[i, j] = image[src_y, src_x]
return result
def adaptive_histogram_equalization(self, image, clip_limit=2.0, grid_size=8):
"""Apply CLAHE (Contrast Limited Adaptive Histogram Equalization)."""
h, w = image.shape[:2]
grid_h, grid_w = h // grid_size, w // grid_size
result = np.zeros_like(image)
for i in range(grid_size):
for j in range(grid_size):
y_start, y_end = i * grid_h, (i + 1) * grid_h
x_start, x_end = j * grid_w, (j + 1) * grid_w
patch = image[y_start:y_end, x_start:x_end]
if len(patch.shape) == 3:
patch = np.mean(patch, axis=2).astype(np.uint8)
# Compute histogram
hist, bins = np.histogram(patch.flatten(), 256, [0, 256])
cdf = hist.cumsum()
cdf_normalized = cdf * 255 / cdf[-1]
# Equalize patch
equalized = np.interp(patch.flatten(), bins[:-1], cdf_normalized)
result[y_start:y_end, x_start:x_end] = equalized.reshape(patch.shape)
return result.astype(np.uint8)
def denoise_bilateral(self, image, sigma_space=10, sigma_color=25):
"""Bilateral filtering (edge-preserving denoising)."""
# Simplified bilateral filter
kernel_size = 5
pad = kernel_size // 2
padded = np.pad(image, ((pad, pad), (pad, pad)), mode='reflect')
result = np.zeros_like(image, dtype=np.float64)
for i in range(image.shape[0]):
for j in range(image.shape[1]):
patch = padded[i:i+kernel_size, j:j+kernel_size].astype(np.float64)
center = float(padded[i + pad, j + pad])
# Spatial weights
y_coords, x_coords = np.mgrid[0:kernel_size, 0:kernel_size]
spatial_w = np.exp(-((y_coords - pad)**2 + (x_coords - pad)**2) /
(2 * sigma_space**2))
# Range weights
range_w = np.exp(-(patch - center)**2 / (2 * sigma_color**2))
# Combined weights
weights = spatial_w * range_w
weights /= weights.sum()
result[i, j] = np.sum(patch * weights)
return result.astype(np.uint8)
def motion_deblur(self, image, kernel_size=5):
"""Simple motion deblurring using Wiener filter approximation."""
# Estimate blur kernel (simplified)
kernel = np.ones((kernel_size, kernel_size)) / (kernel_size**2)
# Frequency domain
img_fft = np.fft.fft2(image.astype(np.float64))
kernel_fft = np.fft.fft2(kernel, s=image.shape)
# Wiener deconvolution
K = 0.01 # Noise-to-signal ratio
wiener = np.conj(kernel_fft) / (np.abs(kernel_fft)**2 + K)
deblurred = np.real(np.fft.ifft2(img_fft * wiener))
return np.clip(deblurred, 0, 255).astype(np.uint8)
# Example usage
np.random.seed(42)
# Simulate drone image (grayscale for simplicity)
image = np.random.randint(50, 200, (100, 100), dtype=np.uint8)
preprocessor = ImagePreprocessor()
equalized = preprocessor.adaptive_histogram_equalization(image)
denoised = preprocessor.denoise_bilateral(image)
print(f"Original: mean={image.mean():.1f}, std={image.std():.1f}")
print(f"Equalized: mean={equalized.mean():.1f}, std={equalized.std():.1f}")
print(f"Denoised: variance reduction={image.var() - denoised.var():.1f}")
Edge Detection
Edge detection identifies boundaries—critical for obstacle avoidance and terrain analysis.
Real-world analogy: Edges are like the outlines in a coloring book. They define where one object ends and another begins, giving structure to the visual scene.
python
class EdgeDetector:
"""Edge detection algorithms for drone vision."""
def __init__(self):
# Sobel kernels
self.sobel_x = np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]], dtype=np.float64)
self.sobel_y = np.array([[-1, -2, -1],
[0, 0, 0],
[1, 2, 1]], dtype=np.float64)
def convolve2d(self, image, kernel):
"""2D convolution."""
h, w = image.shape
kh, kw = kernel.shape
pad_h, pad_w = kh // 2, kw // 2
padded = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode='reflect')
output = np.zeros_like(image, dtype=np.float64)
for i in range(h):
for j in range(w):
region = padded[i:i+kh, j:j+kw]
output[i, j] = np.sum(region * kernel)
return output
def sobel_edges(self, image):
"""Sobel edge detection."""
gray = image if len(image.shape) == 2 else np.mean(image, axis=2)
grad_x = self.convolve2d(gray, self.sobel_x)
grad_y = self.convolve2d(gray, self.sobel_y)
magnitude = np.sqrt(grad_x**2 + grad_y**2)
direction = np.arctan2(grad_y, grad_x)
return magnitude, direction
def canny_edges(self, image, low_threshold=50, high_threshold=150):
"""Canny edge detection pipeline."""
gray = image if len(image.shape) == 2 else np.mean(image, axis=2)
# 1. Gaussian smoothing
gaussian_kernel = np.array([
[1, 2, 1],
[2, 4, 2],
[1, 2, 1]
], dtype=np.float64) / 16
smoothed = self.convolve2d(gray, gaussian_kernel)
# 2. Sobel gradients
magnitude, direction = self.sobel_edges(smoothed.astype(np.uint8))
# 3. Non-maximum suppression
h, w = magnitude.shape
nms = np.zeros_like(magnitude)
for i in range(1, h-1):
for j in range(1, w-1):
angle = direction[i, j] * 180 / np.pi
angle = angle % 180
if angle < 22.5 or angle >= 157.5:
neighbors = [magnitude[i, j-1], magnitude[i, j+1]]
elif 22.5 <= angle < 67.5:
neighbors = [magnitude[i-1, j+1], magnitude[i+1, j-1]]
elif 67.5 <= angle < 112.5:
neighbors = [magnitude[i-1, j], magnitude[i+1, j]]
else:
neighbors = [magnitude[i-1, j-1], magnitude[i+1, j+1]]
if magnitude[i, j] >= max(neighbors):
nms[i, j] = magnitude[i, j]
# 4. Double threshold
strong = 255
weak = 75
result = np.zeros_like(nms)
result[nms >= high_threshold] = strong
result[(nms >= low_threshold) & (nms < high_threshold)] = weak
# 5. Edge tracking by hysteresis
for i in range(1, h-1):
for j in range(1, w-1):
if result[i, j] == weak:
if np.any(result[i-1:i+2, j-1:j+2] == strong):
result[i, j] = strong
else:
result[i, j] = 0
return result
def hough_lines(self, edges, threshold=50):
"""Simple Hough line detection."""
h, w = edges.shape
max_diag = int(np.sqrt(h**2 + w**2))
accumulator = np.zeros((180, 2 * max_diag))
# Get edge points
edge_points = np.argwhere(edges > 0)
for y, x in edge_points:
for theta in range(180):
theta_rad = theta * np.pi / 180
rho = int(x * np.cos(theta_rad) + y * np.sin(theta_rad))
rho += max_diag
accumulator[theta, rho] += 1
# Find peaks
lines = []
for theta in range(180):
for rho in range(2 * max_diag):
if accumulator[theta, rho] >= threshold:
lines.append((accumulator[theta, rho], theta, rho - max_diag))
lines.sort(reverse=True)
return lines[:20] # Return top 20 lines
# Example: Edge detection on drone imagery
np.random.seed(42)
image = np.random.randint(0, 255, (128, 128), dtype=np.uint8)
# Add some structure
image[40:50, 20:100] = 200
image[20:110, 60:70] = 180
detector = EdgeDetector()
magnitude, direction = detector.sobel_edges(image)
canny = detector.canny_edges(image, low_threshold=30, high_threshold=100)
lines = detector.hough_lines(canny, threshold=30)
print(f"Sobel magnitude range: {magnitude.min():.1f} - {magnitude.max():.1f}")
print(f"Canny edges detected: {np.sum(canny > 0)} pixels")
print(f"Hough lines found: {len(lines)}")
Feature Detection and Description
Advertisement
Feature detectors identify keypoints and describe local image patches—enabling visual odometry and SLAM.
python
class ORBDetector:
"""Simplified ORB (Oriented FAST and Rotated BRIEF) feature detector."""
def __init__(self, n_features=500):
self.n_features = n_features
def compute_gradients(self, image):
"""Compute image gradients."""
h, w = image.shape
grad_x = np.zeros_like(image, dtype=np.float64)
grad_y = np.zeros_like(image, dtype=np.float64)
grad_x[:, 1:-1] = image[:, 2:].astype(float) - image[:, :-2].astype(float)
grad_y[1:-1, :] = image[2:, :].astype(float) - image[:-2, :].astype(float)
magnitude = np.sqrt(grad_x**2 + grad_y**2)
direction = np.arctan2(grad_y, grad_x)
return magnitude, direction
def fast_corners(self, image, threshold=20):
"""Simplified FAST corner detection."""
h, w = image.shape
corners = []
# Circle of 16 pixels around center
circle_offsets = [
(0, -3), (1, -3), (2, -2), (3, -1),
(3, 0), (3, 1), (2, 2), (1, 3),
(0, 3), (-1, 3), (-2, 2), (-3, 1),
(-3, 0), (-3, -1), (-2, -2), (-1, -3)
]
for y in range(3, h - 3):
for x in range(3, w - 3):
center = int(image[y, x])
brighter = 0
darker = 0
for dx, dy in circle_offsets:
pixel = int(image[y + dy, x + dx])
if pixel > center + threshold:
brighter += 1
elif pixel < center - threshold:
darker += 1
if brighter >= 12 or darker >= 12:
corners.append((y, x))
return corners[:self.n_features]
def compute_descriptors(self, image, keypoints):
"""Compute BRIEF-like binary descriptors."""
descriptors = []
h, w = image.shape
desc_length = 32 # 256-bit descriptor
np.random.seed(42)
# Generate random test pairs
test_pairs = np.random.randint(-8, 9, (desc_length * 2, 2))
valid_keypoints = []
for y, x in keypoints:
if 16 <= y < h - 16 and 16 <= x < w - 16:
descriptor = []
for i in range(desc_length):
dy1, dx1 = test_pairs[i * 2]
dy2, dx2 = test_pairs[i * 2 + 1]
p1 = image[y + dy1, x + dx1]
p2 = image[y + dy2, x + dx2]
descriptor.append(1 if p1 > p2 else 0)
descriptors.append(np.array(descriptor))
valid_keypoints.append((y, x))
return valid_keypoints, np.array(descriptors)
def match_descriptors(self, desc1, desc2, max_distance=10):
"""Match descriptors using Hamming distance."""
matches = []
for i, d1 in enumerate(desc1):
distances = []
for j, d2 in enumerate(desc2):
distance = np.sum(d1 != d2) # Hamming distance
distances.append((distance, j))
distances.sort()
if distances[0][0] < max_distance:
if distances[0][0] < distances[1][0] * 0.75: # Lowe's ratio test
matches.append((i, distances[0][1], distances[0][0]))
return matches
# Example: Feature detection and matching
np.random.seed(42)
# Create two similar images (simulating consecutive frames)
image1 = np.random.randint(50, 200, (128, 128), dtype=np.uint8)
image2 = image1.copy()
image2 = np.roll(image2, 2, axis=1) # Simulate small translation
orb = ORBDetector(n_features=100)
keypoints1 = orb.fast_corners(image1)
keypoints2 = orb.fast_corners(image2)
kp1, desc1 = orb.compute_descriptors(image1, keypoints1)
kp2, desc2 = orb.compute_descriptors(image2, keypoints2)
matches = orb.match_descriptors(desc1, desc2)
print(f"Image 1: {len(kp1)} keypoints")
print(f"Image 2: {len(kp2)} keypoints")
print(f"Matches found: {len(matches)}")
if matches:
avg_distance = np.mean([m[2] for m in matches])
print(f"Average match distance: {avg_distance:.1f}")
Visual Odometry
Visual odometry estimates camera motion from sequential images—essential for GPS-denied navigation.
python
class VisualOdometry:
"""Monocular visual odometry for drone navigation."""
def __init__(self, focal_length=500, baseline=0.1):
self.focal_length = focal_length
self.baseline = baseline
self.prev_frame = None
self.trajectory = [(0, 0, 0)]
def detect_features(self, image, max_features=200):
"""Detect good features to track."""
h, w = image.shape
# Harris corner detection (simplified)
Ix = np.zeros_like(image, dtype=np.float64)
Iy = np.zeros_like(image, dtype=np.float64)
Ix[:, 1:-1] = image[:, 2:].astype(float) - image[:, :-2].astype(float)
Iy[1:-1, :] = image[2:, :].astype(float) - image[:-2, :].astype(float)
Ixx = Ix ** 2
Iyy = Iy ** 2
Ixy = Ix * Iy
# Gaussian window
window = np.ones((5, 5)) / 25
# Harris response
det = np.zeros_like(image, dtype=np.float64)
trace = np.zeros_like(image, dtype=np.float64)
for i in range(2, h-2):
for j in range(2, w-2):
Sxx = np.sum(Ixx[i-2:i+3, j-2:j+3] * window)
Syy = np.sum(Iyy[i-2:i+3, j-2:j+3] * window)
Sxy = np.sum(Ixy[i-2:i+3, j-2:j+3] * window)
det[i, j] = Sxx * Syy - Sxy ** 2
trace[i, j] = Sxx + Syy
response = det - 0.04 * trace ** 2
# Non-maximum suppression
features = []
for i in range(3, h-3):
for j in range(3, w-3):
if response[i, j] > 0:
local_max = np.max(response[i-3:i+4, j-3:j+4])
if response[i, j] == local_max:
features.append((i, j, response[i, j]))
features.sort(key=lambda x: x[2], reverse=True)
return [(f[0], f[1]) for f in features[:max_features]]
def track_features(self, prev_img, curr_img, prev_points, window_size=15):
"""Track features using template matching."""
tracked = []
h, w = prev_img.shape
for py, px in prev_points:
# Extract template
half_win = window_size // 2
y1 = max(0, py - half_win)
y2 = min(h, py + half_win + 1)
x1 = max(0, px - half_win)
x2 = min(w, px + half_win + 1)
template = prev_img[y1:y2, x1:x2]
if template.size == 0:
continue
# Search in current frame (assume small motion)
search_range = 20
best_score = -1
best_pos = (py, px)
for dy in range(-search_range, search_range + 1, 2):
for dx in range(-search_range, search_range + 1, 2):
ny, nx = py + dy, px + dx
if (half_win <= ny < h - half_win and
half_win <= nx < w - half_win):
patch = curr_img[ny-half_win:ny+half_win+1,
nx-half_win:nx+half_win+1]
if patch.shape == template.shape:
# Normalized cross-correlation
score = np.sum((template - template.mean()) *
(patch - patch.mean()))
denom = np.sqrt(np.sum((template - template.mean())**2) *
np.sum((patch - patch.mean())**2))
if denom > 0:
score /= denom
if score > best_score:
best_score = score
best_pos = (ny, nx)
if best_score > 0.5:
tracked.append((best_pos, (py, px)))
return tracked
def estimate_motion(self, tracked_points):
"""Estimate camera motion from feature correspondences."""
if len(tracked_points) < 4:
return np.zeros(6)
# Compute optical flow
points_curr = np.array([p[0] for p in tracked_points])
points_prev = np.array([p[1] for p in tracked_points])
flow = points_curr - points_prev
# Simple motion estimation
# Translation from mean flow
tx = np.mean(flow[:, 1]) # x-translation
ty = np.mean(flow[:, 0]) # y-translation
# Rotation from flow field
cx = np.mean(points_prev[:, 1])
cy = np.mean(points_prev[:, 0])
dx = points_prev[:, 1] - cx
dy = points_prev[:, 0] - cy
# Estimate rotation angle
angles = np.arctan2(flow[:, 0] - ty, flow[:, 1] - tx)
rz = np.median(angles)
# Scale from flow magnitude
flow_magnitude = np.sqrt(flow[:, 0]**2 + flow[:, 1]**2)
scale = np.median(flow_magnitude) / 10 if np.median(flow_magnitude) > 0 else 1
# Update trajectory
prev_x, prev_y, prev_z = self.trajectory[-1]
new_x = prev_x + tx * scale
new_y = prev_y + ty * scale
new_z = prev_z # Monocular VO doesn't give absolute scale
self.trajectory.append((new_x, new_y, new_z))
return np.array([tx, ty, 0, 0, 0, rz])
# Example: Visual odometry pipeline
np.random.seed(42)
vo = VisualOdometry(focal_length=500)
# Simulate drone frames
for frame_idx in range(5):
frame = np.random.randint(50, 200, (128, 128), dtype=np.uint8)
# Add structure
frame[40+frame_idx*2:50+frame_idx*2, 20:100] = 200
if vo.prev_frame is not None:
prev_features = vo.detect_features(vo.prev_frame)
tracked = vo.track_features(vo.prev_frame, frame, prev_features)
motion = vo.estimate_motion(tracked)
print(f"Frame {frame_idx}: "
f"Tracked {len(tracked)} features, "
f"Translation: ({motion[0]:.2f}, {motion[1]:.2f})")
vo.prev_frame = frame
print(f"\nTrajectory: {len(vo.trajectory)} poses")
print(f"Final position: {vo.trajectory[-1]}")
Hands-On Project: Drone Visual SLAM System
Build a complete visual SLAM (Simultaneous Localization and Mapping) system.
python
class SimpleSLAM:
"""Simplified visual SLAM for drone navigation."""
def __init__(self):
self.map_points = []
self.keyframes = []
self.camera_poses = []
def process_frame(self, frame, frame_id):
"""Process a new frame and update the map."""
features = self._extract_features(frame)
# Check if this is a new keyframe
if self._is_keyframe(features):
self._add_keyframe(frame, features, frame_id)
self._triangulate_points()
self._optimize()
return self.get_current_pose()
def _extract_features(self, frame):
"""Extract features from frame."""
# Simplified feature extraction
h, w = frame.shape[:2]
features = []
for _ in range(50):
y = np.random.randint(10, h-10)
x = np.random.randint(10, w-10)
descriptor = np.random.randint(0, 2, 32)
features.append({
'point': (y, x),
'descriptor': descriptor,
'frame_id': len(self.keyframes)
})
return features
def _is_keyframe(self, features):
"""Determine if current frame should be a keyframe."""
if len(self.keyframes) == 0:
return True
# Check if enough new features
return len(features) > 30
def _add_keyframe(self, frame, features, frame_id):
"""Add a new keyframe to the map."""
self.keyframes.append({
'frame': frame,
'features': features,
'frame_id': frame_id
})
# Add new map points
for feat in features:
point_3d = np.array([
feat['point'][1],
feat['point'][0],
np.random.uniform(5, 15) # Random depth
])
self.map_points.append({
'position': point_3d,
'observations': [(frame_id, feat['point'])]
})
# Camera pose (simplified)
pose = np.eye(4)
pose[0, 3] = frame_id * 0.5 # x translation
pose[2, 3] = 10 # height
self.camera_poses.append(pose)
def _triangulate_points(self):
"""Triangulate 3D points from multiple views."""
# Simplified triangulation
for point in self.map_points:
if len(point['observations']) >= 2:
# Simple midpoint triangulation
point['position'][2] = np.random.uniform(5, 15)
def _optimize(self):
"""Bundle adjustment optimization."""
# Simplified optimization
for point in self.map_points:
point['position'] += np.random.normal(0, 0.01, 3)
def get_current_pose(self):
"""Get current camera pose."""
if self.camera_poses:
return self.camera_poses[-1]
return np.eye(4)
def get_map(self):
"""Get the current map."""
return {
'points': [p['position'] for p in self.map_points],
'keyframes': len(self.keyframes),
'poses': self.camera_poses
}
# Run SLAM simulation
np.random.seed(42)
slam = SimpleSLAM()
print("Running Visual SLAM on drone footage...")
for frame_id in range(20):
frame = np.random.randint(0, 255, (240, 320), dtype=np.uint8)
pose = slam.process_frame(frame, frame_id)
if (frame_id + 1) % 5 == 0:
map_data = slam.get_map()
print(f"Frame {frame_id+1}: "
f"{map_data['keyframes']} keyframes, "
f"{len(map_data['points'])} map points, "
f"Pose: [{pose[0,3]:.1f}, {pose[1,3]:.1f}, {pose[2,3]:.1f}]")
final_map = slam.get_map()
print(f"\nSLAM Complete:")
print(f" Keyframes: {final_map['keyframes']}")
print(f" Map points: {len(final_map['points'])}")
print(f" Camera poses: {len(final_map['poses'])}")
Key Takeaways
- Preprocessing normalizes drone imagery for consistent analysis
- Edge detection identifies boundaries for obstacle avoidance
- Feature detection enables visual odometry and SLAM
- Visual odometry provides GPS-denied navigation capability
- SLAM builds maps while localizing the drone within them
Next, we'll explore object detection and tracking for identifying and following targets.