is it possible to speed-up matlab plotting by calling c / c++ code in matlab?

Question

It is generally very easy to call mex files (written in c/c++) in Matlab to speed up certain calculations. In my experience however, the true bottleneck in Matlab is data plotting. Creating handles is extremely expensive and even if you only update handle data (e.g., XData, YData, ZData), this might take ages. Even worse, since Matlab is a single threaded program, it is impossible to update multiple plots at the same time.

Therefore my question: Is it possible to write a Matlab GUI and call C++ (or some other parallelizable code) which would take care of the plotting / visualization? I'm looking for a cross-platform solution that will work on Windows, Mac and Linux, but any solution that get's me started on either OS is greatly appreciated!

I found a C++ library that seems to use Matlab's plot() syntax but I'm not sure whether this would speed things up, since I'm afraid that if I plot into Matlab's figure() window, things might get slowed down again.

I would appreciate any comments and feedback from people who have dealt with this kind of situation before!

EDIT: obviously, I've already profiled my code and the bottleneck is the plotting (dozen of panels with lots of data).

EDIT2: for you to get the bounty, I need a real life, minimal working example on how to do this - suggestive answers won't help me.

EDIT3: regarding the data to plot: in a most simplistic case, think about 20 line plots, that need to be updated each second with something like 1000000 data points.

EDIT4: I know that this is a huge amount of points to plot but I never said that the problem was easy. I can not just leave out certain data points, because there's no way of assessing what points are important, before actually plotting them (data is sampled a sub-ms time resolution). As a matter of fact, my data is acquired using a commercial data acquisition system which comes with a data viewer (written in c++). This program has no problem visualizing up to 60 line plots with even more than 1000000 data points.

EDIT5: I don't like where the current discussion is going. I'm aware that sub-sampling my data might speeds up things - however, this is not the question. The question here is how to get a c / c++ / python / java interface to work with matlab in order hopefully speed up plotting by talking directly to the hardware (or using any other trick / way)

Answer 1

Did you try the trivial solution of changing the render method to OpenGL ?

opengl hardware;
set(gcf,'Renderer','OpenGL');

Warning! There will be some things that disappear in this mode, and it will look a bit different, but generally plots will runs much faster, especially if you have a hardware accelerator.

By the way, are you sure that you will actually gain a performance increase? For example, in my experience, WPF graphics in C# are considerably slower than Matlabs, especially scatter plot and circles.

Edit: I thought about the fact that the number of points that is actually drawn to the screen can't be that much. Basically it means that you need to interpolate at the places where there is a pixel in the screen. Check out this object:

classdef InterpolatedPlot < handle
    properties(Access=private)
        hPlot;
    end

    methods(Access=public)
        function this = InterpolatedPlot(x,y,varargin)
            this.hPlot = plot(0,0,varargin{:});
            this.setXY(x,y);
        end        
    end    

    methods
        function setXY(this,x,y)
            parent = get(this.hPlot,'Parent');
            set(parent,'Units','Pixels')
            sz = get(parent,'Position');
            width = sz(3); %Actual width in pixels
            subSampleX = linspace(min(x(:)),max(x(:)),width);

            subSampleY = interp1(x,y,subSampleX);
            set(this.hPlot,'XData',subSampleX,'YData',subSampleY);
        end

    end
end

And here is an example how to use it:

function TestALotOfPoints()
    x = rand(10000,1);
    y = rand(10000,1);

    ip = InterpolatedPlot(x,y,'color','r','LineWidth',2);

end

---

Another possible improvement: Also, if your x data is sorted, you can use interp1q instead of interp, which will be much faster.

classdef InterpolatedPlot < handle
    properties(Access=private)
        hPlot;
    end

%     properties(Access=public)
%         XData;
%         YData;      
%     end

    methods(Access=public)
        function this = InterpolatedPlot(x,y,varargin)
            this.hPlot = plot(0,0,varargin{:});
            this.setXY(x,y);
%             this.XData = x;
%             this.YData = y;
        end        
    end    

    methods
        function setXY(this,x,y)
            parent = get(this.hPlot,'Parent');
            set(parent,'Units','Pixels')
            sz = get(parent,'Position');
            width = sz(3); %Actual width in pixels
            subSampleX = linspace(min(x(:)),max(x(:)),width);

            subSampleY = interp1q(x,y,transpose(subSampleX));
            set(this.hPlot,'XData',subSampleX,'YData',subSampleY);
        end

    end
end

And the use case:

function TestALotOfPoints()
    x = rand(10000,1);
    y = rand(10000,1);
    x = sort(x);
    ip = InterpolatedPlot(x,y,'color','r','LineWidth',2);

end

Answer 2

Since you want maximum performance you should consider writing a minimal OpenGL viewer. Dump all the points to a file and launch the viewer using the "system"-command in MATLAB. The viewer can be really simple. Here is one implemented using GLUT, compiled for Mac OS X. The code is cross platform so you should be able to compile it for all the platforms you mention. It should be easy to tweak this viewer for your needs.

If you are able to integrate this viewer more closely with MATLAB you might be able to get away with not having to write to and read from a file (= much faster updates). However, I'm not experienced in the matter. Perhaps you can put this code in a mex-file?

EDIT: I've updated the code to draw a line strip from a CPU memory pointer.

// On Mac OS X, compile using: g++ -O3 -framework GLUT -framework OpenGL glview.cpp
// The file "input" is assumed to contain a line for each point:
// 0.1 1.0
// 5.2 3.0
#include <vector>
#include <sstream>
#include <fstream>
#include <iostream>
#include <GLUT/glut.h>
using namespace std;
struct float2 { float2() {} float2(float x, float y) : x(x), y(y) {} float x, y; };
static vector<float2> points;
static float2 minPoint, maxPoint;
typedef vector<float2>::iterator point_iter;
static void render() {
    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
    glClear(GL_COLOR_BUFFER_BIT);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(minPoint.x, maxPoint.x, minPoint.y, maxPoint.y, -1.0f, 1.0f);
    glColor3f(0.0f, 0.0f, 0.0f);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(2, GL_FLOAT, sizeof(points[0]), &points[0].x);
    glDrawArrays(GL_LINE_STRIP, 0, points.size());
    glDisableClientState(GL_VERTEX_ARRAY);
    glutSwapBuffers();
}
int main(int argc, char* argv[]) {
    ifstream file("input");
    string line;
    while (getline(file, line)) {
        istringstream ss(line);
        float2 p;
        ss >> p.x;
        ss >> p.y;
        if (ss)
            points.push_back(p);
    }
    if (!points.size())
        return 1;
    minPoint = maxPoint = points[0];
    for (point_iter i = points.begin(); i != points.end(); ++i) {
        float2 p = *i;
        minPoint = float2(minPoint.x < p.x ? minPoint.x : p.x, minPoint.y < p.y ? minPoint.y : p.y);
        maxPoint = float2(maxPoint.x > p.x ? maxPoint.x : p.x, maxPoint.y > p.y ? maxPoint.y : p.y);
    }
    float dx = maxPoint.x - minPoint.x;
    float dy = maxPoint.y - minPoint.y;
    maxPoint.x += dx*0.1f; minPoint.x -= dx*0.1f;
    maxPoint.y += dy*0.1f; minPoint.y -= dy*0.1f;
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
    glutInitWindowSize(512, 512);
    glutCreateWindow("glview");
    glutDisplayFunc(render);
    glutMainLoop();
    return 0;
}

EDIT: Here is new code based on the discussion below. It renders a sin function consisting of 20 vbos, each containing 100k points. 10k new points are added each rendered frame. This makes a total of 2M points. The performance is real-time on my laptop.

// On Mac OS X, compile using: g++ -O3 -framework GLUT -framework OpenGL glview.cpp
#include <vector>
#include <sstream>
#include <fstream>
#include <iostream>
#include <cmath>
#include <iostream>
#include <GLUT/glut.h>
using namespace std;
struct float2 { float2() {} float2(float x, float y) : x(x), y(y) {} float x, y; };
struct Vbo {
    GLuint i;
    Vbo(int size) { glGenBuffersARB(1, &i); glBindBufferARB(GL_ARRAY_BUFFER, i); glBufferDataARB(GL_ARRAY_BUFFER, size, 0, GL_DYNAMIC_DRAW); } // could try GL_STATIC_DRAW
    void set(const void* data, size_t size, size_t offset) { glBindBufferARB(GL_ARRAY_BUFFER, i); glBufferSubData(GL_ARRAY_BUFFER, offset, size, data); }
    ~Vbo() { glDeleteBuffers(1, &i); }
};
static const int vboCount = 20;
static const int vboSize = 100000;
static const int pointCount = vboCount*vboSize;
static float endTime = 0.0f;
static const float deltaTime = 1e-3f;
static std::vector<Vbo*> vbos;
static int vboStart = 0;
static void addPoints(float2* points, int pointCount) {
    while (pointCount) {
        if (vboStart == vboSize || vbos.empty()) {
            if (vbos.size() >= vboCount+2) { // remove and reuse vbo
                Vbo* first = *vbos.begin();
                vbos.erase(vbos.begin());
                vbos.push_back(first);
            }
            else { // create new vbo
                vbos.push_back(new Vbo(sizeof(float2)*vboSize));
            }
            vboStart = 0;
        }

        int pointsAdded = pointCount;

        if (pointsAdded + vboStart > vboSize)
            pointsAdded = vboSize - vboStart;

        Vbo* vbo = *vbos.rbegin();
        vbo->set(points, pointsAdded*sizeof(float2), vboStart*sizeof(float2));

        pointCount -= pointsAdded;
        points += pointsAdded;
        vboStart += pointsAdded;
    }
}
static void render() {
    // generate and add 10000 points
    const int count = 10000;
    float2 points[count];
    for (int i = 0; i < count; ++i) {
        float2 p(endTime, std::sin(endTime*1e-2f));
        endTime += deltaTime;
        points[i] = p;
    }
    addPoints(points, count);
    // render
    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
    glClear(GL_COLOR_BUFFER_BIT);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(endTime-deltaTime*pointCount, endTime, -1.0f, 1.0f, -1.0f, 1.0f);
    glColor3f(0.0f, 0.0f, 0.0f);
    glEnableClientState(GL_VERTEX_ARRAY);
    for (size_t i = 0; i < vbos.size(); ++i) {
        glBindBufferARB(GL_ARRAY_BUFFER, vbos[i]->i);
        glVertexPointer(2, GL_FLOAT, sizeof(float2), 0);

        if (i == vbos.size()-1)
            glDrawArrays(GL_LINE_STRIP, 0, vboStart);
        else
            glDrawArrays(GL_LINE_STRIP, 0, vboSize);
    }
    glDisableClientState(GL_VERTEX_ARRAY);
    glutSwapBuffers();
    glutPostRedisplay();
}
int main(int argc, char* argv[]) {
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
    glutInitWindowSize(512, 512);
    glutCreateWindow("glview");
    glutDisplayFunc(render);
    glutMainLoop();
    return 0;
}