绘制等高线图的算法

最新推荐文章于 2024-04-29 20:24:00 发布

还是叫明

最新推荐文章于 2024-04-29 20:24:00 发布

阅读量1.7w

点赞数 10

分类专栏：算法文章标签：等高线图 Marching squares 算法

本文链接：https://blog.csdn.net/lweiyue/article/details/91490460

版权

算法专栏收录该内容

23 篇文章 18 订阅

订阅专栏

假设我们在一个矩形范围内，有一些离散的高度数据，形如(x,y,height)这样的数据集，然后我们要得到一张这样的等高线图：

我们描述一下生成等高线图的算法。

一、图形概述

实际地图上实用的单位是米，但我们显示时使用的单位是像素，这里面有一个转换的关系。而且，显示的图形有可能需要缩放。所以我们收集的数据，x和y统一使用相对值。我们把总长和总宽都定为1，实际的坐标相对于单位1来定。例如地图长8km，然后我们的坐标点x为1km，那么x就是1/8=0.125。

实际上，等高线图是分成一个一个网格的，我所使用的网格大小是10*10。每个网格里的线和颜色组成了整张的等高线图。

二、插值

我们的数据都是形如(x,y,height)这样的格式，如(0.85f, 0.05f, 0.49f)、(0.75f, 0.1f, 0.7f)、(0.5f, 0.65f, 0.04f)。

假如我们要显示的图形大小是800*600，切成网格后，就会有80*60=4800个网格。一般情况下，我们收集到的数据量远远小于这个数，所以我们就需要插值。

我们用比较简单的距离反比策略进行插值。我们知道，跟某个点的距离越近，受它的影响就越大。所以，跟点V(x1,y1)相距d的点P(x,y)，它受V的影响可按以下公式计算：

P-V=V/d

其中d=(x-x1)^2+(y-y1)^2

而P的最终值，由平面上所有已知值的点确定。

代码如下：

double D = 0;
double DV = 0;
foreach (IntMeasureData imd in measure_data)
{
    double d = 1.0 / ((imd.X - i) * (imd.X - i) + (imd.Y - j) * (imd.Y - j));
    D += d;
    DV += imd.Z * d;
}
value = (float)(DV / D);

三、Marching squares算法

这就是生成等高线图的一种算法，我们来描述一下。

我们上面提到过，把800*600的图切成了80*60个网格，而在第二步里面，我们把这些网格都填满数据了。我们假设数据如下：

图中的数据就是高度。我们现在要画一条高度值为5的线。我们把5跟表格中的数据比较，比5大的为1，用实心黑点表示，比5小的为0，化为一个二值的表格，如下图所示：

然后如果我们把每个格子取出来，就会发现有以下的16种情况：

只要我们把上图的线画出来，就会拼出一张等高线图。

四、算法细节

问题1，我们可能需要的是有颜色的图，而不是只有线条的图。

我们只要把上图中格子里的线和实心黑点围成的区域填充颜色即可。

问题2，不同的高度值，如何绘图。

先画低阈值的，再画高阈值的，一层层往上叠即可。

问题3，有没有办法把所有格子的线条围成一个区域，然后对这个区域进行一次性颜色的填充。

如果有，就不会用这种分成一格一格的方法去做了。只能一格一格地填充颜色。

五、代码解析

代码是C#、WPF。

我们先定义测量数据的数据结构：

/// <summary>
/// 测量数据
/// </summary>
public struct MeasureData
{
    /// <summary>
    /// 初始化测量数据
    /// </summary>
    /// <param name="x">坐标x</param>
    /// <param name="y">坐标y</param>
    /// <param name="z">高度</param>
    public MeasureData(float x, float y, float z)
    {
        X = x;
        Y = y;
        Z = z;
    }

    /// <summary>
    /// 坐标X
    /// </summary>
    public float X;

    /// <summary>
    /// 坐标Y
    /// </summary>
    public float Y;

    /// <summary>
    /// 高度
    /// </summary>
    public float Z;
}

