Sample program to query an aerial drone image:

# filename: vlm_scene_query.py

"""

Download an aerial image from an Azure SAS URL, extract XFIF/GPS metadata,

and query a vision-language model (e.g., Qwen2.5VL-7B VLM) to answer a

scene-level question such as estimating area based on parking spot counts.

Usage:

    export MODEL_ID="your-qwen-model-id-or-hf-repo"

    export HF_API_TOKEN="..." # if required by the model host

    python vlm_scene_query.py --sas-url "<SAS_URL>" --question "Estimate area in square meters"

Notes:

- This script uses the Hugging Face transformers pipeline as a generic interface.

  Qwen VLM may require a provider-specific SDK or a different pipeline name.

  Replace the model loading section with the provider-specific code if needed.

- The script extracts GPS EXIF if present and returns it with the model response.

"""

import os

import sys

import argparse

import tempfile

import json

import math

import logging

from typing import Optional, Dict, Any, Tuple

import requests

from PIL import Image

import exifread

# Optional: transformers pipeline for vision-language models

try:

    from transformers import pipeline, AutoTokenizer, AutoModelForSeq2SeqLM

    HF_AVAILABLE = True

except Exception:

    HF_AVAILABLE = False

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")

logger = logging.getLogger("vlm_scene_query")

def download_image_from_sas(sas_url: str, dest_path: str, timeout: int = 30) -> None:

    """Download an image from an Azure SAS URL to dest_path."""

    logger.info("Downloading image from SAS URL")

    resp = requests.get(sas_url, stream=True, timeout=timeout)

    resp.raise_for_status()

    with open(dest_path, "wb") as f:

        for chunk in resp.iter_content(chunk_size=8192):

            if chunk:

                f.write(chunk)

    logger.info("Downloaded image to %s", dest_path)

def extract_gps_from_exif(image_path: str) -> Dict[str, Any]:

    """Extract GPS EXIF data (if present) using exifread and return a dict."""

    logger.info("Extracting EXIF metadata")

    with open(image_path, "rb") as f:

        tags = exifread.process_file(f, details=False)

    gps = {}

    def _get(tag):

        return tags.get(tag)

    # Common EXIF GPS tags

    lat_ref = _get("GPS GPSLatitudeRef")

    lat = _get("GPS GPSLatitude")

    lon_ref = _get("GPS GPSLongitudeRef")

    lon = _get("GPS GPSLongitude")

    alt = _get("GPS GPSAltitude")

    if lat and lon and lat_ref and lon_ref:

        def _to_deg(value):

            # value is like [Rational(37,1), Rational(46,1), Rational(0,1)]

            try:

                parts = [float(x.num) / float(x.den) for x in value.values]

                deg = parts[0] + parts[1] / 60.0 + parts[2] / 3600.0

                return deg

            except Exception:

                return None

        lat_deg = _to_deg(lat)

        lon_deg = _to_deg(lon)

        if lat_deg is not None and lon_deg is not None:

            if str(lat_ref).upper().startswith("S"):

                lat_deg = -lat_deg

            if str(lon_ref).upper().startswith("W"):

                lon_deg = -lon_deg

            gps["latitude"] = lat_deg

            gps["longitude"] = lon_deg

    if alt:

        try:

            gps["altitude"] = float(alt.values[0].num) / float(alt.values[0].den)

        except Exception:

            pass

    # XFIF or other tags may be present; include raw tags for inspection

    gps["raw_tags"] = {k: str(v) for k, v in tags.items() if k.startswith("GPS")}

    return gps

def default_area_estimate_from_parking_count(parking_count: int,

                                             spot_length_m: float = 4.5,

                                             spot_width_m: float = 1.8,

                                             spacing_factor: float = 1.2) -> Tuple[float, float]:

    """

    Estimate area in square meters and square feet given a count of parking spots.

    Default sedan footprint: 4.5m x 1.8m = 8.1 m^2. spacing_factor accounts for drive lanes and spacing.

    Returns (area_m2, area_ft2).

    """

    single_spot_area = spot_length_m * spot_width_m * spacing_factor

    total_m2 = parking_count * single_spot_area

    total_ft2 = total_m2 * 10.7639

    return total_m2, total_ft2

