Render filled complex polygons with large number of vertices with OpenGL

Question

I use OpenGL to render 2D map and in the process I need to render filled polygons with large number of vertices(100,000+). To do this, I tessellated the polygons to triangles using glu tessellator and rendered the triangles with VBO.

The polygons are rendered successfully. The problem is that the tessellation process turns out to be extremely slow. For some charts with 500,000 vertices, it will take nearly 2 mins on my laptop(i5-3230M 2.6GHz, 8G RAM). This is unacceptable for my application.

Are there any other tessellation algorithm faster than glu tessellator?

Or I have done it wrong?

The following two images are the the rendering results with

glPolygonMode(GL_FRONT, GL_LINE)

EDIT : The map data is static and the original polygon data is in latitude-longitude format. I've already saved the tessellated polygon data (those triangles) in separate file.

To be more clear(Not directly related to the issue),for rendering on screen, a projection is needed to transform the LL format to screen coordinates.

The problem is that the user may have thousands of charts to install(in which the tessellation will be done). Although the tessellation will be run only once, it still takes too long.

Answer 1

A couple of thoughts that I have... can you split the tesselation into organized chunks, ... kind of like a grid? Then if things move you could intelligently only re-tesselate the parts that changed.

Answer 2

is the map static or dynamic?

For static maps

why not store tesselated polygons in some file and not tesselate it again ...

For dynamic maps

would be may be faster using different rendering approach that do not need tesselation to convex polygons like this:

1. clear screen with island color

2. render islands outlines

   not filled primitives like GL_LINE_LOOP do not need to tesselate at all.

3. fill in the watter

   simply start from point outside any polygon and flood fill the map with watter. If the flood fill is coded right (no recursion and fill with lines instead of pixels) then it should take just few [ms]. The problem with this approach is that you need to access the rendered stuff so you need at least 2 passes of rendering. Also implementing flood filling on GPU is not easy.

   There are also alternatives like storing edge points on CPU side and pre-compute the watter filling on CPU side. In that case you need to have list of x coordinates for every y scan line of image that will hold the start and end points for each land. then just fill in the gaps in single render pass...

   This should be rendered in RT easily

[Edit] Grow Fill test Demo

Did some testing with iterative grow filling on your data. There are some problems with the dataset like your polygons overlap which is possibly just holes but as I do not have filling color info but just object ID instead so it is hard to say. Anyway that can be repaired too. Here small win 32 VCL/OpenGL/SW demo with approach I mentioned above (Dynamic maps):

• Win32 Demo use Slow download which is free and no need for any registration just input the code.

It is Win32 stand alone no install using OpenGL+VCL

• mouse wheel zooms
• Shift+mouse wheel selects different polygon
• mouse+Left button pans

There are few issues which can be repaired but as proof of concept it works well. I compiled your ASCII map into binary form (so it loads faster but the form is the same just count of polygons, then count of points per polygon and the points x,y as 64 bit doubles. Counts are 32bit ints)

I am using my own OpenGL engine (which I can not share) so you would need to encode the stuff (like OpenGL,FBO and Texture init/set/usage). Anyway here the C++ code for this VCL app:

//$$---- Form CPP ----
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "win_main.h"
#include "gl/OpenGL3D_double.cpp"
#include "performance.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
// VCL
TMain *Main;
// OpenGL
OpenGLtime tim;
OpenGLscreen scr;
OpenGL_FBO fbo;
GLuint txr_map=-1;
// miscel
int pn=0;                       // vertex count
double px0,px1,py0,py1;         // bbox
double mx,my;                   // mouse
double view[16],iview[16];      // direct and inverse Modelview matrix
double zoom=1.0,dzoom=1.1,viewx=0.0,viewy=0.0;  // view
int index=0;                    // selected polygon
bool _redraw=true;
DWORD cl_water=0xFFEE9040;
DWORD cl_land =0xFF70A0B0;
DWORD cl_edge =0xFF000000;
DWORD cl_sel  =0xFF00FFFF;
AnsiString tcpu,tgpu;
// map
List< List<double> > polygon;   // loaded polygons
      List<double>   water;     // points with water from last frame
//---------------------------------------------------------------------------
void view_compute()
    {
    double x,y;
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glLoadIdentity();
    x=divide(1.0,px1-px0)*scr.aspect;
    y=divide(1.0,py1-py0)*scr._aspect;
    if (x>y) x=y;
    x*=zoom;
    glTranslated(viewx,viewy,0.0);
    glScaled(x,x,1.0);
    glTranslated(-0.5*(px0+px1),-0.5*(py0+py1),0.0);
    glGetDoublev(GL_MODELVIEW_MATRIX,view);
    glPopMatrix();
    matrix_inv(iview,view);
    }
