R-tree 与空间索引：PostGIS 的底层结构

#define MAX_ENTRIES 50   /* M：每个节点最多条目数 */
#define MIN_ENTRIES 20   /* m：每个节点最少条目数 */

typedef struct rtree_node {
    int count;                       /* 当前条目数 */
    int is_leaf;                     /* 是否为叶节点 */
    Rect rects[MAX_ENTRIES];         /* MBR 数组 */
    struct rtree_node *children[MAX_ENTRIES]; /* 子节点指针（叶节点为 NULL） */
    int data_ids[MAX_ENTRIES];       /* 数据 ID（仅叶节点使用） */
} RTreeNode;

/* 计算二维坐标的 Hilbert 值 */
unsigned int xy_to_hilbert(unsigned int x, unsigned int y, int order)
{
    unsigned int d = 0;
    for (int s = order / 2; s > 0; s /= 2) {
        unsigned int rx = (x & s) > 0 ? 1 : 0;
        unsigned int ry = (y & s) > 0 ? 1 : 0;
        d += s * s * ((3 * rx) ^ ry);
        /* 旋转象限 */
        if (ry == 0) {
            if (rx == 1) {
                x = s - 1 - x;
                y = s - 1 - y;
            }
            unsigned int t = x;
            x = y;
            y = t;
        }
    }
    return d;
}

-- 创建空间表
CREATE TABLE restaurants (
    id SERIAL PRIMARY KEY,
    name TEXT,
    geom GEOMETRY(Point, 4326)  -- WGS84 坐标系
);

-- 创建 GiST 空间索引
CREATE INDEX idx_restaurants_geom ON restaurants USING gist(geom);

-- 范围查询：找出边界框内的餐厅
SELECT name FROM restaurants
WHERE geom && ST_MakeEnvelope(116.3, 39.9, 116.5, 40.0, 4326);

-- 距离查询：找出 1km 内的餐厅
SELECT name, ST_Distance(
    geom::geography,
    ST_SetSRID(ST_MakePoint(116.4, 39.95), 4326)::geography
) AS dist
FROM restaurants
WHERE ST_DWithin(
    geom::geography,
    ST_SetSRID(ST_MakePoint(116.4, 39.95), 4326)::geography,
    1000  -- 1000 米
)
ORDER BY dist;

-- k-NN 查询：找最近的 5 家餐厅（利用 GiST 的 KNN 支持）
SELECT name, geom <-> ST_SetSRID(ST_MakePoint(116.4, 39.95), 4326) AS dist
FROM restaurants
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(116.4, 39.95), 4326)
LIMIT 5;

-- 对按空间顺序插入的数据，BRIN 索引极小但有效
CREATE INDEX idx_tracks_brin ON gps_tracks USING brin(geom);

typedef struct quad_node {
    double cx, cy;           /* 中心点 */
    double half_w, half_h;   /* 半宽和半高 */
    int point_count;
    double points_x[4];     /* 最多 4 个点 */
    double points_y[4];
    struct quad_node *nw, *ne, *sw, *se;
} QuadNode;

/* 计算 Morton code（Z-order value） */
unsigned int morton_encode(unsigned int x, unsigned int y)
{
    x = (x | (x << 8)) & 0x00FF00FF;
    x = (x | (x << 4)) & 0x0F0F0F0F;
    x = (x | (x << 2)) & 0x33333333;
    x = (x | (x << 1)) & 0x55555555;

    y = (y | (y << 8)) & 0x00FF00FF;
    y = (y | (y << 4)) & 0x0F0F0F0F;
    y = (y | (y << 2)) & 0x33333333;
    y = (y | (y << 1)) & 0x55555555;

    return x | (y << 1);
}

/*
 * rtree.c — 完整的二维 R-tree 实现
 *
 * 支持：插入、范围查询、k-NN 查询
 * 分裂策略：二次分裂（Quadratic Split）
 *
 * 编译：gcc -O2 -o rtree rtree.c -lm
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>

/* ===== 配置参数 ===== */
#define DIMS        2
#define MAX_ENTRIES 8     /* 为便于演示，使用较小的 M */
#define MIN_ENTRIES 3     /* m = ceil(M * 0.4) */
#define KNN_HEAP_CAP 256