def build_prompt_for_vlm(question: str, guidance: Optional[str] = None) -> str:

    """

    Build a clear prompt for the vision-language model. Guidance can include

    assumptions to make (e.g., sedan footprint).

    """

    base = (

        "You are given an aerial image. Answer the user's question precisely. "

        "If you need to make reasonable assumptions, state them explicitly. "

        "Return a JSON object with keys: 'answer_text', 'parking_spot_count' (int or null), "

        "'assumptions' (list of strings), and 'computed' (object with numeric fields). "

    )

    if guidance:

        base += guidance + " "

    base += "User question: " + question

    return base

def query_vlm_with_image(model_id: str, image_path: str, prompt: str, hf_token: Optional[str] = None) -> Dict[str, Any]:

    import torch

    from PIL import Image

    from transformers import AutoProcessor, AutoModelForCausalLM

    if hf_token:

        os.environ["HUGGINGFACEHUB_API_TOKEN"] = hf_token

    device = "cuda" if torch.cuda.is_available() else "cpu"

    processor = AutoProcessor.from_pretrained(model_id, trust_remote_code=True)

    model = AutoModelForCausalLM.from_pretrained(model_id, trust_remote_code=True).to(device)

    image = Image.open(image_path).convert("RGB")

    # Qwen expects both text + image in the processor

    inputs = processor(text=prompt, images=image, return_tensors="pt").to(device)

    generated_ids = model.generate(

        **inputs,

        max_new_tokens=512

    )

    output_text = processor.batch_decode(generated_ids, skip_special_tokens=True)[0]

    try:

        return json.loads(output_text)

    except Exception:

        return {"answer_text": output_text, "raw": output_text}

def torch_cuda_available() -> bool:

    try:

        import torch

        return torch.cuda.is_available()

    except Exception:

        return False

def parse_parking_count_from_vlm_response(vlm_resp: Dict[str, Any]) -> Optional[int]:

    """

    Extract parking_spot_count if present in the VLM response dict.

    """

    try:

        count = vlm_resp.get("parking_spot_count")

        if count is None:

            # Try to parse from answer_text heuristically

            text = vlm_resp.get("answer_text", "")

            # naive heuristic: find first integer in text

            import re

            m = re.search(r"\b(\d{1,4})\b", text)

            if m:

                return int(m.group(1))

            return None

        return int(count)

    except Exception:

        return None