internal struct IntMeasureData
{
    public IntMeasureData(MeasureData md, int x_num, int y_num)
    {
        X = (int)(md.X * x_num);
        if (X >= x_num)
        {
            X = x_num - 1;
        }
        Y = (int)(md.Y * y_num);
        if (Y >= y_num)
        {
            Y = y_num - 1;
        }
        Z = md.Z;
    }

    public int X;
    public int Y;
    public float Z;
}

IntMeasureData是存放在网格里的数据。

然后我们做第一步：数据插值：

/// <summary>
/// 数据插值
/// </summary>
private void InitData()
{
    x_num = (int)(w / grid_w);
    y_num = (int)(h / grid_h);
    data = new float[x_num, y_num];

    IntMeasureData[] measure_data = new IntMeasureData[HeightDots.Length];
    for (int i = HeightDots.Length - 1; i >= 0; i--)
    {
        measure_data[i] = new IntMeasureData(HeightDots[i], x_num, y_num);
    }

    min = float.MaxValue;
    max = float.MinValue;
    for (int i = 0; i < x_num; i++)
    {
        for (int j = 0; j < y_num; j++)
        {
            float value = 0;
            bool find = false;
            foreach (IntMeasureData imd in measure_data)
            {
                if (i == imd.X && j == imd.Y)
                {
                    value = imd.Z;
                    find = true;
                    break;
                }
            }

            if (!find)
            {
                double D = 0;
                double DV = 0;
                foreach (IntMeasureData imd in measure_data)
                {
                    double d = 1.0 / ((imd.X - i) * (imd.X - i) + (imd.Y - j) * (imd.Y - j));
                    D += d;
                    DV += imd.Z * d;
                }
                value = (float)(DV / D);
            }

            data[i, j] = value;
            min = Math.Min(min, value);
            max = Math.Max(max, value);
        }
    }
}

核心是Marching squares算法的实现：