//---------------------------------------------------------------------------
void map_load_csv(AnsiString filename)
    {
    BYTE *dat;
    AnsiString lin,s,s0;
    int ix,i,l,hnd,siz,adr;
    double x,y;
    List< AnsiString > id;

         id.allocate(128);      id.num=0;
    polygon.allocate(128); polygon.num=0;

    hnd=FileOpen(filename,fmOpenRead); if (hnd<0) return;
    siz=FileSeek(hnd,0,2);
        FileSeek(hnd,0,0);
    dat=new BYTE[siz]; if (dat==NULL) { FileClose(hnd); return; }
    siz=FileRead(hnd,dat,siz);
    FileClose(hnd);

    adr=0; txt_load_lin(dat,siz,adr,true);
    for (ix=-1,s0="";adr<siz;)
        {
        lin=txt_load_lin(dat,siz,adr,true);
        if (lin=="") continue;
        i=1; l=lin.Length();
        s=str_load_str(lin,i,true); s=s.SubString(2,s.Length()-2);
        if (s0!=s)
            {
            for (ix=0;ix<id.num;ix++) if (id[ix]==s) break;
            if (ix>=id.num)
                {
                ix=id.num;
                id.add(s);
                polygon.add();
                polygon[ix].allocate(256);
                polygon[ix].num=0;
                }
            s0=s;
            }
        s=str_load_str(lin,i,true); s=s.SubString(2,s.Length()-2); x=str2flt(s);
        s=str_load_str(lin,i,true); s=s.SubString(2,s.Length()-2); y=str2flt(s);
        polygon[ix].add(x);
        polygon[ix].add(y);
        }
    }
//---------------------------------------------------------------------------
void map_save_bin(AnsiString filename)
    {
    int hnd,i;
    hnd=FileCreate(filename); if (hnd<0) return;
    FileWrite(hnd,&polygon.num,4);
    for (i=0;i<polygon.num;i++)
        {
        FileWrite(hnd,&polygon[i].num,4);
        FileWrite(hnd,polygon[i].dat,polygon[i].num*8);
        }
    FileClose(hnd);
    }
//---------------------------------------------------------------------------
void map_load_bin(AnsiString filename)
    {
    int hnd,i,n,m;
    hnd=FileOpen(filename,fmOpenRead); if (hnd<0) return;
    FileRead(hnd,&n,4);
    polygon.allocate(n); polygon.num=n;
    for (i=0;i<n;i++)
        {
        FileRead(hnd,&m,4);
        polygon[i].allocate(m); polygon[i].num=m;
        FileRead(hnd,polygon[i].dat,m*8);
        }
    FileClose(hnd);
    }
//---------------------------------------------------------------------------
void map_bbox()
    {
    int ix,i,n;
    double *p,a;
    pn=0;
    px0=px1=polygon[0][0];
    py0=py1=polygon[0][1];
    for (ix=0;ix<polygon.num;ix++)
        {
        p=polygon[ix].dat;
        n=polygon[ix].num; pn+=n>>1;
        for (i=0;i<n;i+=2)
            {
            a=*p; p++; if (px0>a) px0=a; if (px1<a) px1=a;
            a=*p; p++; if (py0>a) py0=a; if (py1<a) py1=a;
            }
        }
    }
//---------------------------------------------------------------------------
void map_draw()
    {
    int ix,i,n;
    double *p,a;
//  glLineWidth(2.0);
    for (ix=0;ix<polygon.num;ix++)
        {
        p=polygon[ix].dat;
        n=polygon[ix].num;
        if (ix==index) glColor4ubv((BYTE*)&cl_sel);
         else          glColor4ubv((BYTE*)&cl_edge);
        glBegin(GL_LINE_LOOP);
        for (i=0;i<n;i+=2,p+=2) glVertex2dv(p);
        glEnd();
        }
//  glLineWidth(1.0);
    }
