osgearth如何实现heatmap热力图

liyihongcug

数据：
var points = [
        {"lng":116.418261, "lat" : 39.921984, "count" : 50},
        { "lng":116.423332,"lat" : 39.916532,"count" : 51 },
        { "lng":116.419787,"lat" : 39.930658,"count" : 15 },
        { "lng":116.418455,"lat" : 39.920921,"count" : 40 },
        { "lng":116.418843,"lat" : 39.915516,"count" : 100 },
        { "lng":116.42546,"lat" : 39.918503,"count" : 6 },
        { "lng":116.423289,"lat" : 39.919989,"count" : 18 },
        { "lng":116.418162,"lat" : 39.915051,"count" : 80 },
        { "lng":116.422039,"lat" : 39.91782,"count" : 11 },
        { "lng":116.41387,"lat" : 39.917253,"count" : 7 },
        { "lng":116.41773,"lat" : 39.919426,"count" : 42 },
        { "lng":116.421107,"lat" : 39.916445,"count" : 4 },
        { "lng":116.417521,"lat" : 39.917943,"count" : 27 },
        { "lng":116.419812,"lat" : 39.920836,"count" : 23 },
        { "lng":116.420682,"lat" : 39.91463,"count" : 60 },
        { "lng":116.415424,"lat" : 39.924675,"count" : 8 },
        { "lng":116.419242,"lat" : 39.914509,"count" : 15 },
        { "lng":116.422766,"lat" : 39.921408,"count" : 25 },
        { "lng":116.421674,"lat" : 39.924396,"count" : 21 },
        { "lng":116.427268,"lat" : 39.92267,"count" : 1 },
        { "lng":116.417721,"lat" : 39.920034,"count" : 51 },
        { "lng":116.412456,"lat" : 39.92667,"count" : 7 },
        { "lng":116.420432,"lat" : 39.919114,"count" : 11 },
        { "lng":116.425013,"lat" : 39.921611,"count" : 35 },
        { "lng":116.418733,"lat" : 39.931037,"count" : 22 },
        { "lng":116.419336,"lat" : 39.931134,"count" : 4 },
        { "lng":116.413557,"lat" : 39.923254,"count" : 5 },
        { "lng":116.418367,"lat" : 39.92943,"count" : 3 },
        { "lng":116.424312,"lat" : 39.919621,"count" : 100 },
        { "lng":116.423874,"lat" : 39.919447,"count" : 87 },
        { "lng":116.424225,"lat" : 39.923091,"count" : 32 },
        { "lng":116.417801,"lat" : 39.921854,"count" : 44 },
        { "lng":116.417129,"lat" : 39.928227,"count" : 21 },
        { "lng":116.426426,"lat" : 39.922286,"count" : 80 },
        { "lng":116.421597,"lat" : 39.91948,"count" : 32 },
        { "lng":116.423895,"lat" : 39.920787,"count" : 26 },
        { "lng":116.423563,"lat" : 39.921197,"count" : 17 },
        { "lng":116.417982,"lat" : 39.922547,"count" : 17 },
        { "lng":116.426126,"lat" : 39.921938,"count" : 25 },
        { "lng":116.42326,"lat" : 39.915782,"count" : 100 },
        { "lng":116.419239,"lat" : 39.916759,"count" : 39 },
        { "lng":116.417185,"lat" : 39.929123,"count" : 11 },
        { "lng":116.417237,"lat" : 39.927518,"count" : 9 },
        { "lng":116.417784,"lat" : 39.915754,"count" : 47 },
        { "lng":116.420193,"lat" : 39.917061,"count" : 52 },
        { "lng":116.422735,"lat" : 39.915619,"count" : 100 },
        { "lng":116.418495,"lat" : 39.915958,"count" : 46 },
        { "lng":116.416292,"lat" : 39.931166,"count" : 9 },
        { "lng":116.419916,"lat" : 39.924055,"count" : 8 },
        { "lng":116.42189,"lat" : 39.921308,"count" : 11 },
        { "lng":116.413765,"lat" : 39.929376,"count" : 3 },
        { "lng":116.418232,"lat" : 39.920348,"count" : 50 },
        { "lng":116.417554,"lat" : 39.930511,"count" : 15 },
        { "lng":116.418568,"lat" : 39.918161,"count" : 23 },
        { "lng":116.413461,"lat" : 39.926306,"count" : 3 },
        { "lng":116.42232,"lat" : 39.92161,"count" : 13 },
        { "lng":116.4174,"lat" : 39.928616,"count" : 6 },
        { "lng":116.424679,"lat" : 39.915499,"count" : 21 },
        { "lng":116.42171,"lat" : 39.915738,"count" : 29 },
        { "lng":116.417836,"lat" : 39.916998,"count" : 99 },
        { "lng":116.420755,"lat" : 39.928001,"count" : 10 },
        { "lng":116.414077,"lat" : 39.930655,"count" : 14 },
        { "lng":116.426092,"lat" : 39.922995,"count" : 16 },
        { "lng":116.41535,"lat" : 39.931054,"count" : 15 },
        { "lng":116.413022,"lat" : 39.921895,"count" : 13 },
        { "lng":116.415551,"lat" : 39.913373,"count" : 17 },
        { "lng":116.421191,"lat" : 39.926572,"count" : 1 },
        { "lng":116.419612,"lat" : 39.917119,"count" : 9 },
        { "lng":116.418237,"lat" : 39.921337,"count" : 54 },
        { "lng":116.423776,"lat" : 39.921919,"count" : 26 },
        { "lng":116.417694,"lat" : 39.92536,"count" : 17 },
        { "lng":116.415377,"lat" : 39.914137,"count" : 19 },
        { "lng":116.417434,"lat" : 39.914394,"count" : 43 },
        { "lng":116.42588,"lat" : 39.922622,"count" : 27 },
        { "lng":116.418345,"lat" : 39.919467,"count" : 8 },
        { "lng":116.426883,"lat" : 39.917171,"count" : 3 },
        { "lng":116.423877,"lat" : 39.916659,"count" : 34 },
        { "lng":116.415712,"lat" : 39.915613,"count" : 14 },
        { "lng":116.419869,"lat" : 39.931416,"count" : 12 },
        { "lng":116.416956,"lat" : 39.925377,"count" : 11 },
        { "lng":116.42066,"lat" : 39.925017,"count" : 38 },
        { "lng":116.416244,"lat" : 39.920215,"count" : 91 },
        { "lng":116.41929,"lat" : 39.915908,"count" : 54 },
        { "lng":116.422116,"lat" : 39.919658,"count" : 21 },
        { "lng":116.4183,"lat" : 39.925015,"count" : 15 },
        { "lng":116.421969,"lat" : 39.913527,"count" : 3 },
        { "lng":116.422936,"lat" : 39.921854,"count" : 24 },
        { "lng":116.41905,"lat" : 39.929217,"count" : 12 },
        { "lng":116.424579,"lat" : 39.914987,"count" : 57 },
        { "lng":116.42076,"lat" : 39.915251,"count" : 70 },
        { "lng":116.425867,"lat" : 39.918989,"count" : 8 }];