/// <summary>
/// 获取某个阈值下的图形数据
/// </summary>
/// <param name="threshold">阈值</param>
/// <returns>图形数据</returns>
private List<List<PointF>> CreateMapData(float threshold)
{
    byte[,] binary_data = new byte[x_num, y_num];
    for (int i = 0; i < x_num; i++)
    {
        for (int j = 0; j < y_num; j++)
        {
            binary_data[i, j] = (byte)((data[i, j] >= threshold) ? 1 : 0);
        }
    }

    List<List<PointF>> shapes = new List<List<PointF>>();
    for (int i = 1; i < x_num; i++)
    {
        for (int j = 1; j < y_num; j++)
        {
            int num = (binary_data[i - 1, j - 1] << 3) + (binary_data[i, j - 1] << 2) + (binary_data[i, j] << 1) + binary_data[i - 1, j];
            List<PointF> points = new List<PointF>();
            switch (num)
            {
                case 0:
                    break;
                case 1:
                    AddLeft(points, i, j, threshold);
                    AddLeftBottom(points, i, j);
                    AddBottom(points, i, j, threshold);
                    break;
                case 2:
                    AddRight(points, i, j, threshold);
                    AddRightBottom(points, i, j);
                    AddBottom(points, i, j, threshold);
                    break;
                case 3:
                    AddLeft(points, i, j, threshold);
                    AddLeftBottom(points, i, j);
                    AddRightBottom(points, i, j);
                    AddRight(points, i, j, threshold);
                    break;
                case 4:
                    AddTop(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddRight(points, i, j, threshold);
                    break;
                case 5:
                    AddLeft(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddRight(points, i, j, threshold);
                    AddBottom(points, i, j, threshold);
                    AddLeftBottom(points, i, j);
                    AddLeft(points, i, j, threshold);
                    break;
                case 6:
                    AddTop(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddRightBottom(points, i, j);
                    AddBottom(points, i, j, threshold);
                    break;
                case 7:
                    AddTop(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddRightBottom(points, i, j);
                    AddLeftBottom(points, i, j);
                    AddLeft(points, i, j, threshold);
                    break;
                case 8:
                    AddTop(points, i, j, threshold);
                    AddLeftTop(points, i, j);
                    AddLeft(points, i, j, threshold);
                    break;
                case 9:
                    AddTop(points, i, j, threshold);
                    AddLeftTop(points, i, j);
                    AddLeftBottom(points, i, j);
                    AddBottom(points, i, j, threshold);
                    break;
                case 10:
                    AddTop(points, i, j, threshold);
                    AddLeftTop(points, i, j);
                    AddLeft(points, i, j, threshold);
                    AddBottom(points, i, j, threshold);
                    AddRightBottom(points, i, j);
                    AddRight(points, i, j, threshold);
                    break;
                case 11:
                    AddTop(points, i, j, threshold);
                    AddLeftTop(points, i, j);
                    AddLeftBottom(points, i, j);
                    AddRightBottom(points, i, j);
                    AddRight(points, i, j, threshold);
                    break;
                case 12:
                    AddLeft(points, i, j, threshold);
                    AddRight(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddLeftTop(points, i, j);
                    break;
                case 13:
                    AddRight(points, i, j, threshold);
                    AddRightTop(points, i, j);
                    AddLeftTop(points, i, j);
                    AddLeftBottom(points, i, j);
                    AddBottom(points, i, j, threshold);
                    break;
                case 14:
                    AddLeft(points, i, j, threshold);
                    AddBottom(points, i, j, threshold);
                    AddRightBottom(points, i, j);
                    AddRightTop(points, i, j);
                    AddLeftTop(points, i, j);
                    break;
                case 15:
                    AddLeftTop(points, i, j);
                    AddRightTop(points, i, j);
                    AddRightBottom(points, i, j);
                    AddLeftBottom(points, i, j);
                    break;
            }
            if (num != 0)
            {
                shapes.Add(points);
            }
        }
    }

    return shapes;
}

private void AddLeft(List<PointF> list, int x, int y, float threshold)
{
    list.Add(new PointF((x - 1) * grid_w, (y - 1 + V(data[x - 1, y - 1], data[x - 1, y], threshold)) * grid_h));
}

private void AddRight(List<PointF> list, int x, int y, float threshold)
{
    list.Add(new PointF(x * grid_w, (y - 1 + V(data[x, y - 1], data[x, y], threshold)) * grid_h));
}

private void AddTop(List<PointF> list, int x, int y, float threshold)
{
    list.Add(new PointF((x - 1 + V(data[x - 1, y - 1], data[x, y - 1], threshold)) * grid_w, (y - 1) * grid_h));
}

private void AddBottom(List<PointF> list, int x, int y, float threshold)
{
    list.Add(new PointF((x - 1 + V(data[x - 1, y], data[x, y], threshold)) * grid_w, y * grid_h));
}

private void AddLeftTop(List<PointF> list, int x, int y)
{
    list.Add(new PointF((x - 1) * grid_w, (y - 1) * grid_h));
}

private void AddRightTop(List<PointF> list, int x, int y)
{
    list.Add(new PointF(x * grid_w, (y - 1) * grid_h));
}

private void AddLeftBottom(List<PointF> list, int x, int y)
{
    list.Add(new PointF((x - 1) * grid_w, y * grid_h));
}

private void AddRightBottom(List<PointF> list, int x, int y)
{
    list.Add(new PointF(x * grid_w, y * grid_h));
}

private float V(float min, float max, float threshold)
{
    return (threshold - min) / (max - min);
}

结果是我们得到了每个网格要填充颜色的端点列表。把这些端点围成一个多边形，再绘制成图形即可。

完整代码下载

还是叫明

关注

10
点赞
踩
91

收藏

觉得还不错? 一键收藏
17
评论
绘制等高线图的算法

假设我们在一个矩形范围内，有一些离散的高度数据，形如(x,y,height)这样的数据集，然后我们要得到一张这样的等高线图：我们描述一下生成等高线图的算法。一、图形概述实际地图上实用的单位是米，但我们显示时使用的单位是像素，这里面有一个转换的关系。而且，显示的图形有可能需要缩放。所以我们收集的数据，x和y统一使用相对值。我们把总长和总宽都定为1，实际的坐标相对于单位1来定。例如地图...
复制链接

扫一扫