def main():

    parser = argparse.ArgumentParser(description="Query a VLM about an aerial image from an Azure SAS URL")

    parser.add_argument("--sas-url", required=True, help="Azure SAS URL to the JPEG image")

    parser.add_argument("--question", required=True, help="Natural language question to ask the VLM")

    parser.add_argument("--model-id", default=os.environ.get("MODEL_ID", "Qwen/Qwen-2.5V-L-7B"), help="VLM model id or repo")

    parser.add_argument("--hf-token", default=os.environ.get("HF_API_TOKEN"), help="Hugging Face API token if required")

    parser.add_argument("--assume-spot-length-m", type=float, default=4.5, help="Assumed parking spot length in meters")

    parser.add_argument("--assume-spot-width-m", type=float, default=1.8, help="Assumed parking spot width in meters")

    parser.add_argument("--spacing-factor", type=float, default=1.2, help="Factor to account for drive lanes and spacing")

    parser.add_argument("--no-vlm", action="store_true", help="Skip VLM and use deterministic heuristic only")

    args = parser.parse_args()

    with tempfile.TemporaryDirectory() as tmpdir:

        img_path = os.path.join(tmpdir, "scene.jpg")

        try:

            download_image_from_sas(args.sas_url, img_path)

        except Exception as e:

            logger.error("Failed to download image: %s", e)

            sys.exit(1)

        gps = extract_gps_from_exif(img_path)

        logger.info("Extracted GPS metadata: %s", gps)

        # Build prompt

        guidance = (

            f"Assume a typical U.S. sedan footprint of {args.assume_spot_length_m}m x {args.assume_spot_width_m}m "

            f"and a spacing factor of {args.spacing_factor} to account for drive lanes. "

            "Count visible parking spots if possible and compute total area in square meters and square feet."

        )

        prompt = build_prompt_for_vlm(args.question, guidance=guidance)

        vlm_response = None

        parking_count = None

        if not args.no_vlm:

            try:

                vlm_response = query_vlm_with_image(args.model_id, img_path, prompt, hf_token=args.hf_token)

                logger.info("VLM response received")

                parking_count = parse_parking_count_from_vlm_response(vlm_response)

            except Exception as e:

                logger.warning("VLM query failed or not available: %s", e)

                vlm_response = {"error": str(e)}

                parking_count = None

        # If VLM didn't provide a parking count, fall back to asking user assumption or using a heuristic

        if parking_count is None:

            # Heuristic fallback: try to detect cars using a very small, dependency-free heuristic is not reliable.

            # Instead, we will ask the model's textual output for a number if available; otherwise, default to 10 spots.

            if vlm_response and isinstance(vlm_response, dict):

                parking_count = parse_parking_count_from_vlm_response(vlm_response)

            if parking_count is None:

                logger.info("No parking count from VLM; using fallback default of 10 spots for estimation")

                parking_count = 10 # conservative default; in production, prefer human-in-the-loop

        area_m2, area_ft2 = default_area_estimate_from_parking_count(

            parking_count,

            spot_length_m=args.assume_spot_length_m,

            spot_width_m=args.assume_spot_width_m,

            spacing_factor=args.spacing_factor

        )

        # Build final structured response

        response = {

            "run_id": f"run-{os.urandom(6).hex()}",

            "agent_id": "qwen-vlm-estimator",

            "start_time": None,

            "end_time": None,

            "model_version": args.model_id,

            "gps": gps,

            "question": args.question,

            "vlm_raw_response": vlm_response,

            "parking_spot_count_used": parking_count,

            "assumptions": [

                f"sedan footprint {args.assume_spot_length_m}m x {args.assume_spot_width_m}m",

                f"spacing factor {args.spacing_factor}"

            ],

            "computed": {

                "area_m2": round(area_m2, 2),

                "area_ft2": round(area_ft2, 2),

                "spot_area_m2": round(args.assume_spot_length_m * args.assume_spot_width_m * args.spacing_factor, 2)

            },

            "answer_text": (

                f"Estimated total area ≈ {round(area_m2,2)} m² ({round(area_ft2,2)} ft²) "

                f"based on {parking_count} parking spots and assumed sedan footprint "

                f"{args.assume_spot_length_m}m x {args.assume_spot_width_m}m with spacing factor {args.spacing_factor}."

            )

        }

        # Print JSON response

        print(json.dumps(response, indent=2))

if __name__ == "__main__":

    main()

Results:

{

  "run_id": "run-186a8c43b7e7",

  "agent_id": "qwen-vlm-estimator",

  "start_time": null,

  "end_time": null,

  "model_version": "Qwen/Qwen2.5-VL-7B-Instruct",

  "gps": {

    "raw_tags": {}

  },

  "question": "Estimate area of unoccupied spots in square meters",

  "parking_spot_count_used": 10,

  "assumptions": [

    "sedan footprint 4.5m x 1.8m",

    "spacing factor 1.2"

  ],

  "computed": {

    "area_m2": 97.2,

    "area_ft2": 1046.25,

    "spot_area_m2": 9.72

  },

  "answer_text": "Estimated total area \u2248 97.2 m\u00b2 (1046.25 ft\u00b2) based on 10 parking spots and assumed sedan footprint 4.5m x 1.8m with spacing factor 1.2."

}

 Agentic judges for drone image analytics

Andrew Ng’s agentic workflow pattern—reflection, tool use, planning, and multi-agent collaboration—applies to drone-vision benchmarking. Reflection lets an agent critique and revise detections, captions, SQL answers, or mission reports. Tool use grounds reasoning in retrieval, code execution, geospatial operators, detector APIs, and database queries. Planning decomposes a workload into explicit steps before execution and supports replanning when a probe fails. Multi-agent orchestration assigns specialized roles: one agent checks geospatial consistency, another evaluates temporal coherence, another audits semantic alignment, and an arbiter aggregates evidence into a score.

Memory has been a design constraint. Loops let agents think; graphs let agents remember. A reflection loop can improve a single answer, but without persistent state the agent forgets why it inspected a frame, which detector disagreed, or which workload constraint failed. A graph turns those transient observations into reusable memory: frames, objects, captions, detections, SQL results, tool calls, critiques, and final judgments become linked evidence rather than buried transcript text. In ezbenchmark, this converts an agentic judge from a one-pass evaluator into a stateful audit system.

