Center Of Mass Algorithm

Question

I am implementing a simple center of mass/centroid algorithm for 2D rasters. This is rather trivial on the CPU but has proven difficult to port to the GPU. My CPU version is something like:

float sumX = 0., sumY = 0., c = 0.
for y in imageH:
    for x in imageW:
        if image[x, y] > threshold:
            sumX += x
            sumY += y
            c += 1.
 Point2D centroid = Point2D( sumX / c, sumY / c )

My GPU version (GLSL with backbuffers/Ping Pong) that currently does not work and basically should sum neighbor pixels if they meet the condition (the input texture could have a bitmask image consisting of black and white shapes, for now one shape would suffice and I will start from there).

vec2 pixSize = 1. / resolution;

vec4 originalTexture = texture( myTex, uv );

if( frame == 0 && originalTex.r > treshold )
{

    gl_FragColor = vec4( gl_FragCoord, 0., 0. );

else
{

    for( int y = -1; y <= 1; ++y )
    for( int x = -1; x <= 1; ++x )
    {

        vec4 data = texture( backbuffer, uv + pixSize * vec2( x, y );
        if( originalTexture.r > threshold )
        {

            sumCoord += data.xy;
            counter += 1.;
   
         }

    }

    gl_FragColor = vec4( sumCoord / counter, 0., 0. );

}

EDIT: On second thoughts I see that I am not adding the current fragCoord after the first frame so that is something I must take into account.

Answer 1

1. render ROI 2D texture with unique color for each object

2. for each object (unique color) render pixel at specified location

   Using GLSL shaders compute the CoM for all pixels with the "same" color.

3. glReadPixels at specified locations for each object

4. convert color back to 2D position

Here naive GLSL implementation of mine I just busted for fun:

//---------------------------------------------------------------------------
// Vertex
//---------------------------------------------------------------------------
#version 420 core
//---------------------------------------------------------------------------
layout(location=0) in vec4 vertex;
layout(location=3) in vec4 color;
flat out vec4 col;          // object color
uniform mat4 mvp;       // texture resolution
void main()
    {
    col=color;
    gl_Position=mvp*vertex;
    }
//---------------------------------------------------------------------------

CPU sise C++/OpenGL/VCL code:

//---------------------------------------------------------------------------
// Fragment
//---------------------------------------------------------------------------
#version 420 core
//---------------------------------------------------------------------------
flat in vec4 col;       // object color
out vec4 gl_FragColor;  // fragment output color
uniform sampler2D txr;  // texture to blur
uniform float xs,ys;    // texture resolution
//---------------------------------------------------------------------------
void main()
    {
    float x,y,dx,dy,xx,yy,nn;
    vec3 dc;
    dx=1.0/(xs-1.0); xx=0.0;
    dy=1.0/(ys-1.0); yy=0.0; nn=0.0;
    for (x=0.0;x<=1.0;x+=dx)
     for (y=0.0;y<=1.0;y+=dy)
      if (length(texture(txr,vec2(x,y)).rgb-col.rgb)<0.05)
       { xx+=x; yy+=y; nn++; }
    if (nn!=0){ xx/=nn; yy/=nn; }
    gl_FragColor=vec4(xx,yy,1.0,1.0);
    }
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
GLuint  txr_img;
//---------------------------------------------------------------------------
struct _obj
    {
    DWORD col;  // color
    int x,y;    // CoM
    };
const int _max_objs=32;
_obj obj[_max_objs];
int objs=0;
//---------------------------------------------------------------------------
void gl_draw()
    {
    glClear(GL_COLOR_BUFFER_BIT);
    glDisable(GL_DEPTH_TEST);
    glDisable(GL_CULL_FACE);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(-1.0,-1.0,0.0);
    glScalef(2.0/xs,2.0/ys,1.0);
    glEnable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D,txr_img);

    // compute CoM using GLSL
    float m[16];
    int i,x,y,r;
    glUseProgram(prog_id);
    glUniform1i(glGetUniformLocation(prog_id,"txr" ),0);
    glUniform1f(glGetUniformLocation(prog_id,"xs"  ),xs);
    glUniform1f(glGetUniformLocation(prog_id,"ys"  ),ys);
    glGetFloatv(GL_MODELVIEW_MATRIX,m);
    glUniformMatrix4fv(glGetUniformLocation(prog_id,"mvp" ),1,false,m);
    glBegin(GL_POINTS);
    for (x=0,y=0,i=0;i<objs;i++,x++)
        {
        glColor4ubv((BYTE*)&(obj[i].col));
        glVertex2i(x,y);
        }
    glEnd();
    glUseProgram(0);
    glFlush();

    // obtain data
    DWORD dat[_max_objs];
    glReadPixels(0,0,_max_objs,1,GL_BGRA,GL_UNSIGNED_BYTE,dat);
    for (i=0;i<objs;i++)
        {
        obj[i].x=(xs*(    ((dat[i]>>16)&255)))>>8;
        obj[i].y=(ys*(255-((dat[i]>> 8)&255)))>>8;
        }

    // render image (old API)
    glBegin(GL_QUADS);
    glColor3f(1.0,1.0,1.0);
    glTexCoord2f(0.0,1.0); glVertex2i( 0, 0);
    glTexCoord2f(0.0,0.0); glVertex2i( 0,ys);
    glTexCoord2f(1.0,0.0); glVertex2i(xs,ys);
    glTexCoord2f(1.0,1.0); glVertex2i(xs, 0);
    glEnd();

    glBindTexture(GL_TEXTURE_2D,0);
    glDisable(GL_TEXTURE_2D);

    // render computed CoM
    r=10;
    glLineWidth(5.0);
    glColor3f(0.1,0.1,0.1);
    glBegin(GL_LINES);
    for (i=0;i<objs;i++)
        {
        x=obj[i].x;
        y=obj[i].y;
        glVertex2i(x-r,y-r);
        glVertex2i(x+r,y+r);
        glVertex2i(x-r,y+r);
        glVertex2i(x+r,y-r);
        }
    glEnd();
    r=9;
    glLineWidth(2.0);
    glBegin(GL_LINES);
    for (i=0;i<objs;i++)
        {
        x=obj[i].x;
        y=obj[i].y;
        glColor4ubv((BYTE*)&obj[i].col);
        glVertex2i(x-r,y-r);
        glVertex2i(x+r,y+r);
        glVertex2i(x-r,y+r);
        glVertex2i(x+r,y-r);
        }
    glEnd();

    glFlush();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    gl_init(Handle);

    // textures
    Byte q;
    unsigned int *pp;
    int i,xs,ys,x,y,adr,*txr;
    union { unsigned int c32; Byte db[4]; } c;
    Graphics::TBitmap *bmp=new Graphics::TBitmap;   // new bmp

    objs=0;

    // image texture
    bmp->LoadFromFile("texture.bmp");   // load from file
    bmp->HandleType=bmDIB;      // allow direct access to pixels
    bmp->PixelFormat=pf32bit;   // set pixel to 32bit so int is the same size as pixel
    xs=bmp->Width;              // resolution should be power of 2
    ys=bmp->Height;
    txr=new int[xs*ys];         // create linear framebuffer
    for(adr=0,y=0;y<ys;y++)
        {
        pp=(unsigned int*)bmp->ScanLine[y];
        for(x=0;x<xs;x++,adr++)
            {
            // rgb2bgr and copy bmp -> txr[]
            c.c32=pp[x]; c.db[3]=0;
            q      =c.db[2];
            c.db[2]=c.db[0];
            c.db[0]=q;
            txr[adr]=c.c32;
            for (i=0;i<objs;i++)
             if (obj[i].col==c.c32)
              { i=-1; break; }
            if (i>=0)
                {
                obj[objs].col=c.c32;
                obj[objs].x=0;
                obj[objs].y=0;
                objs++;
                }
            }
        }
    glGenTextures(1,&txr_img);
    glEnable(GL_TEXTURE_2D);    // copy it to gfx card
    glBindTexture(GL_TEXTURE_2D,txr_img);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 4);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S,GL_CLAMP);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T,GL_CLAMP);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER,GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER,GL_NEAREST);
    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,GL_MODULATE);
