为什么我用GDAL矢量求交很慢

本来计算16000条道路和自身（16000条道路）交点，随着程序运行，C盘空间会逐渐减少，最后会因为内存不足运行中止，故只能改为1000条道路和16000条道路的交点。这还要话15分钟。

#include <stdio.h>
#include <iostream>
#include<string>
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
using namespace std;
void main()
{
    //———————————————————————打开文件——————————————————————————
    //读取矢量
    //为支持中文路径
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
    string strFilePath1 = "D:\wuhan_osm_ln.shp";
    //注册所有驱动
    GDALAllRegister();
    GDALDataset* poSrcDS = (GDALDataset*)GDALOpenEx(strFilePath1.c_str(), GA_Update, NULL, NULL, NULL);
    if (poSrcDS == NULL)
    {
        cout << "源文件打开失败" << endl;
        return;
    }
    cout << "文件打开成功" << endl;
    OGRLayer* poSrcLayer= poSrcDS->GetLayer(0);//打开图层
    if (poSrcLayer == NULL)
    {
        cout << "打开失败" << endl;
        return;
    }
    //———————————————————————创建交点shp文件————————————————————
    const char *pszDriverName = "ESRI Shapefile";
    //调用对Shape文件读写的Driver
    GDALDriver *poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName);
    if (poDriver == NULL)
    {
        cout << pszDriverName << "驱动不可用！" << endl;
        return;
    }
    //创建数据源
    string stroutputPath = "D:\intersetPoint0-999.shp";
    GDALDataset *poDstDs = poDriver->Create(stroutputPath.c_str(), 0, 0, 0, GDT_Byte, NULL);
    if (poDstDs == NULL)
    {
        cout << "DataSource Creation Error" << endl;
        return;
    }
    //创建图层Layer
    //参数说明：新图层名称，坐标系，图层的几何类型，创建选项，与驱动有关
    OGRSpatialReference spatialReference;
    spatialReference.SetWellKnownGeogCS("WGS84");
    OGRLayer* poDstLayer = poDstDs->CreateLayer("intersetPoint0-999", &spatialReference, wkbPoint, NULL);
    //先创建一个“ID”整型属性
    OGRFieldDefn oFieldId("Interset1", OFTInteger);
    oFieldId.SetWidth(7);
    OGRFieldDefn oFieldId2("InTerset2", OFTInteger);
    oFieldId.SetWidth(7);
    poDstLayer->CreateField(&oFieldId);
    poDstLayer->CreateField(&oFieldId2);
    //———————————————————————————————————————————————————
    OGRFeature *poDstFeature = NULL;
    OGRGeometry* poDstGeometry = NULL;//16552
    int nfeatureCount = poSrcLayer->GetFeatureCount();//获得要素数量
    //遍历要素求交
//本来n为nfeatureCount，但是因为跑的太慢，而且越跑C盘剩余空间越小，会因为C盘占满，内存不足，代码中止。
    for (int n = 0; n < 999; n++)
    {
        OGRFeature*poSrcFeature = poSrcLayer->GetFeature(n);//获得要素
        OGRGeometry*poSrcGeometry = poSrcFeature->GetGeometryRef();//获得该要素几何体
        for (int m = 0;  m <nfeatureCount; m++)
        {
            if (m == n)continue;//跳过相同线

            OGRFeature*poOtherFeature = poSrcLayer->GetFeature(m);//获得另一要素
            OGRGeometry*poOtherGeometry = poOtherFeature->GetGeometryRef();//获得另一要素几何信息
            poDstGeometry = poSrcGeometry->Intersection(poOtherGeometry);//获取相交几何信息
            if (poDstGeometry!=NULL&&(!poDstGeometry->IsEmpty()))
            {
                poDstFeature = new OGRFeature(poDstLayer->GetLayerDefn());//GetLayerDefn()获取当前图层的属性表结构，新建要素
                poDstFeature->SetField("Interset1", poSrcFeature->GetFieldAsInteger("OBJECTID"));//第一个相交的线ID；
                poDstFeature->SetField("Interset2", poOtherFeature->GetFieldAsInteger("OBJECTID"));//第二个相交的线ID；
                poDstFeature->SetGeometry(poDstGeometry);//设置几何信息
                poDstLayer->CreateFeature(poDstFeature);//设置几何要素
            }
        }
    }


}