A practical build path is incremental. First, add one critique call after every generated answer or score; this is the highest-return change because it catches unsupported claims, missing visual evidence, and weak workload alignment before output is finalized. Second, expose the judge to tools: vector retrieval over frame evidence, SQL over the scene catalog, detector re-runs, spatial predicates, and temporal-neighbor comparisons. Third, require the judge to emit a structured plan before execution, such as JSON steps with expected evidence, tool calls, and failure conditions. Fourth, split evaluation across specialized agents and connect them through a shared graph store. The result is not merely a stronger prompt; it is an architecture in which weak models can outperform stronger single-pass models because the workflow supplies iteration, grounding, decomposition, and durable memory.

Agents are usually dedicated to perception, reasoning, and control in different ways. Sapkota et al. introduce the term “Agentic UAVs” to describe systems that integrate perception, cognition, control, and communication into layered, goal-driven agents that operate with contextual reasoning and memory, rather than fixed scripts or reactive control loops [1]. In their framework, aerial image understanding is only one layer in a broader cognitive stack: perception agents extract structure from imagery and other sensors; cognitive agents plan and replan missions; control agents execute trajectories; and communication agents coordinate with humans and other UAVs. This layered view is useful when we start thinking about agentic frameworks as “judges” for benchmarking: the judging capability can itself be an agent, sitting in the cognition layer, consuming outputs from perception agents and workload metadata rather than raw pixels alone [1].

Vision–language–driven agents are a distinct subclass. Sapkota et al. explicitly highlight vision–language models and multimodal sensing as key enabling technologies for Agentic UAVs, noting that they allow agents to parse complex scenes, follow natural-language instructions, and ground symbolic goals in visual context [1]. These agents differ from traditional planners in that they can reason over image and text jointly, which makes them natural candidates for roles like “mission explainer,” “anomaly triager,” or, in our case, “benchmark judge” for aerial analytics workloads. Instead of judging purely from numeric metrics, a vision–language agent can look at a drone scene, read a workload description, inspect candidate outputs, and form a qualitative judgment about which pipeline better captures the intended analytic semantics [1].

 Drone Fleet motion propagation

Autonomous drone fleet does not have a centralized controller. Motion of the forward section of the fleet must be followed by the backward section of the fleet via propagation of the direction and speed information. This is like linear heat propagation and is best described by Fourier Transforms. Study of the propagation of fleet movement information is necessary to enable the fleet to behave as coordinated as a whole unit. The movement information is mere payload as it can be enriched with many forms of sensor data that can help the fleet to perform actions such as avoiding obstacles, optimizing flight paths and speed, and achieving formation objectives or cost functions. While these actions can be computed locally or in co-ordination by one or more members of the fleet, the propagation of messages is independent of the computation and involves phenomena like heat propagation. This article discusses the wave propagation by Fourier Transforms assuming peer processing of sensor data can be coordinated with consensus algorithms.

A Fast Fourier Transform converts wave form data in the time domain into the frequency domain. It achieves this by breaking down the original time-based waveform into a series of sinusoidal terms, each with a unique magnitude, frequency, and phase. This process converts a waveform in the time domain into a series of sinusoidal functions which when added together reconstruct the original waveform. Plotting the amplitude of each sinusoidal term versus its frequency creates a power spectrum, which is the response of the original waveform in the frequency domain.

The Fourier transform is a generalization of the Fourier series and allows taking any function as a total of simple sinusoids. A functions’ Fourier transform is a complex-valued function denoting the constituent complex sinusoid’s that contain the original function. For every frequency, the magnitude of the complex value denotes the constituent complex sinusoid’s amplitude with that frequency, and the complex value’s argument denotes the phase offset of the complex sinusoid. If a frequency does not exist, the transform possesses a value of zero for that frequency. Functions that are localized in the domain of time have Fourier transforms that extend out across the domain of frequency and vice versa. The Fourier transform of a Gaussian function is always another Gaussian function. The solutions to heat equations are the Gaussian functions.

The Fourier transform can be generalized to functions of various variables on Euclidean space, forwarding a function of three-dimensional position space to a three-dimensional momentum function (or a space and time function to a 4-momentum function). This helps spatial Fourier Transform to be solutions of waves as functions of either momentum or position or both.

When Fourier transforms are applicable, it means the “earth response” now is the same as the “earth response” later. Switching our point of view from time to space, the applicability of the Fourier transformation means that the “impulse response” here is the same as the “impulse response” there. An impulse is a column vector full of zeros with somewhere a one. An impulse response is a column from the matrix q = Bp The collection of impulse responses in q=Bp defines the convolution operation.

