**Real-world analogy:** Semantic segmentation is like a city planner coloring a map—every inch of land gets classified as residential, commercial, park, or road. Nothing is left unclassified.

## U-Net Architecture

U-Net's symmetric encoder-decoder with skip connections is the gold standard for drone image segmentation.

```python

class ConvBlock:
    """Double convolution block for U-Net."""

    def __init__(self, in_channels, out_channels):
        self.weights1 = np.random.randn(out_channels, in_channels, 3, 3) * np.sqrt(2.0 / (in_channels * 9))
        self.bias1 = np.zeros(out_channels)
        self.weights2 = np.random.randn(out_channels, out_channels, 3, 3) * np.sqrt(2.0 / (out_channels * 9))
        self.bias2 = np.zeros(out_channels)

    def conv2d(self, x, weights, bias, stride=1, padding=1):
        """2D convolution."""
        batch, in_ch, h, w = x.shape
        out_ch = weights.shape[0]
        k = weights.shape[2]

        # Padding
        x_padded = np.pad(x, ((0,0), (0,0), (padding,padding), (padding,padding)), mode='reflect')

        out_h = (h + 2*padding - k) // stride + 1
        out_w = (w + 2*padding - k) // stride + 1
        output = np.zeros((batch, out_ch, out_h, out_w))

        for b in range(batch):
            for oc in range(out_ch):
                for i in range(out_h):
                    for j in range(out_w):
                        h_start = i * stride
                        h_end = h_start + k
                        w_start = j * stride
                        w_end = w_start + k
                        output[b, oc, i, j] = np.sum(x_padded[b, :, h_start:h_end, w_start:w_end] * weights[oc]) + bias[oc]

        return output

    def relu(self, x):
        return np.maximum(0, x)

    def forward(self, x):
        x = self.conv2d(x, self.weights1, self.bias1)
        x = self.relu(x)
        x = self.conv2d(x, self.weights2, self.bias2)
        x = self.relu(x)
        return x

class UNet:
    """U-Net architecture for semantic segmentation."""

    def __init__(self, in_channels=3, num_classes=6):
        # Encoder
        self.enc1 = ConvBlock(in_channels, 64)
        self.enc2 = ConvBlock(64, 128)
        self.enc3 = ConvBlock(128, 256)
        self.enc4 = ConvBlock(256, 512)

        # Bottleneck
        self.bottleneck = ConvBlock(512, 1024)

        # Decoder
        self.dec4 = ConvBlock(1024 + 512, 512)
        self.dec3 = ConvBlock(512 + 256, 256)
        self.dec2 = ConvBlock(256 + 128, 128)
        self.dec1 = ConvBlock(128 + 64, 64)

        # Output
        self.output_conv = np.random.randn(num_classes, 64, 1, 1) * 0.01
        self.output_bias = np.zeros(num_classes)

    def max_pool(self, x, pool_size=2):
        """Max pooling."""
        batch, ch, h, w = x.shape
        out_h = h // pool_size
        out_w = w // pool_size
        output = np.zeros((batch, ch, out_h, out_w))

        for b in range(batch):
            for c in range(ch):
                for i in range(out_h):
                    for j in range(out_w):
                        region = x[b, c, i*pool_size:(i+1)*pool_size,
                                  j*pool_size:(j+1)*pool_size]
                        output[b, c, i, j] = np.max(region)
        return output

    def upsample(self, x, target_size):
        """Bilinear upsampling."""
        batch, ch, h, w = x.shape
        output = np.zeros((batch, ch, target_size, target_size))

        for b in range(batch):
            for c in range(ch):
                for i in range(target_size):
                    for j in range(target_size):
                        src_i = i * h / target_size
                        src_j = j * w / target_size
                        i0, j0 = int(src_i), int(src_j)
                        i1, j1 = min(i0+1, h-1), min(j0+1, w-1)
                        di, dj = src_i - i0, src_j - j0

                        output[b, c, i, j] = (
                            x[b, c, i0, j0] * (1-di) * (1-dj) +
                            x[b, c, i1, j0] * di * (1-dj) +
                            x[b, c, i0, j1] * (1-di) * dj +
                            x[b, c, i1, j1] * di * dj
                        )
        return output

    def forward(self, x):
        """Forward pass through U-Net."""
        # Encoder
        e1 = self.enc1.forward(x)
        e2 = self.enc2.forward(self.max_pool(e1))
        e3 = self.enc3.forward(self.max_pool(e2))
        e4 = self.enc4.forward(self.max_pool(e3))

        # Bottleneck
        b = self.bottleneck.forward(self.max_pool(e4))

        # Decoder with skip connections
        d4 = self.upsample(b, e4.shape[2])
        d4 = np.concatenate([d4, e4], axis=1)
        d4 = self.dec4.forward(d4)

        d3 = self.upsample(d4, e3.shape[2])
        d3 = np.concatenate([d3, e3], axis=1)
        d3 = self.dec3.forward(d3)

        d2 = self.upsample(d3, e2.shape[2])
        d2 = np.concatenate([d2, e2], axis=1)
        d2 = self.dec2.forward(d2)

        d1 = self.upsample(d2, e1.shape[2])
        d1 = np.concatenate([d1, e1], axis=1)
        d1 = self.dec1.forward(d1)

        # Output
        out = np.zeros((x.shape[0], self.output_conv.shape[0], x.shape[2], x.shape[3]))
        for b in range(x.shape[0]):
            for c in range(self.output_conv.shape[0]):
                out[b, c] = np.sum(d1[b] * self.output_conv[c].squeeze(), axis=0) + self.output_bias[c]

        # Softmax
        out_exp = np.exp(out - np.max(out, axis=1, keepdims=True))
        out = out_exp / np.sum(out_exp, axis=1, keepdims=True)

        return out

# Example: Run U-Net inference
np.random.seed(42)
model = UNet(in_channels=3, num_classes=6)

# Simulate drone image batch
batch = np.random.randn(2, 3, 64, 64)  # Smaller for demo
output = model.forward(batch)

print(f"Input shape: {batch.shape}")
print(f"Output shape: {output.shape}")
print(f"Classes: Building, Road, Vegetation, Water, Vehicle, Person")

# Per-pixel predictions
predictions = np.argmax(output, axis=1)
for i in range(2):
    unique, counts = np.unique(predictions[i], return_counts=True)
    print(f"\nImage {i+1} class distribution:")
    for cls, count in zip(unique, counts):
        pct = count / predictions[i].size * 100
        print(f"  Class {cls}: {count} pixels ({pct:.1f}%)")