liyihongcug

方法复杂main.cpp   osg::ref_ptr<osg::Node> dt = osgDB::readRefNodeFile("dumptruck.osg");

    // display a solid version of the dump truck
    osg::ref_ptr<osg:ositionAttitudeTransform> solidModel = new osg::PositionAttitudeTransform;
    solidModel->setPosition(osg::Vec3f(7.0f, -2.0f, 7.0f));
    solidModel->addChild(dt);

    // generate the 3D heatmap surface to display
    osg::ref_ptr<Heatmap> hm = new Heatmap(30, 30, 10, 30, 30, 1.0, 0.25);
    float data[30][30];
    for (int x=0; x < 30; ++x)
        for (int y=0; y < 30; ++y)
            data[y][x] = (double)rand() / RAND_MAX;
    hm->setData((float*)data, 10.0, 1.0, 0.25);

    // add a transparent version of the truck to the scene also
    osg::ref_ptr<osg::PositionAttitudeTransform> transparentTruck = new osg::PositionAttitudeTransform;
    transparentTruck->setPosition(osg::Vec3f(7.0f, -25.0f, 7.0f));

    // set the states of the truck so that it actually appears transparently and nicely lit.
    osg::StateSet *state = transparentTruck->getOrCreateStateSet();
    state->setMode(GL_BLEND, osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE);
    state->setAttribute(new osg::BlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA), osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE);
    osg::Material* material = new osg::Material;
    material->setAmbient(osg::Material::FRONT_AND_BACK,osg::Vec4(0.2f,0.2f,0.2f,0.3f));
    material->setDiffuse(osg::Material::FRONT_AND_BACK,osg::Vec4(0.8f,0.8f,0.8f,0.3f));
    state->setAttribute(material,osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE);
    osg:ightModel *lm = new osg::LightModel();
    lm->setTwoSided(true);
    state->setAttribute(lm, osg::StateAttribute::ON | osg::StateAttribute::OVERRIDE);
    (transparentTruck.get())->addChild(dt.get());

    // place the heatmap and a transparent dump truck in the transparent geometry group
    osg::ref_ptr<osg::Group> transparentModel = new osg::Group;
    (transparentModel.get())->addChild(hm.get());
    (transparentModel.get())->addChild(transparentTruck.get());

liyihongcug

#ifndef HEATMAP_H
#define HEATMAP_H

#include <osg/Geode>
#include <osg/Uniform>
#include <osg/Texture2D>
#include <osg/Texture1D>

class Heatmap : public osg::Geode
{
public:
    Heatmap(float width, float depth, float maxheight, unsigned int K, unsigned int N, float maximum, float transparency);
    ~Heatmap();

    void setData(float *buffer, float maxheight, float maximum, float transparency);

protected:
    unsigned int m_K;
    unsigned int m_N;
    float *m_data;
    osg::ref_ptr<osg::Image> m_img2;
    osg::ref_ptr<osg::Texture2D> m_tex2;

    osg::ref_ptr<osg::Image> colorimg;
    osg::ref_ptr<osg::Texture1D> colortex;

    osg::Uniform *maximumUniform;
    osg::Uniform *maxheightUniform;
    osg::Uniform *transparencyUniform;
};

#endif // #ifndef HEATMAP_H



/*
* 3D Heat map using vertex displacement mapping
* Rendered using depth peeling in fragment shader.
*/

#include <osg/Geode>
#include <osg/Geometry>
#include <osg/Vec3>
#include <osg/MatrixTransform>
#include <osg/PositionAttitudeTransform>
#include <osg/LightModel>
#include <osg/io_utils>
#include <osg/Material>
#include <osgDB/ReadFile>
#include <osgDB/WriteFile>
#include <osgGA/TrackballManipulator>
#include <osgViewer/Viewer>
#include <osg/Math>
#include <iostream>

#define _USE_MATH_DEFINES
#include <math.h>

#include <osg/TexEnv>
#include <osg/TexMat>
#include <osg/Depth>
#include <osg/ShapeDrawable>
#include <osg/Texture1D>
#include <osg/Texture2D>
#include <osg/PolygonMode>
#include <osg/PolygonOffset>

#include "HeatMap.h"
#include "DepthPeeling.h"

///////////////////////////////////////////////////////////////////////////
// in-line GLSL source code

static const char *VertexShader = {
    "#version 120\n"
    "uniform float maximum;\n"
    "uniform float maxheight;\n"
    "uniform float transparency;\n"
    "uniform sampler1D colortex;\n"
    "uniform sampler2D datatex;\n"
    "in vec2  xypos;\n"
    "void main(void)\n"
    "{\n"
    "    float foo;\n"
    "    float tmp     = min(texture2D(datatex, xypos).x / maximum, 1.0);\n"
    "    gl_Position   = gl_ModelViewProjectionMatrix * (gl_Vertex + vec4(0.0, 0.0, maxheight * tmp, 0.0));\n"
    "    vec4 color    = texture1D(colortex, tmp);\n"
    "    color.w       = color.w * transparency;\n"
    "    gl_FrontColor = color;\n"
    "    gl_BackColor  = color;\n"
    "}\n"
};

static const char *FragmentShader =
{
    "#version 120\n"
    "bool depthpeeling();\n"
    "void main(void)\n"
    "{\n"
    "  if( depthpeeling() ) discard;\n"
    "  gl_FragColor = gl_Color;\n"
    "}\n"
};

/**
* Overloaded Geometry class to return predefined bounds
*/
class MyGeometry : public osg::Geometry
{
public:
    MyGeometry(osg::BoundingBox bounds)
    {
        m_bounds = bounds;
        m_bsphere = osg::BoundingSphere(bounds);
    }

    // an attempt to return a reasonable bounding box. Still does not prevent clipping of the heat map.
    virtual const osg::BoundingBox& getBoundingBox() const {return m_bounds;}
    virtual osg::BoundingBox computeBoundingBox() const {return m_bounds;}
    virtual const osg::BoundingSphere& getBound() const {return m_bsphere;}

protected:
    osg::BoundingBox m_bounds;
    osg::BoundingSphere m_bsphere;
};

Heatmap::Heatmap(float width, float depth, float maxheight, unsigned int K, unsigned int N, float maximum, float transparency)
{
    m_K = K;
    m_N = N;
    const int O = 4;

    // create Geometry object to store all the vertices primitives.
    osg::Geometry *meshGeom  = new MyGeometry(osg::BoundingBox(osg::Vec3(-width/2, -depth/2, 0), osg::Vec3(width/2, depth/2, maxheight)));

    // we use a float attribute array storing texcoords
    osg::Vec2Array* xypositions = new osg::Vec2Array();
    xypositions->setName("xypos");

    // create vertex coordinates
    osg::Vec3Array* vertices = new osg::Vec3Array();
    osg::Vec3 off(-width/2, -depth/2, 0);
    for (unsigned int y=0; y < O*N; y++) {
        if (y % 2 == 0)
        {
            for (unsigned int x=0; x < O*K; x++) {
                vertices->push_back(osg::Vec3(width*x/(O*K-1), depth*y/(O*N-1), 0.0)+off);
                xypositions->push_back(osg::Vec2(((float)x+0.5f)/(O*K),((float)y+0.5f)/(O*N)));
            }
        }
        else
        {
            vertices->push_back(osg::Vec3(0, depth*y/(O*N-1), 0.0)+off);
            xypositions->push_back(osg::Vec2(0.5f/(O*K),((float)y+0.5f)/(O*N)));
            for (unsigned int x=0; x < O*K-1; x++) {
                vertices->push_back(osg::Vec3(width*(0.5+x)/(O*K-1), depth*y/(O*N-1), 0.0)+off);
                xypositions->push_back(osg::Vec2(((float)(x+0.5f)+0.5f)/(O*K),((float)y+0.5f)/(O*N)));
            }
            vertices->push_back(osg::Vec3(width, depth*y/(O*N-1), 0.0)+off);
            xypositions->push_back(osg::Vec2(1.0f-0.5f/(O*K),((float)y+0.5f)/(O*N)));
        }
    }
    xypositions->setBinding(osg::Array::BIND_PER_VERTEX);
    xypositions->setNormalize(false);

    meshGeom->setVertexAttribArray(6, xypositions);
    meshGeom->setVertexArray(vertices);

    // generate several tri strips to form a mesh
    GLuint *indices = new GLuint[4*O*K];
    for (unsigned int y=0; y < O*N-1; y++) {
        if (y % 2 == 0)
        {
            int base  = (y/2) * (O*K+O*K+1);
            int base2 = (y/2) * (O*K+O*K+1) + O*K;
            int i=0; for (unsigned int x=0; x < O*K; x++) { indices[i++] = base2+x; indices[i++] = base+x;}
            indices[i++] = base2+O*K;
            meshGeom->addPrimitiveSet(new osg:rawElementsUInt(osg:rimitiveSet::TRIANGLE_STRIP, i, indices));
        }
        else
        {
            int base = (y/2) * (O*K+O*K+1) + O*K;
            int base2 = (y/2) * (O*K+O*K+1) + O*K + O*K+1;
            int i=0; for (unsigned int x=0; x < O*K; x++) { indices[i++] = base+x; indices[i++] = base2+x;}
            indices[i++] = base+O*K;
            meshGeom->addPrimitiveSet(new osg::DrawElementsUInt(osg::PrimitiveSet::TRIANGLE_STRIP, i, indices));
        }
    }
    delete[] indices;

    // create vertex and fragment shader
    osg::Program* program = new osg::Program;
    program->setName( "mesh" );
    program->addBindAttribLocation("xypos", 6);
    program->addShader( new osg::Shader( osg::Shader::VERTEX, VertexShader ) );
    program->addShader( new osg::Shader( osg::Shader::FRAGMENT, DepthPeeling::PeelingShader ) );
    program->addShader( new osg::Shader( osg::Shader::FRAGMENT, FragmentShader ) );

    // create a 1D texture for color lookups
    colorimg = new osg::Image();
    colorimg->allocateImage(5, 1, 1, GL_RGBA, GL_UNSIGNED_BYTE);
    unsigned char *data = colorimg->data();
    *data++ =   0; *data++ =   0; *data++ = 255; *data++ =   0;  // fully transparent blue
    *data++ =   0; *data++ = 255; *data++ = 255; *data++ = 255;  // turquoise
    *data++ =   0; *data++ = 255; *data++ =   0; *data++ = 255;  // green
    *data++ = 255; *data++ = 255; *data++ =   0; *data++ = 255;  // yellow
    *data++ = 255; *data++ =   0; *data++ =   0; *data++ = 255;  // red
    colortex = new osg::Texture1D(colorimg.get());
    colortex->setFilter(osg::Texture::MIN_FILTER, osg::Texture:INEAR);
    colortex->setFilter(osg::Texture::MAG_FILTER, osg::Texture::LINEAR);
    colortex->setWrap(osg::Texture::WRAP_S, osg::Texture::CLAMP_TO_EDGE);
    colortex->setResizeNonPowerOfTwoHint(false);

    // create a 2D texture for data lookups
    m_img2 = new osg::Image();
    m_img2->allocateImage(K, N, 1, GL_LUMINANCE, GL_FLOAT);
    m_img2->setInternalTextureFormat(GL_RGB32F_ARB);
    m_data = (float*)m_img2.get()->data();
    m_tex2 = new osg::Texture2D(m_img2.get());
    m_tex2.get()->setResizeNonPowerOfTwoHint(false);
    m_tex2.get()->setFilter(osg::Texture::MIN_FILTER, osg::Texture::LINEAR);
    m_tex2.get()->setFilter(osg::Texture::MAG_FILTER, osg::Texture::LINEAR);
    m_tex2.get()->setWrap(osg::Texture::WRAP_S, osg::Texture::CLAMP_TO_EDGE);
    m_tex2.get()->setWrap(osg::Texture::WRAP_T, osg::Texture::CLAMP_TO_EDGE);

    // set render states
    osg::StateSet *meshstate = meshGeom->getOrCreateStateSet();
    meshstate->setMode(GL_BLEND,  osg::StateAttribute::ON);
    meshstate->setMode(GL_LIGHTING, osg::StateAttribute::OFF);
    meshstate->setAttributeAndModes(program, osg::StateAttribute::ON);
    meshstate->setTextureAttributeAndModes(0,colortex.get(),osg::StateAttribute::ON);
    meshstate->setTextureAttributeAndModes(1,m_tex2.get(),osg::StateAttribute::ON);

    // uniforms for height and color scaling
    maximumUniform = new osg::Uniform( "maximum", (float)maximum );
    maxheightUniform = new osg::Uniform( "maxheight", (float)maxheight );
    transparencyUniform = new osg::Uniform( "transparency", (float)transparency);

    osg::Uniform* texUniform = new osg::Uniform(osg::Uniform::SAMPLER_1D, "colortex");
    texUniform->set(0);
    osg::Uniform* texUniform2 = new osg::Uniform(osg::Uniform::SAMPLER_2D, "datatex");
    texUniform2->set(1);
    meshstate->addUniform(texUniform);
    meshstate->addUniform(texUniform2);
    meshstate->addUniform(maximumUniform);
    meshstate->addUniform(maxheightUniform);
    meshstate->addUniform(transparencyUniform);

    // add the geometries to the geode.
    meshGeom->setUseDisplayList(false);
    addDrawable(meshGeom);
}

void Heatmap::setData(float *buffer, float maxheight, float maximum, float transparency)
{
    memcpy(m_data, buffer, m_N*m_K*sizeof(float));

    maximumUniform->set( maximum );
    maxheightUniform->set( maxheight );
    transparencyUniform->set ( transparency );

    m_img2.get()->dirty();
}

Heatmap::~Heatmap()
{
}

liyihongcug

还是有较大问题， osgearth实现应该在球体一个矩形平面绘图 。应该输入数据，还是有差别.
反馈思考还是搞定，借助插件。
解决这个的思路：如果从原生来做，画圆调色等比较慢速度----  参考Baidu地图确实Js都不见得快（heatmap.js）
反复考虑 采用空间数据库来做比较合适i
空间数据库可以更新地理数据和数值大小样式自己定义。
Osgearth本身是支持OGC的，这样只需找打heatmap的插件就可以。完美解决。
唯一动态设置的需要对样式的间隔进行手动设置