The difference between Fourier Transform and Fourier Series is that the Fourier Transform is applicable for non-periodic signals, while the Fourier Series is applicable to periodic signals. The properties of Fourier transform are duality, linear transform, modulation and Parseval’s theorem. Duality implies that if h(t) possesses a Fourier transform H(f), then the Fourier transform related to H(t) is H(-f). Linear Transform implies that if g(t) and h(t) are two Fourier transforms also represented by G(f) and H(f) respectively, then the linear combination of h and g also has a Fourier transform. The modulation property implies that functions are modulated by other functions if they are multiplied in time. Parseval’s theorem states that the Fourier transform is unitary and the sum of the squares of the H(f) equals that of h(t)

While the dominant interest in the application of wave propagation transforms to fleet movement is the duration it takes for the entire fleet to respond, the drone traffic can be taken as a Gaussian distribution and one that can even be treated as an invariant through different formations. This makes it easier to predict the fleet movements.

Sample FFT application:

import numpy as nm

import scipy

import scipy.fftpack

import pylab

def lowpass_cosine( y, tau, f_3db, width, padd_data=True):

    # padd_data = True means we are going to symmetric copies of the data to the start and stop

    # to reduce/eliminate the discontinuities at the start and stop of a dataset due to filtering

    #

    # False means we're going to have transients at the start and stop of the data

    # kill the last data point if y has an odd length

    if nm.mod(len(y),2):

        y = y[0:-1]

    # add the weird padd

    # so, make a backwards copy of the data, then the data, then another backwards copy of the data

    if padd_data:

        y = nm.append( nm.append(nm.flipud(y),y) , nm.flipud(y) )

    # take the FFT

    ffty=scipy.fftpack.fft(y)

    ffty=scipy.fftpack.fftshift(ffty)

    # make the companion frequency array

    delta = 1.0/(len(y)*tau)

    nyquist = 1.0/(2.0*tau)

    freq = nm.arange(-nyquist,nyquist,delta)

    # turn this into a positive frequency array

    pos_freq = freq[(len(ffty)/2):]

    # make the transfer function for the first half of the data

    i_f_3db = min( nm.where(pos_freq >= f_3db)[0] )

    f_min = f_3db - (width/2.0)

    i_f_min = min( nm.where(pos_freq >= f_min)[0] )

    f_max = f_3db + (width/2);

    i_f_max = min( nm.where(pos_freq >= f_max)[0] )

    transfer_function = nm.zeros(len(y)/2)

    transfer_function[0:i_f_min] = 1

    transfer_function[i_f_min:i_f_max] = (1 + nm.sin(-nm.pi * ((freq[i_f_min:i_f_max] - freq[i_f_3db])/width)))/2.0

    transfer_function[i_f_max:(len(freq)/2)] = 0

    # symmetrize this to be [0 0 0 ... .8 .9 1 1 1 1 1 1 1 1 .9 .8 ... 0 0 0] to match the FFT

    transfer_function = nm.append(nm.flipud(transfer_function),transfer_function)

    # plot up the transfer function

    # since "freq" is only the positive frequencies, select out

    pylab.figure(1)

    pylab.clf()

    pylab.plot(freq,transfer_function)

    pylab.xlabel('Frequency [Hz]')

    pylab.ylabel('Filter Transfer Function')

    pylab.xlim([-10.0,10.0])

    pylab.ylim([-0.05,1.05])

    # apply the filter, undo the fft shift, and invert the fft

    filtered=nm.real(scipy.fftpack.ifft(scipy.fftpack.ifftshift(ffty*transfer_function)))

    # remove the padd, if we applied it

    if padd_data:

        filtered = filtered[(len(y)/3):(2*(len(y)/3))]

    # return the filtered data

    return filtered

# do an example of lowpass filtering

# first make some fake data

# a sine wave fluctuating once every pi seconds

# samples 1000 times per second

fakedata = nm.sin(nm.arange(0,11,0.001)) + nm.random.randn(len(nm.arange(0,11,0.001)))/4.0

# run the filter

# lowpass at 5 Hz, with a 1 Hz width of its roll-off

filtered = lowpass_cosine(fakedata,0.001,5.0,1.0,padd_data=True)