/* ===== MBR（最小边界矩形） ===== */
typedef struct {
    double lo[DIMS];
    double hi[DIMS];
} Rect;

static Rect rect_empty(void)
{
    Rect r;
    for (int d = 0; d < DIMS; d++) {
        r.lo[d] = DBL_MAX;
        r.hi[d] = -DBL_MAX;
    }
    return r;
}

static Rect rect_from_point(double x, double y)
{
    Rect r;
    r.lo[0] = r.hi[0] = x;
    r.lo[1] = r.hi[1] = y;
    return r;
}

static int rect_overlaps(const Rect *a, const Rect *b)
{
    for (int d = 0; d < DIMS; d++) {
        if (a->lo[d] > b->hi[d] || a->hi[d] < b->lo[d])
            return 0;
    }
    return 1;
}

static int rect_contains(const Rect *a, const Rect *b)
{
    for (int d = 0; d < DIMS; d++) {
        if (a->lo[d] > b->lo[d] || a->hi[d] < b->hi[d])
            return 0;
    }
    return 1;
}

static double rect_area(const Rect *r)
{
    double area = 1.0;
    for (int d = 0; d < DIMS; d++) {
        double side = r->hi[d] - r->lo[d];
        if (side <= 0.0) return 0.0;
        area *= side;
    }
    return area;
}

static Rect rect_expand(const Rect *a, const Rect *b)
{
    return rect_union(a, b);
}

static double rect_enlargement(const Rect *a, const Rect *b)
{
    Rect u = rect_union(a, b);
    return rect_area(&u) - rect_area(a);
}

/* 点到 MBR 的最小距离（用于 k-NN） */
static double rect_min_dist(const Rect *r, double px, double py)
{
    double dx = 0.0, dy = 0.0;

    if (px < r->lo[0])      dx = r->lo[0] - px;
    else if (px > r->hi[0]) dx = px - r->hi[0];

    if (py < r->lo[1])      dy = r->lo[1] - py;
    else if (py > r->hi[1]) dy = py - r->hi[1];

    return sqrt(dx * dx + dy * dy);
}

/* ===== R-tree 节点 ===== */
typedef struct rtree_node {
    int count;
    int is_leaf;
    int height;      /* 叶节点 height=0 */
    Rect rects[MAX_ENTRIES + 1];   /* +1 用于溢出时暂存 */
    union {
        struct rtree_node *children[MAX_ENTRIES + 1];
        int data_ids[MAX_ENTRIES + 1];
    };
} Node;

typedef struct {
    Node *root;
    int size;        /* 数据对象总数 */
    int height;      /* 树高 */
} RTree;

static Node *node_new(int is_leaf, int height)
{
    Node *n = (Node *)calloc(1, sizeof(Node));
    n->is_leaf = is_leaf;
    n->height = height;
    return n;
}

static Rect node_cover(const Node *n)
{
    Rect cover = n->rects[0];
    for (int i = 1; i < n->count; i++)
        cover = rect_union(&cover, &n->rects[i]);
    return cover;
}

/* ===== 二次分裂（Quadratic Split） ===== */
static void pick_seeds(Node *n, int *s1, int *s2)
{
    double worst = -DBL_MAX;
    *s1 = 0;
    *s2 = 1;

    for (int i = 0; i < n->count; i++) {
        for (int j = i + 1; j < n->count; j++) {
            Rect combined = rect_union(&n->rects[i], &n->rects[j]);
            double waste = rect_area(&combined)
                         - rect_area(&n->rects[i])
                         - rect_area(&n->rects[j]);
            if (waste > worst) {
                worst = waste;
                *s1 = i;
                *s2 = j;
            }
        }
    }
}