//  glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE,GL_DECAL);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, xs, ys, 0, GL_RGBA, GL_UNSIGNED_BYTE, txr);
    glDisable(GL_TEXTURE_2D);
    delete[] txr;
    delete bmp;

    int hnd,siz; char vertex[4096],fragment[4096];
    hnd=FileOpen("center_of_mass.glsl_vert",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,vertex  ,siz); vertex  [siz]=0; FileClose(hnd);
    hnd=FileOpen("center_of_mass.glsl_frag",fmOpenRead); siz=FileSeek(hnd,0,2); FileSeek(hnd,0,0); FileRead(hnd,fragment,siz); fragment[siz]=0; FileClose(hnd);
    glsl_init(vertex,NULL,fragment); // old version of gl_simple is without the NULL !!!
    hnd=FileCreate("GLSL.txt"); FileWrite(hnd,glsl_log,glsl_logs); FileClose(hnd);

    ClientWidth=xs;
    ClientHeight=ys;
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    glDeleteTextures(1,&txr_img);
    gl_exit();
    glsl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    gl_resize(ClientWidth,ClientHeight);
    gl_draw();
    }
//---------------------------------------------------------------------------

it uses gl_simple.h

This is input ROI texture:

And output:

for better precision you can use more pixels per single object output, or use framebuffer with more than 8bit per channel ...

Answer 2

This is a straight-forward reduction problem.

Create a three-channel texture T with the values

if(image[x,y] > threshold) {
    T[x,y,0] = x;
    T[x,y,1] = y;
    T[x,y,2] = 1;
} else {
    T[x,y,0] = T[x,y,1] = T[x,y,2] = 0;
}

Then reduce the texture by a constant factor till you get a single 1x1 pixel. This takes O(log(max(w,h))) GPU friendly steps, thus is extremely efficient. Here are some slides about how to do reduction on GPU. In fact even something like glGenerateMipmap() will work, though the exact filter it uses is implementation-specefic.

Next, after you reduced the texture to a single pixel, you extract the center of mass by

x = T[0,0,0]/T[0,0,2];
y = T[0,0,1]/T[0,0,2];

Done.

EDIT: a simple factor-two reduction with a fragment shader can be done as simple as follows:

layout(binding = 0) uniform sampler2D TEX;
void main() {
    vec2 uv = 2*gl_FragCoord.xy/textureSize(TEX, 0);
    gl_FragColor = texture(TEX, uv);
}

You need to bind the previous texture to TEX with a bilinear minification filter and (0,0,0,0) clamp-to-border wrapping mode. Setup the next texture layer as the render target.