# plot the noisy data, with the filtered data on top

pylab.figure(2)

pylab.clf()

pylab.plot(nm.arange(0,11,0.001),fakedata,label='Noisy Data')

pylab.plot(nm.arange(0,11,0.001),filtered,label='Lowpass Filtered Data')

pylab.xlabel('Time [s]')

pylab.ylabel('Voltage')

pylab.legend()

pylab.ion()

pylab.show()

#codingexercise: CodingExercise-07-29-2026.docx

 Agentic judges for drone image analytics

Andrew Ng’s agentic workflow pattern—reflection, tool use, planning, and multi-agent collaboration—applies to drone-vision benchmarking. Reflection lets an agent critique and revise detections, captions, SQL answers, or mission reports. Tool use grounds reasoning in retrieval, code execution, geospatial operators, detector APIs, and database queries. Planning decomposes a workload into explicit steps before execution and supports replanning when a probe fails. Multi-agent orchestration assigns specialized roles: one agent checks geospatial consistency, another evaluates temporal coherence, another audits semantic alignment, and an arbiter aggregates evidence into a score.

Memory has been a design constraint. Loops let agents think; graphs let agents remember. A reflection loop can improve a single answer, but without persistent state the agent forgets why it inspected a frame, which detector disagreed, or which workload constraint failed. A graph turns those transient observations into reusable memory: frames, objects, captions, detections, SQL results, tool calls, critiques, and final judgments become linked evidence rather than buried transcript text. In ezbenchmark, this converts an agentic judge from a one-pass evaluator into a stateful audit system.

A practical build path is incremental. First, add one critique call after every generated answer or score; this is the highest-return change because it catches unsupported claims, missing visual evidence, and weak workload alignment before output is finalized. Second, expose the judge to tools: vector retrieval over frame evidence, SQL over the scene catalog, detector re-runs, spatial predicates, and temporal-neighbor comparisons. Third, require the judge to emit a structured plan before execution, such as JSON steps with expected evidence, tool calls, and failure conditions. Fourth, split evaluation across specialized agents and connect them through a shared graph store. The result is not merely a stronger prompt; it is an architecture in which weak models can outperform stronger single-pass models because the workflow supplies iteration, grounding, decomposition, and durable memory.

Agents are usually dedicated to perception, reasoning, and control in different ways. Sapkota et al. introduce the term “Agentic UAVs” to describe systems that integrate perception, cognition, control, and communication into layered, goal-driven agents that operate with contextual reasoning and memory, rather than fixed scripts or reactive control loops [1]. In their framework, aerial image understanding is only one layer in a broader cognitive stack: perception agents extract structure from imagery and other sensors; cognitive agents plan and replan missions; control agents execute trajectories; and communication agents coordinate with humans and other UAVs. This layered view is useful when we start thinking about agentic frameworks as “judges” for benchmarking: the judging capability can itself be an agent, sitting in the cognition layer, consuming outputs from perception agents and workload metadata rather than raw pixels alone [1].

Vision–language–driven agents are a distinct subclass. Sapkota et al. explicitly highlight vision–language models and multimodal sensing as key enabling technologies for Agentic UAVs, noting that they allow agents to parse complex scenes, follow natural-language instructions, and ground symbolic goals in visual context [1]. These agents differ from traditional planners in that they can reason over image and text jointly, which makes them natural candidates for roles like “mission explainer,” “anomaly triager,” or, in our case, “benchmark judge” for aerial analytics workloads. Instead of judging purely from numeric metrics, a vision–language agent can look at a drone scene, read a workload description, inspect candidate outputs, and form a qualitative judgment about which pipeline better captures the intended analytic semantics [1].

#codingexercise: CodingExercise-07-29-2026.docx

 

A minimal reproduction: tracking a vehicle by its Fourier signature

Let’s track a red vehicle across a short sequence of aerial drone frames purely from its Fourier-descriptor "shape signature," adapted from the public timfeirg/Fourier-Descriptors approach. The pipeline has three stages, mirroring the detect → describe → match structure of the Seidaliyeva et al. system:

import cv2
import numpy as np