static void quadratic_split(Node *n, Node **out_a, Node **out_b)
{
    int total = n->count;  /* 此时 count == MAX_ENTRIES + 1 */
    int assigned[MAX_ENTRIES + 1];
    memset(assigned, 0, sizeof(assigned));

    int s1, s2;
    pick_seeds(n, &s1, &s2);

    Node *a = node_new(n->is_leaf, n->height);
    Node *b = node_new(n->is_leaf, n->height);

    /* 将种子分配到两组 */
    a->rects[0] = n->rects[s1];
    if (n->is_leaf) a->data_ids[0] = n->data_ids[s1];
    else            a->children[0] = n->children[s1];
    a->count = 1;
    assigned[s1] = 1;

    b->rects[0] = n->rects[s2];
    if (n->is_leaf) b->data_ids[0] = n->data_ids[s2];
    else            b->children[0] = n->children[s2];
    b->count = 1;
    assigned[s2] = 1;

    int remaining = total - 2;

    while (remaining > 0) {
        /* 检查是否需要将剩余全部分配给某组 */
        if (a->count + remaining == MIN_ENTRIES) {
            for (int i = 0; i < total; i++) {
                if (assigned[i]) continue;
                a->rects[a->count] = n->rects[i];
                if (n->is_leaf) a->data_ids[a->count] = n->data_ids[i];
                else            a->children[a->count] = n->children[i];
                a->count++;
                assigned[i] = 1;
            }
            break;
        }
        if (b->count + remaining == MIN_ENTRIES) {
            for (int i = 0; i < total; i++) {
                if (assigned[i]) continue;
                b->rects[b->count] = n->rects[i];
                if (n->is_leaf) b->data_ids[b->count] = n->data_ids[i];
                else            b->children[b->count] = n->children[i];
                b->count++;
                assigned[i] = 1;
            }
            break;
        }

        /* 选择下一个条目：最大偏好差异 */
        Rect cover_a = node_cover(a);
        Rect cover_b = node_cover(b);

        int best_idx = -1;
        double best_diff = -DBL_MAX;
        int best_group = 0;

        for (int i = 0; i < total; i++) {
            if (assigned[i]) continue;
            double ea = rect_enlargement(&cover_a, &n->rects[i]);
            double eb = rect_enlargement(&cover_b, &n->rects[i]);
            double diff = fabs(ea - eb);
            if (diff > best_diff) {
                best_diff = diff;
                best_idx = i;
                best_group = (ea < eb) ? 0 : 1;
            }
        }

        Node *target = (best_group == 0) ? a : b;
        target->rects[target->count] = n->rects[best_idx];
        if (n->is_leaf) target->data_ids[target->count] = n->data_ids[best_idx];
        else            target->children[target->count] = n->children[best_idx];
        target->count++;
        assigned[best_idx] = 1;
        remaining--;
    }

    *out_a = a;
    *out_b = b;
}

/* ===== 插入 ===== */
static int choose_subtree(Node *n, const Rect *r)
{
    int best = 0;
    double best_enlarge = rect_enlargement(&n->rects[0], r);
    double best_area = rect_area(&n->rects[0]);

    for (int i = 1; i < n->count; i++) {
        double enlarge = rect_enlargement(&n->rects[i], r);
        double area = rect_area(&n->rects[i]);
        if (enlarge < best_enlarge ||
            (enlarge == best_enlarge && area < best_area)) {
            best = i;
            best_enlarge = enlarge;
            best_area = area;
        }
    }
    return best;
}

/*
 * 递归插入，返回分裂产生的新节点（如果没有分裂则返回 NULL）
 * split_rect 返回新节点的 MBR
 */
static Node *insert_internal(Node *node, const Rect *r, int data_id,
                             int target_height, Rect *split_rect)
{
    if (node->height == target_height) {
        /* 在当前节点插入条目 */
        node->rects[node->count] = *r;
        if (node->is_leaf)
            node->data_ids[node->count] = data_id;
        else
            node->children[node->count] = NULL; /* 由调用者设置 */
        node->count++;

        if (node->count <= MAX_ENTRIES)
            return NULL;

        /* 需要分裂 */
        Node *a, *b;
        quadratic_split(node, &a, &b);

        /* 复制 a 的内容到当前节点 */
        node->count = a->count;
        memcpy(node->rects, a->rects, sizeof(Rect) * a->count);
        if (node->is_leaf)
            memcpy(node->data_ids, a->data_ids, sizeof(int) * a->count);
        else
            memcpy(node->children, a->children, sizeof(Node *) * a->count);
        free(a);

        *split_rect = node_cover(b);
        return b;
    }

    /* 选择子树 */
    int idx = choose_subtree(node, r);
    Rect child_split_rect;
    Node *new_child = insert_internal(node->children[idx], r, data_id,
                                       target_height, &child_split_rect);

    /* 更新该子树的 MBR */
    node->rects[idx] = node_cover(node->children[idx]);

    if (new_child == NULL)
        return NULL;

    /* 子节点分裂了，需要在当前节点加入新子节点 */
    node->rects[node->count] = child_split_rect;
    node->children[node->count] = new_child;
    node->count++;

    if (node->count <= MAX_ENTRIES)
        return NULL;

    /* 当前节点也需要分裂 */
    Node *a, *b;
    quadratic_split(node, &a, &b);

    node->count = a->count;
    memcpy(node->rects, a->rects, sizeof(Rect) * a->count);
    memcpy(node->children, a->children, sizeof(Node *) * a->count);
    free(a);

    *split_rect = node_cover(b);
    return b;
}