//---------------------------------------------------------------------------
void TMain::draw()
    {
    tbeg();
    tim.tbeg();

    // [ render outline to texture ]
    fbo.bind(scr);
    glClearColor(divide((cl_land)&255,255),divide((cl_land>>8)&255,255),divide((cl_land>>16)&255,255),1.0);
    scr.cls();
    glMatrixMode(GL_MODELVIEW);
    glLoadMatrixd(view);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

    if (water.num)  // water start points for grow fill
        {
        // add water around txr border
        glBegin(GL_POINTS);
        glColor4ubv((BYTE*)&cl_water);
        for (int i=0;i<water.num;i+=2)
         glVertex2dv(water.dat+i);
        glEnd();
        }

    map_draw();
    scr.exe();

    fbo.unbind(scr);

    // [ copy GL texture to CPU image ]
    scr.txrs.txr_ld(txr_map);
    // [ create ScanLines for direct pixel access pyx[y][x] ]
    int e,x,y,xs,ys; DWORD **pyx,*p,c0,c1; double a[3];
    xs=scr.txrs.txr.xs;                     // texture resolution (rounded up to power of 2)
    ys=scr.txrs.txr.ys;
    pyx=new DWORD*[ys];
    p=(DWORD*)scr.txrs.txr.txr;             // CPU image pixel data
    for (y=0;y<ys;y++,p+=xs) pyx[y]=p;      // scan line pointers

    // [ Grow Fill water ]
    c0=rgb2bgr(cl_land);
    c1=rgb2bgr(cl_water);
    if (water.num==0)   // first frame view must be set so water is on all borders
        {
        // add water around txr border
        for (x=   1,y=0;y<ys;y++) pyx[y][x]=c1;
        for (x=xs-2,y=0;y<ys;y++) pyx[y][x]=c1;
        for (y=   1,x=0;x<xs;x++) pyx[y][x]=c1;
        for (y=ys-2,x=0;x<xs;x++) pyx[y][x]=c1;
        }

    for (e=1;e;)                            // grow it
    for (e=0,y=1;y<ys-1;y++)
     for (   x=1;x<xs-1;x++)
      if (pyx[y][x]==c0)
       if ((pyx[y-1][x]==c1)
         ||(pyx[y+1][x]==c1)
         ||(pyx[y][x-1]==c1)
         ||(pyx[y][x+1]==c1)) { e=1; pyx[y][x]=c1; }

    // create water start points for next frame
    water.num=0;
    e=4;    // step
    for (y=1;y<ys-2;y+=e)
     for (x=1;x<xs-2;x+=e)
       if ((pyx[y-1][x-1]==c1) // enough water around (x,y)?
         &&(pyx[y-1][x  ]==c1)
         &&(pyx[y-1][x+1]==c1)
         &&(pyx[y  ][x-1]==c1)
         &&(pyx[y  ][x  ]==c1)
         &&(pyx[y  ][x+1]==c1)
         &&(pyx[y+1][x-1]==c1)
         &&(pyx[y+1][x  ]==c1)
         &&(pyx[y+1][x+1]==c1))
            {
            // convert pixel(x,y) -> World(x,y)
            a[0]=divide(2.0*x,xs)-1.0;
            a[1]=divide(2.0*y,ys)-1.0;
            a[2]=0.0;
            matrix_mul_vector(a,iview,a);
            water.add(a[0]);
            water.add(a[1]);
            }

    // [ copy CPU image back to GL texture ]
    delete[] pyx;                           // release ScanLines no need for them anymore
    scr.txrs.txr.rgb2bgr();                 // I got RGB/BGR mismatch somewhere
    scr.txrs.txr_st(txr_map);               // scr.txrs.txr.txr holds pointer to 32bit pixel data
    scr.exe();

    // [ render texture to screen ]
    scr.cls();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    scr.txrs.bind(txr_map);
    glColor3f(1.0,1.0,1.0);
    glBegin(GL_QUADS);
    glTexCoord2f(0.0,0.0); glVertex2f(-1.0,-1.0);
    glTexCoord2f(1.0,0.0); glVertex2f(+1.0,-1.0);
    glTexCoord2f(1.0,1.0); glVertex2f(+1.0,+1.0);
    glTexCoord2f(0.0,1.0); glVertex2f(-1.0,+1.0);
    glEnd();
    scr.txrs.unbind();
    // [info]
    glColor3f(1.0,1.0,1.0);
    scr.text_init_pix(1.0);
    scr.text(tcpu);
    scr.text(tgpu);
    scr.text_exit();

    scr.exe();
    scr.rfs();

    tend(); tcpu=" CPU time: "+tstr(1);
    tim.tend();
    }