def extract_red_car_contour(image):
    # Convert to HSV and threshold to isolate red color (works for most red cars)
    hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)
    # Lower and upper bounds for red in HSV (red wraps around hue 0/180)
    lower_red1 = np.array([0, 70, 50])
    upper_red1 = np.array([10, 255, 255])
    lower_red2 = np.array([170, 70, 50])
    upper_red2 = np.array([180, 255, 255])
    mask1 = cv2.inRange(hsv, lower_red1, upper_red1)
    mask2 = cv2.inRange(hsv, lower_red2, upper_red2)
    mask = mask1 | mask2

    # Morphological operations to clean the mask
    kernel = np.ones((5, 5), np.uint8)
    mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel)
    mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, kernel)

    # Find contours; assume the single largest red blob is the target vehicle
    contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    if contours:
        return max(contours, key=cv2.contourArea), mask
    return None, mask

def compute_fourier_descriptors(contour, num_descriptors=16):
    contour = contour.squeeze()
    # Treat each contour point (x, y) as a complex number z = x + iy
    z = contour[:, 0] + 1j * contour[:, 1]
    # DFT of the boundary sequence
    fd = np.fft.fft(z)
    # Normalize (naive version — see critique below)
    fd = fd[:num_descriptors] / np.abs(fd)
    return fd

def match_descriptors(fd1, fd2):
    return np.linalg.norm(fd1 - fd2)

Each frame's red-vehicle silhouette is boiled down to a 1D closed boundary curve, which is then treated as a discrete complex-valued signal sampled at contour points as the boundary is traced. Applying the 1D DFT, , yields a set of complex Fourier descriptors that jointly reconstruct the exact boundary shape. Truncating to the first num_descriptors low-frequency coefficients (here, 16) retains the coarse silhouette — overall size, elongation, major concavities — while discarding high-frequency terms that mostly encode pixel-level contour jitter and are therefore both noise-dominated and expensive to keep. Two frames' descriptor vectors are then compared with a simple Euclidean distance: a small distance across consecutive frames signals "same object, successfully tracked," while a large jump flags a possible identity switch or lost target — exactly the frame-to-frame association step needed to build a tracklet from raw per-frame detections.

A normalization subtlety worth getting right

The compute_fourier_descriptors function above uses a compact one-line normalization, fd[:num_descriptors] / np.abs(fd), dividing each retained coefficient by its own magnitude. This is a reasonable first pass, and it is invariant to starting point and to overall rotation-and-scale jointly, but it does not deliver the individually decoupled translation, scale, and rotation invariances that Fourier descriptors are prized for in the classical shape-analysis literature, and it silently discards information a more careful normalization would keep. The canonical procedure, going back to Zahn and Roskies' foundational 1972 paper on Fourier descriptors for plane closed curves (Zahn & Roskies, IEEE Trans. Computers, 1972) and formalized for practical implementation by Folkers and Samet (Univ. of Maryland / ICPR 2002), separates the three invariances explicitly:

·        Translation invariance: the DC term is the centroid of the contour (the mean of all boundary points); it carries no shape information at all, only the object's position in the frame. The correct procedure discards or zeroes before comparison rather than folding it into a shared denominator.

·        Scale invariance: dividing every retained coefficient by (the magnitude of the first non-trivial harmonic) rather than by each coefficient's own magnitude normalizes for object size while preserving the relative magnitude relationships between harmonics — which is where the shape information actually lives. Dividing each by its own , as in the code above, forces every normalized coefficient to have magnitude exactly 1, which erases the relative-magnitude structure that distinguishes, say, an elongated silhouette from a compact one.

·        Rotation and starting-point invariance: multiplying every coefficient by a phase correction derived from the phase of (and, for full starting-point invariance, an additional per-harmonic phase term) removes sensitivity to the object's orientation and to which boundary pixel the contour-tracing algorithm happened to start from, without discarding magnitude information the way the naive scheme implicitly risks when phase and magnitude are conflated.

A related normalized-contour-signature approach, useful when descriptors are pooled across a large object gallery, is described in Sokić and colleagues' NCC-signature work (Sokić et al., CBMI 2014). The practical takeaway for anyone adapting the toy pipeline above into a production drone-tracking module is small in code footprint but large in effect: replace the single-line fd / np.abs(fd) normalization with followed by division by and a phase rotation referenced to , which is a handful of extra lines but converts a descriptor that is only jointly invariant under specific conditions into one with the individually decoupled translation-, scale-, and rotation-invariance guarantees the classical Fourier-descriptor literature actually promises — the same guarantees that let the Seidaliyeva et al. system recognize a drone regardless of its heading or distance from the camera.

Suggested fix:

import numpy as np

def compute_fourier_descriptors(contour, num_descriptors=16):

    """

    Compute translation-, scale-, rotation-, and starting-point-invariant

    Fourier descriptors for a closed 2D contour.

 

    Fixes the naive normalization `fd[:num_descriptors] / np.abs(fd)`, which

    only cancels a *combined* rotation+scale ambiguity and leaves the

    translation-encoding DC term and the starting-point ambiguity untouched.

    """

    contour = np.asarray(contour, dtype=np.float64).squeeze()

    if contour.ndim != 2 or contour.shape[1] != 2:

        raise ValueError(f"Expected an (N, 2) contour, got shape {contour.shape}")

 

    N = contour.shape[0]

    if N < 4:

        raise ValueError(f"Contour has only {N} points; need at least a few to be meaningful")

 

    # Treat each contour point (x, y) as a complex number z = x + iy

    z = contour[:, 0] + 1j * contour[:, 1]

 

    # DFT of the boundary sequence: a(k) for k = 0 .. N-1

    a = np.fft.fft(z)

 

    # --- Translation invariance ---

    # a(0) = sum(z) is N times the contour's centroid -- it encodes only

    # position, never shape, so it's discarded rather than folded into a

    # shared denominator (Zahn & Roskies, 1972). a[0] is never referenced below.

 

    # --- Scale invariance ---

    # Divide by |a(1)| (the fundamental harmonic's magnitude), not by each

    # coefficient's own magnitude -- this preserves the *relative* magnitude

    # relationships between harmonics, which is where the shape lives

    # (Folkers & Samet, ICPR 2002).

    a1_mag = np.abs(a[1])

    if a1_mag < 1e-9:

        raise ValueError(

            "First harmonic |a(1)| is ~0 (degenerate or point-symmetric contour); "

            "cannot scale-normalize."

        )

 

    # --- Rotation invariance ---

    # A global rotation by theta multiplies EVERY coefficient by the same

    # constant e^{i*theta} (DFT linearity); that constant phase cancels once

    # we take magnitudes.

    #

    # --- Starting-point invariance ---

    # A different starting point shifts a(k) by a k-dependent, unit-magnitude

    # phase e^{-i*2*pi*k*n0/N} (DFT shift theorem); it too vanishes under

    # magnitude. Together, |a(k)| / |a(1)| is simultaneously rotation- and

    # starting-point-invariant, on top of translation/scale invariance.

    num_harmonics = min(num_descriptors, N - 2)

    descriptors = np.abs(a[2:2 + num_harmonics]) / a1_mag

 

    if num_harmonics < num_descriptors:

        descriptors = np.pad(descriptors, (0, num_descriptors - num_harmonics))

 

    return descriptors

 

 

def match_descriptors(fd1, fd2):

    """Euclidean distance between two (now real-valued) descriptor vectors."""

    return np.linalg.norm(np.asarray(fd1) - np.asarray(fd2))

 

Correction Test Results:

Transform applied to the same shape	Naive distance	Fixed distance
Translated (+137, -52)	3.07	0
Scaled x3.7	0.97	0
Rotated 63°	4.08	0
Scale x0.4 + rotate 205° + translate	6.45	0
Starting point cyclically shifted by 37 samples	5.49	0
All four combined at once	5.82	0
A genuinely different shape (sanity check)	N/A	0.1

 

Sample invocation:

Import cv2
import numpy as np

image_paths = ['frame11.jpg'] # Replace with actual filepaths

red_car_descriptors = []

for image_path in image_paths:

    image = cv2.imread(image_path)

    contour = extract_red_car_contour(image)

    if contour is not None and len(contour) > 10:

        fd = compute_fourier_descriptors(contour)

        red_car_descriptors.append(fd)

        # Compute bounding box

        x, y, w, h = cv2.boundingRect(contour)

        cv2.rectangle(image, (x, y), (x + w, y + h), (0, 255, 0), 2)

        cv2.imwrite(image_path.replace(".jpg","")+"-tracked.jpg", image)

        cv2.imshow('Red Car Tracking', image)

        cv2.waitKey(3000)

 

cv2.waitKey(10000)

cv2.destroyAllWindows()

 

#Codingexercise: Codingexercise-07-28-2026.docx