void rtree_insert(RTree *tree, double x, double y, int data_id)
{
    Rect r = rect_from_point(x, y);

    if (tree->root == NULL) {
        tree->root = node_new(1, 0);
        tree->root->rects[0] = r;
        tree->root->data_ids[0] = data_id;
        tree->root->count = 1;
        tree->size = 1;
        tree->height = 1;
        return;
    }

    Rect split_rect;
    Node *split_node = insert_internal(tree->root, &r, data_id,
                                        0, &split_rect);
    tree->size++;

    if (split_node != NULL) {
        /* 根节点分裂，创建新根 */
        Node *new_root = node_new(0, tree->root->height + 1);
        new_root->rects[0] = node_cover(tree->root);
        new_root->children[0] = tree->root;
        new_root->rects[1] = split_rect;
        new_root->children[1] = split_node;
        new_root->count = 2;
        tree->root = new_root;
        tree->height++;
    }
}

/* ===== 范围查询 ===== */
int rtree_search(Node *node, const Rect *query, int *results, int max_results)
{
    if (node == NULL) return 0;

    int found = 0;
    for (int i = 0; i < node->count && found < max_results; i++) {
        if (!rect_overlaps(&node->rects[i], query))
            continue;

        if (node->is_leaf) {
            results[found++] = node->data_ids[i];
        } else {
            found += rtree_search(node->children[i], query,
                                  results + found, max_results - found);
        }
    }
    return found;
}

/* ===== k-NN 查询（优先级队列方法） ===== */
typedef struct {
    double dist;
    int is_leaf_entry;
    union {
        Node *node;
        int data_id;
    };
    Rect rect;
} KNNEntry;

typedef struct {
    KNNEntry *items;
    int size;
    int cap;
} MinHeap;

static MinHeap *heap_new(int cap)
{
    MinHeap *h = (MinHeap *)malloc(sizeof(MinHeap));
    h->items = (KNNEntry *)malloc(sizeof(KNNEntry) * cap);
    h->size = 0;
    h->cap = cap;
    return h;
}

static void heap_push(MinHeap *h, KNNEntry e)
{
    if (h->size >= h->cap) {
        h->cap *= 2;
        h->items = (KNNEntry *)realloc(h->items, sizeof(KNNEntry) * h->cap);
    }
    int i = h->size++;
    h->items[i] = e;
    while (i > 0) {
        int parent = (i - 1) / 2;
        if (h->items[parent].dist <= h->items[i].dist) break;
        KNNEntry tmp = h->items[parent];
        h->items[parent] = h->items[i];
        h->items[i] = tmp;
        i = parent;
    }
}

static KNNEntry heap_pop(MinHeap *h)
{
    KNNEntry top = h->items[0];
    h->items[0] = h->items[--h->size];
    int i = 0;
    while (1) {
        int left = 2 * i + 1;
        int right = 2 * i + 2;
        int smallest = i;
        if (left < h->size && h->items[left].dist < h->items[smallest].dist)
            smallest = left;
        if (right < h->size && h->items[right].dist < h->items[smallest].dist)
            smallest = right;
        if (smallest == i) break;
        KNNEntry tmp = h->items[smallest];
        h->items[smallest] = h->items[i];
        h->items[i] = tmp;
        i = smallest;
    }
    return top;
}

/*
 * k-NN 搜索：使用"最佳优先搜索"策略
 * 将节点和数据点统一放入优先级队列，按到查询点的最小距离排序
 */