//---------------------------------------------------------------------------
void TMain::mouse(double x,double y,TShiftState sh)
    {
    x=divide(2.0*x,scr.xs)-1.0;
    y=1.0-divide(2.0*y,scr.ys);
    if (sh.Contains(ssLeft))
        {
        viewx+=x-mx;
        viewy+=y-my;
        view_compute();
        _redraw=true;
        }
    mx=x;
    my=y;
    }
//---------------------------------------------------------------------------
__fastcall TMain::TMain(TComponent* Owner) : TForm(Owner)
    {
    scr.init(this);
    txr_map=fbo.add(scr);
//  map_load_csv("map.csv");
//  map_save_bin("map.bin");
    map_load_bin("map.bin");
    map_bbox();
    view_compute();
    draw();
    _redraw=true;
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormDestroy(TObject *Sender)
    {
    scr.exit();
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormPaint(TObject *Sender)
    {
    _redraw=true;
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormResize(TObject *Sender)
    {
    scr.resize();
    fbo.resize(scr);
    _redraw=true;
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormMouseWheel(TObject *Sender, TShiftState Shift, int WheelDelta, TPoint &MousePos, bool &Handled)
    {
    if (Shift.Contains(ssShift))
        {
        if (WheelDelta>0) index++; else index--;
        if (index>=polygon.num) index=polygon.num-1;
        if (index<0) index=0;
        _redraw=true;
        }
    else{
        double p[3]={ mx,my,0.0 };
        view_compute();
        matrix_mul_vector(p,iview,p);
        if (WheelDelta>0) zoom*=dzoom; else zoom/=dzoom;
        view_compute();
        matrix_mul_vector(p,view,p);
        viewx-=p[0]-mx;
        viewy-=p[1]-my;
        view_compute();
        _redraw=true;
        }
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormMouseMove(TObject *Sender, TShiftState Shift, int X,int Y) { mouse(X,Y,Shift); }
void __fastcall TMain::FormMouseUp(TObject *Sender, TMouseButton Button, TShiftState Shift, int X, int Y) { mouse(X,Y,Shift); }
void __fastcall TMain::FormMouseDown(TObject *Sender, TMouseButton Button, TShiftState Shift, int X, int Y) { mouse(X,Y,Shift); }
//---------------------------------------------------------------------------
void __fastcall TMain::Timer1Timer(TObject *Sender)
    {
    tgpu=AnsiString().sprintf(" GPU time: [%8.3lf ms]",tim.time());
    if (_redraw) { draw(); _redraw=false; }
    }
//---------------------------------------------------------------------------
void __fastcall TMain::FormDblClick(TObject *Sender)
    {
    Width+=10; // ignore this had some bug in resize FBO texture and this was for debugging it
    }
//---------------------------------------------------------------------------

I also use mine dynamic list template so:

List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items

If you need help with the matrix and vector math routines see this:

• Understanding 4x4 homogenous transform matrices

In the linked answers at the bottom you can find even the C++ implementations I use...

You can ignore the VCL stuff the app just have single Timer on it with interval 40 ms to repaint if needed and fetch measuret GL time if ready...

The important stuff for you is just the draw() routine.

It works like this:

1. bind FBO to enable rendering to texture
2. clear it with land color

3. render polygons outline with edge color

   if you got holes render them with watter color and after filling render them again with edge color

4. render watter start points with watter color

   in the first frame you should have view un-zoomed so all land is surrounded with watter. So the first watter points are the border rectangle of texture.

5. unbind FBO and copy texture pixeldata to CPU side memory

6. grow fill all watter to land color pixels (stop on edge color or any other)

   You can use any fill like Flood fill, segmented line fill etc but beware stack overflow for recursive approach. I decided to use:

   • grow fill

   As it is iterative. It is not as fast (hence the big CPU times but big portion of the CPU time is due to sync while transfering texture between GPU/CPU) but can be speed-up significantly by subdividing image into "square" areas and propagate filling if needed.

7. create watter start points from this image

   so scan whole image (with some step do not need to scan all point) and if found watter add it as watter start point watter for next frame (in world coordinates). This works well until your view will not change too much from frame to frame so limit zoom change and pan step ...

   Here comes an gfx issue when new frame contains watter without any start points and not accessible to other watter. This need some thinking/testing but I think it should be solvable by some static predefined watter start points (few for each polygon) obtained from the first frame.

8. render CPU side image directly or pass it back to GL texture and render it.

Here preview: