最大公共子图算法

下面给出一个示例 C++ 程序，用“图论+最大团(或最大公共子图)”的思想来在两组 3D 点中寻找最大对应集合。原理简述如下：

1. 读入两组 3D 点
   分别从两个 txt 文件中读入三维坐标，保存到 vector<Point3D> 容器中。

2. 构造“兼容图”(Compatibility Graph)

   • 我们有两组点：
     • 第一组 $\{S_i\}$ (大小为 $N$)
     • 第二组 $\{T_j\}$ (大小为 $M$)
   • 在“兼容图”中，每个“结点”对应一个可能的配对 $(S_i, T_j)$。因此兼容图总共有 $N\times M$ 个结点。
   • 对兼容图中任意两个结点 $(S_{i1}, T_{j1})$ 和 $(S_{i2}, T_{j2})$，如果满足 $$\bigl|\|S_{i1} - S_{i2}\| - \|T_{j1} - T_{j2}\|\bigr| < \text{阈值}$$ （例如阈值为 1.0f），则在这两个结点之间连一条边，表示它们“兼容”。

3. 寻找最大团 (Maximum Clique)

   • 由于我们要找的是最大公共子图，也就是最多的“一一配对”，并且配对之间的距离关系要一致(在容忍度内)，这可以转换成在上述“兼容图”中找最大团的问题：
   • 一个团(clique)就是一个完全子图：其中任何两个结点都两两相连。换言之，一个团对应了一组两两“距离匹配”的点对。
   • 找到“兼容图”中的最大团后，它所包含的结点个数，就是最大公共子图里匹配到的点对数；对应的结点则是具体的匹配关系。

4. 输出结果

   • 最后把最大团里的那些配对 $(S_i, T_j)$ 输出，即可得到对应点对列表，并可计算每对之间的真实距离误差。

下面的示例代码演示了这一流程，供参考（为方便演示，代码写得比较直白，未做大规模优化或工程封装；若点数较多，可考虑更高效的数据结构或剪枝策略）。

cpp
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <string>
#include <algorithm>

//--------------------------------------
// 数据结构定义
//--------------------------------------
struct Point3D {
    double x, y, z;
};

static double distance3D(const Point3D& p1, const Point3D& p2) {
    double dx = p1.x - p2.x;
    double dy = p1.y - p2.y;
    double dz = p1.z - p2.z;
    return std::sqrt(dx*dx + dy*dy + dz*dz);
}

//--------------------------------------
// 从TXT文件读取点集，每行格式： x y z
//--------------------------------------
std::vector<Point3D> readPointsFromFile(const std::string& filename)
{
    std::vector<Point3D> points;
    std::ifstream fin(filename);
    if(!fin.is_open()){
        std::cerr << "Failed to open file: " << filename << std::endl;
        return points;
    }
    while(!fin.eof()){
        Point3D p;
        fin >> p.x >> p.y >> p.z;
        if(!fin.fail()){
            points.push_back(p);
        }
    }
    fin.close();
    return points;
}

//--------------------------------------
// 兼容图：节点用 (i, j) 表示第一个点集里的 i 与第二个点集里的 j
// adjacency[u] 存放与 u 相连(兼容)的所有节点 v
//--------------------------------------
struct PairNode {
    int i; // 对应第一组点的下标
    int j; // 对应第二组点的下标
};