int rtree_knn(RTree *tree, double qx, double qy, int k,
              int *result_ids, double *result_dists)
{
    if (tree->root == NULL || k <= 0) return 0;

    MinHeap *pq = heap_new(KNN_HEAP_CAP);
    int found = 0;

    /* 将根节点的所有条目加入优先队列 */
    for (int i = 0; i < tree->root->count; i++) {
        KNNEntry e;
        e.dist = rect_min_dist(&tree->root->rects[i], qx, qy);
        e.rect = tree->root->rects[i];
        if (tree->root->is_leaf) {
            e.is_leaf_entry = 1;
            e.data_id = tree->root->data_ids[i];
        } else {
            e.is_leaf_entry = 0;
            e.node = tree->root->children[i];
        }
        heap_push(pq, e);
    }

    while (pq->size > 0 && found < k) {
        KNNEntry top = heap_pop(pq);

        if (top.is_leaf_entry) {
            result_ids[found] = top.data_id;
            result_dists[found] = top.dist;
            found++;
        } else {
            /* 展开内部节点或叶节点 */
            Node *child = top.node;
            for (int i = 0; i < child->count; i++) {
                KNNEntry e;
                e.dist = rect_min_dist(&child->rects[i], qx, qy);
                e.rect = child->rects[i];
                if (child->is_leaf) {
                    e.is_leaf_entry = 1;
                    e.data_id = child->data_ids[i];
                } else {
                    e.is_leaf_entry = 0;
                    e.node = child->children[i];
                }
                heap_push(pq, e);
            }
        }
    }

    free(pq->items);
    free(pq);
    return found;
}

/* ===== 释放内存 ===== */
void rtree_free(Node *node)
{
    if (node == NULL) return;
    if (!node->is_leaf) {
        for (int i = 0; i < node->count; i++)
            rtree_free(node->children[i]);
    }
    free(node);
}

/* ===== 测试代码 ===== */
int main(void)
{
    RTree tree = {NULL, 0, 0};

    /* 插入一些二维点 */
    struct { double x, y; } points[] = {
        {1.0, 2.0}, {3.0, 4.0}, {5.0, 1.0}, {7.0, 3.0},
        {2.0, 6.0}, {8.0, 8.0}, {4.0, 5.0}, {6.0, 7.0},
        {9.0, 2.0}, {1.5, 8.5}, {3.5, 3.5}, {5.5, 5.5},
        {7.5, 1.5}, {2.5, 4.5}, {4.5, 7.5}, {6.5, 6.5},
        {8.5, 4.5}, {0.5, 0.5}, {9.5, 9.5}, {5.0, 5.0},
    };
    int n = sizeof(points) / sizeof(points[0]);

    printf("=== 插入 %d 个点 ===\n", n);
    for (int i = 0; i < n; i++)
        rtree_insert(&tree, points[i].x, points[i].y, i);

    printf("树高: %d, 数据量: %d\n\n", tree.height, tree.size);

    /* 范围查询 */
    Rect query;
    query.lo[0] = 2.0; query.lo[1] = 2.0;
    query.hi[0] = 6.0; query.hi[1] = 6.0;

    int results[100];
    int count = rtree_search(tree.root, &query, results, 100);

    printf("=== 范围查询 [2,2]-[6,6] ===\n");
    printf("找到 %d 个点:\n", count);
    for (int i = 0; i < count; i++) {
        int id = results[i];
        printf("  点 %d: (%.1f, %.1f)\n", id, points[id].x, points[id].y);
    }

    /* k-NN 查询 */
    printf("\n=== k-NN 查询: 最近 5 个点（查询点 5.0, 5.0） ===\n");
    int knn_ids[5];
    double knn_dists[5];
    int knn_count = rtree_knn(&tree, 5.0, 5.0, 5, knn_ids, knn_dists);

    for (int i = 0; i < knn_count; i++) {
        int id = knn_ids[i];
        printf("  #%d: 点 %d (%.1f, %.1f), 距离 = %.3f\n",
               i + 1, id, points[id].x, points[id].y, knn_dists[i]);
    }

    rtree_free(tree.root);
    printf("\n测试完成。\n");
    return 0;
}

gcc -O2 -o rtree rtree.c -lm
./rtree

-- 多维属性查询：找出价格在 200-500 且面积在 80-120 的房源
-- 如果在 (price, area) 上建立多维索引
SELECT * FROM houses
WHERE price BETWEEN 200 AND 500
  AND area BETWEEN 80 AND 120;