//--------------------------------------
// 在所有 (i, j) 与 (k, l) 之间构图
// 如果距离匹配(误差 < tolerance)，则连一条边
//--------------------------------------
std::vector<std::vector<int>> buildCompatibilityGraph(
    const std::vector<Point3D>& srcPoints,
    const std::vector<Point3D>& tgtPoints,
    std::vector<PairNode>& pairNodes, 
    double tolerance)
{
    // 先构造所有可能的配对节点
    // pairNodes[u] = (i, j) 表示 srcPoints[i] <-> tgtPoints[j]
    int N = (int)srcPoints.size();
    int M = (int)tgtPoints.size();

    pairNodes.clear();
    pairNodes.reserve(N*M);
    for(int i = 0; i < N; ++i) {
        for(int j = 0; j < M; ++j) {
            PairNode pn;
            pn.i = i;
            pn.j = j;
            pairNodes.push_back(pn);
        }
    }

    // 构造邻接表 adjacency，大小为 N*M
    std::vector<std::vector<int>> adjacency(N*M);

    // 两层循环遍历所有节点对 (u, v)，检查距离是否匹配
    // 为了减少重复计算，可以先把 srcPoints、tgtPoints 的距离预先缓存下来
    // 这里演示就直接现算
    for(int u = 0; u < (int)pairNodes.size(); ++u) {
        for(int v = u + 1; v < (int)pairNodes.size(); ++v) {
            // pairNodes[u] = (i1, j1)
            // pairNodes[v] = (i2, j2)
            const PairNode& pu = pairNodes[u];
            const PairNode& pv = pairNodes[v];
            double dSrc = distance3D(srcPoints[pu.i], srcPoints[pv.i]);
            double dTgt = distance3D(tgtPoints[pu.j], tgtPoints[pv.j]);
            if(std::fabs(dSrc - dTgt) < tolerance) {
                // 说明这两个配对点对 在距离上“兼容”
                adjacency[u].push_back(v);
                adjacency[v].push_back(u);
            }
        }
    }
    return adjacency;
}

//--------------------------------------
// 最大团搜索(回溯/Bron–Kerbosch 等方法)
// 这里用一个简单的回溯示例，点数较多时可能效率欠佳
//--------------------------------------
void backtrackMaxClique(
    const std::vector<std::vector<int>>& adjacency,
    std::vector<int>& currentClique,
    std::vector<int>& bestClique,
    const std::vector<int>& candidates,
    int startPos)
{
    // 剪枝条件：剩余候选加上 currentClique 的大小都不可能超过 bestClique，就可以停止
    if((int)currentClique.size() + ((int)candidates.size() - startPos) <= (int)bestClique.size()) {
        return;
    }

    for(int i = startPos; i < (int)candidates.size(); ++i) {
        int v = candidates[i];
        // 检查 v 是否与 currentClique 里的所有点都相连
        bool canJoin = true;
        for(int c : currentClique) {
            // 如果 c 与 v 不相连，无法继续
            // adjacency[c] 应该包含 v
            auto &adjList = adjacency[c];
            if(std::find(adjList.begin(), adjList.end(), v) == adjList.end()) {
                canJoin = false;
                break;
            }
        }
        if(canJoin) {
            // 加入 v 到 currentClique
            currentClique.push_back(v);
            // 继续向下搜索
            backtrackMaxClique(adjacency, currentClique, bestClique, candidates, i + 1);
            // 回溯
            if(currentClique.size() > bestClique.size()) {
                bestClique = currentClique;
            }
            currentClique.pop_back();
        }
    }
}

std::vector<int> findMaxClique(const std::vector<std::vector<int>>& adjacency)
{
    // 简单做法：候选就是所有点 0...(n-1)
    std::vector<int> candidates(adjacency.size());
    for(int i = 0; i < (int)adjacency.size(); ++i){
        candidates[i] = i;
    }

    std::vector<int> bestClique;
    std::vector<int> currentClique;
    backtrackMaxClique(adjacency, currentClique, bestClique, candidates, 0);

    return bestClique;
}

//--------------------------------------
// 主流程示例
//--------------------------------------
int main()
{
    // 假设我们有两个文件，分别存放两组点
    std::string fileA = "pointsA.txt";
    std::string fileB = "pointsB.txt";
    
    // 1) 读入两组点
    auto srcPoints = readPointsFromFile(fileA);
    auto tgtPoints = readPointsFromFile(fileB);
    std::cout << "Read " << srcPoints.size() << " points from " << fileA << std::endl;
    std::cout << "Read " << tgtPoints.size() << " points from " << fileB << std::endl;

    if(srcPoints.empty() || tgtPoints.empty()){
        std::cerr << "One of the point sets is empty, abort." << std::endl;
        return -1;
    }

    // 2) 构造兼容图
    double tolerance = 1.0; // 距离容忍度
    std::vector<PairNode> pairNodes; // (i,j) 列表
    auto adjacency = buildCompatibilityGraph(srcPoints, tgtPoints, pairNodes, tolerance);

    // 3) 寻找最大团
    std::vector<int> maxClique = findMaxClique(adjacency);

    // 4) 输出结果
    std::cout << "Max clique size = " << maxClique.size() << std::endl;
    // 对应点对
    for(int idx : maxClique) {
        const auto& pn = pairNodes[idx];
        const Point3D& spt = srcPoints[pn.i];
        const Point3D& tpt = tgtPoints[pn.j];
        // 计算两点的绝对误差（只做一个距离比较，也可以单独输出每对点坐标的差）
        double err = std::fabs(distance3D(spt, spt) - distance3D(tpt, tpt)); 
        // 事实上 distance3D(spt, spt) = 0, distance3D(tpt, tpt) = 0，这里仅演示，可以改为其它误差度量
        // 通常更有意义的是把 S_i 到 T_j 的坐标差或其他误差算出来
        // 这里可以改成：err = distance3D(spt, tpt)，看它们相差多大（若你想看对应点本身是否也相近等等）
        // 或者仅打印出 i, j 即可。
        
        // 演示：直接打印配对以及它们各自坐标
        std::cout << "Source point [" << pn.i << "] = (" 
                  << spt.x << ", " << spt.y << ", " << spt.z << ")  "
                  << "Target point [" << pn.j << "] = (" 
                  << tpt.x << ", " << tpt.y << ", " << tpt.z << ")  "
                  << "Error(demo) = " << err << std::endl;
    }

    return 0;
}

代码要点说明

1. 数据结构

   • Point3D 用来存储点的坐标。
   • PairNode (i, j) 代表第一组点里下标为 i 的点、与第二组点里下标为 j 的点构成的一对。

2. 构造“兼容图”

   • 我们有 N*M 个“候选节点”，遍历两两之间 $(S_{i1}, T_{j1})$ 和 $(S_{i2}, T_{j2})$：
     • 若 $\bigl|\|S_{i1}-S_{i2}\| - \|T_{j1}-T_{j2}\|\bigr| < \text{tolerance}$，则它们之间连边。
   • 这样就得到了邻接表 adjacency。

3. 最大团搜索

   • 代码中使用了一个简单的回溯来搜索最大团。
   • 若点的数量更大，需要更高效的搜索或剪枝(例如 Bron–Kerbosch 算法、或引用一些现成的最大团/最大匹配库)。

4. 输出结果

   • 在找到的最大团 maxClique 中，每个元素是一个索引 idx，对应 pairNodes[idx] = (i, j)。
   • 这就表示：源点集中第 i 个点，与目标点集中第 j 个点，存在一致的距离关系，从而匹配上了。
   • 你可以根据实际需要，打印坐标、计算对应点之间的误差等。

如果你想像示例结果那样，输出每一对匹配点的坐标、并计算它们之间的空间误差(例如直接看它们的欧式距离是否也在一定范围内)，可在最终输出的时候做更多计算。

备注

• 上述方法是基于“所有点两两间距都要匹配”，才能称为同一个公共子图(对距离容忍度内)。在点数较少、或对精度要求较高的情况下，这种方法可行。如果点很多，计算量会随 $N^2 \times M^2$ 增加，需要一定的算法优化(如预计算距离、适当剪枝等)。
• 如果你需要真正去做“刚体变换”下的最优配准(比如 ICP、RANSAC、或基于 SIFT-like 特征的 3D 特征匹配)，则会是另一种思路。但题目里所述是基于“最大公共子图”的距离关系匹配，所以就按照图论思路来做。
• 在示例输出里，你也可以进一步计算“Source point - Target point”的距离作为“Error”，而不是像上面示范里用 distance3D(spt, spt) 这类占位写法。具体请根据你要衡量的“误差定义”作调整。

这样，你就可以根据这段代码，针对两组 3D 点(文本文件输入)，构建兼容图并找出最大公共子图对应的匹配关系了。祝开发顺利!