3
\$\begingroup\$

The authors of this paper compute support_points of a 900×700 image in 118 ms. I have implemented their algorithm below in Halide. In my algorithm, the nested for loops over length and width iterate over xi and yi, which are points in output_x and output_y. Over each iteration of the nested for loops, a vector top_k is computed and pushed_back into support_points. Computing this pipeline even for left_buffer.width() == 20 and left_buffer.height() == 20 takes 500 ms. Thus this implementation is several orders of magnitude slower.

// To compile and run:
// export LD_LIBRARY_PATH=../bin
// g++ support.cpp -g -I ../include -I ../tools -L ../bin -lHalide `libpng-config --cflags --ldflags` -ljpeg -lpthread -ldl -o support -std=c++11 && ./support



#include<stdio.h>
#include<string>
#include "Halide.h"
#include "halide_image_io.h"
#include <math.h>
#include "clock.h"
#include <vector>

using namespace std;
using namespace Halide::Tools;
using namespace Halide::ConciseCasts;
using namespace Halide;

int main(int argc, char **argv) {
    double t1 = current_time();

    Var x("x"), y("y");

    Buffer<uint8_t> left_buffer = load_image("images/stereo/bike_smallest.jpg");

    Expr clamped_x = clamp(x, 0, left_buffer.width() - 1);
    Expr clamped_y = clamp(y, 0, left_buffer.height() - 1);

    Func left_original("left_original");
    left_original(x, y) = f32(left_buffer(clamped_x, clamped_y));
    left_original.compute_root();

    // 3x3 sobel filter
    Buffer<float_t> sobel_1(3);
    sobel_1(0) = -1;
    sobel_1(1) = 0;
    sobel_1(2) = 1;

    Buffer<float_t> sobel_2(3);
    sobel_2(0) = 1;
    sobel_2(1) = 2;
    sobel_2(2) = 1;

    Func output_x_inter("output_x_inter");
    output_x_inter(x, y) = f32(left_original(x + 1, y) * sobel_1(0) + left_original(x, y) * sobel_1(1) + left_original(x - 1, y) * sobel_1(2));
    output_x_inter.compute_root();

    Func sobel_x("sobel_x");
    sobel_x(x, y) = f32(output_x_inter(x, y + 1) * sobel_2(0) + output_x_inter(x,  y) * sobel_2(1) + output_x_inter(x, y - 1) * sobel_2(2));
//    sobel_x.compute_root(); didn't work
    sobel_x.trace_stores();
    Buffer<float> output_x = sobel_x.realize(left_buffer.width(), left_buffer.height());
    output_x.set_name("output_x");

    Func sobel_y_inter("sobel_y_inter");
    sobel_y_inter(x, y) = f32(left_original(x + 1, y) * sobel_2(0) + left_original(x, y) * sobel_2(1) + left_original(x - 1, y) * sobel_2(2));

    Func sobel_y("sobel_y");
    sobel_y(x, y) = f32(sobel_y_inter(x, y + 1) * sobel_1(0) + sobel_y_inter(x, y) * sobel_1(1) + sobel_y_inter(x, y - 1) * sobel_1(2));
    sobel_y.compute_root();
    sobel_y.trace_stores();
    Buffer<float> output_y = sobel_y.realize(left_buffer.width(), left_buffer.height());
    output_y.set_name("output_y");


    int k = 4; // # of support points
    vector<pair<Expr, Expr>> support_points(k * left_buffer.width() * left_buffer.height());
    // Calculate support pixel for each
    Func support("support");
    support(x, y) = Tuple(i32(0), i32(0), f32(0));



    for (int yi = 0; yi < left_buffer.height(); yi++) {
        for (int xi = 0; xi < left_buffer.width() - 2; xi++) {
            bool left = xi < left_buffer.width() / 4;
            bool center = (xi >= left_buffer.width() / 4 && xi < left_buffer.width() * 3 / 4);
            bool right = xi >= left_buffer.width() * 3 / 4;

            vector <pair<Expr, Expr>> scan_range;
            pair <Expr, Expr> scan_height(0, (Expr) left_buffer.height());
            pair <Expr, Expr> scan_width;
            int which_pred = 0;
            if (left) {

                    scan_width = make_pair((Expr) 0, (Expr) left_buffer.width() / 2);
                    which_pred = 0;
            }
            else if (center) {
                    scan_width = make_pair((Expr) xi - left_buffer.width() / 4, (Expr) left_buffer.width() / 2);
                    which_pred = 1;
            }
            else if (right) {
                    scan_width = make_pair((Expr) left_buffer.width() / 2, (Expr) left_buffer.width() / 2);
                    which_pred = 2;
            }
            else {
                cout<<"Error"<<endl;
            }

            scan_range = {scan_width, scan_height};
//            cout<<"xi "<<xi<<endl;
//            cout<<"yi "<<yi<<endl;
//            cout<<"scan_width= "<<scan_width.first<<" "<<scan_width.second<<endl;
//            cout<<"scan_height= "<<scan_height.first<<" "<<scan_height.second<<endl;


            RDom scanner(scan_range);
            Expr predicate[3] = {scanner.x != xi && scanner.y != yi, scanner.x != 0 && scanner.y != 0, scanner.x != xi && scanner.y != yi};
            scanner.where(predicate[which_pred]);
            std::vector<Expr> top_k(k * 3);
            for (int i = 0; i < k; i++) { // say we want top 4 support points.
                top_k[3*i] = 10000.0f;
                top_k[3*i+1] = 0;
                top_k[3*i+2] = 0;
            }

            Func argmin("argmin");
            argmin() = Tuple(top_k);
            Expr next_val = abs(output_x(xi, yi) - output_x(scanner.x, scanner.y)) + abs(output_y(xi, yi) - output_y(scanner.x, scanner.y));
            Expr next_x = scanner.x;
            Expr next_y = scanner.y;

            top_k = Tuple(argmin()).as_vector();
            // Insert a single element into a sorted list without actually branching
            top_k.push_back(next_val);
            top_k.push_back(next_x);
            top_k.push_back(next_y);
            for (int i = k; i > 0; i--) {
                Expr prev_val = top_k[(i-1)*3];
                Expr prev_x = top_k[(i-1)*3 + 1];
                Expr prev_y = top_k[(i-1)*3 + 2];
                Expr should_swap = top_k[i*3] < prev_val;

                top_k[(i-1)*3] = select(should_swap, top_k[i*3], prev_val);
                top_k[(i-1)*3 + 1] = select(should_swap, top_k[i*3 + 1], prev_x);
                top_k[(i-1)*3 + 2] = select(should_swap, top_k[i*3 + 2], prev_y);
                top_k[i*3] = select(should_swap, prev_val, top_k[i*3]);
                top_k[i*3 + 1] = select(should_swap, prev_x, top_k[i*3 + 1]);
                top_k[i*3 + 2] = select(should_swap, prev_y, top_k[i*3 + 2]);
            }
            // Discard the k+1th element
            top_k.pop_back(); top_k.pop_back(); top_k.pop_back();

            bool cond = xi == 10 && yi == 10;
            cout << xi << " "<< yi << " " << cond << endl;

            Expr e = argmin()[0];

            e = print_when(cond, e, "<- argmin() val");
            argmin() = Tuple(top_k);
            argmin.compute_root();
//            argmin.trace_stores();


            argmin.compile_to_lowered_stmt("argmin.html", {}, HTML);
            Realization real = argmin.realize();
            for (int i = 0; i < k; i++) {
                pair<Expr, Expr> c(top_k[3*i+1], top_k[3*i+2]);
                support_points.push_back(c);
            }
        }
    }
    double t2 = current_time();

    cout<<(t2-t1)/100<<" ms"<<endl;
    cout<<"executed"<<endl;
}

How can I make it more efficient?

----- EDIT -----

As an answer suggested, I pulled out all the declarations out of the iterations, but this has had only a marginal effect, if any, on the performance. This is what my code looks like now:

#include<stdio.h>
#include<string>
#include "Halide.h"
#include "halide_image_io.h"
#include <math.h>
#include "clock.h"
#include <vector>

#define R 0
#define G 1
#define B 2
using namespace std;
using namespace Halide::Tools;
using namespace Halide::ConciseCasts;
using namespace Halide;

int main(int argc, char **argv) {
    double t1 = current_time();

    Var x("x"), y("y");

    Buffer<uint8_t> left_buffer = load_image("images/stereo/bike_smallest.jpg");

    Expr clamped_x = clamp(x, 0, left_buffer.width() - 1);
    Expr clamped_y = clamp(y, 0, left_buffer.height() - 1);

    Func left_original("left_original");
    left_original(x, y) = f32(left_buffer(clamped_x, clamped_y));
    left_original.compute_root();

    // 3x3 sobel filter
    Buffer<float_t> sobel_1(3);
    sobel_1(0) = -1;
    sobel_1(1) = 0;
    sobel_1(2) = 1;

    Buffer<float_t> sobel_2(3);
    sobel_2(0) = 1;
    sobel_2(1) = 2;
    sobel_2(2) = 1;

    Func output_x_inter("output_x_inter");
    output_x_inter(x, y) = f32(left_original(x + 1, y) * sobel_1(0) + left_original(x, y) * sobel_1(1) + left_original(x - 1, y) * sobel_1(2));
    output_x_inter.compute_root();

    Func sobel_x("sobel_x");
    sobel_x(x, y) = f32(output_x_inter(x, y + 1) * sobel_2(0) + output_x_inter(x,  y) * sobel_2(1) + output_x_inter(x, y - 1) * sobel_2(2));
//    sobel_x.compute_root(); didn't work
//    sobel_x.trace_stores();
    Buffer<float> output_x = sobel_x.realize(left_buffer.width(), left_buffer.height());
    output_x.set_name("output_x");

    Func sobel_y_inter("sobel_y_inter");
    sobel_y_inter(x, y) = f32(left_original(x + 1, y) * sobel_2(0) + left_original(x, y) * sobel_2(1) + left_original(x - 1, y) * sobel_2(2));

    Func sobel_y("sobel_y");
    sobel_y(x, y) = f32(sobel_y_inter(x, y + 1) * sobel_1(0) + sobel_y_inter(x, y) * sobel_1(1) + sobel_y_inter(x, y - 1) * sobel_1(2));
    sobel_y.compute_root();
//    sobel_y.trace_stores();
    Buffer<float> output_y = sobel_y.realize(left_buffer.width(), left_buffer.height());
    output_y.set_name("output_y");


    int k = 4; // # of support points
    vector<pair<Expr, Expr>> support_points(k * left_buffer.width() * left_buffer.height());
    // Calculate support pixel for each
    Func support("support");
    support(x, y) = Tuple(i32(0), i32(0), f32(0));


    vector <pair<Expr, Expr>> scan_range(2);
    pair <Expr, Expr> scan_height(0, (Expr) left_buffer.height());
    pair <Expr, Expr> scan_width;

    std::vector<Expr> top_k(k * 3);
    Expr predicate[3];
    int count = 0;
    RDom scanner;
    bool left;
    bool center;
    bool right;
    int which_pred = 0;
    Expr next_val, next_x, next_y, prev_val, prev_x, prev_y;
    Expr should_swap;
    pair<Expr, Expr> c;


    for (int yi = 0; yi < left_buffer.height(); yi++) {
        for (int xi = 0; xi < left_buffer.width() - 2; xi++) {

            left = xi < left_buffer.width() / 4;
            center = (xi >= left_buffer.width() / 4 && xi < left_buffer.width() * 3 / 4);
            right = xi >= left_buffer.width() * 3 / 4;

            if (left) {

                    scan_width = make_pair((Expr) 0, (Expr) left_buffer.width() / 2);
                    which_pred = 0;
            }
            else if (center) {
                    scan_width = make_pair((Expr) xi - left_buffer.width() / 4, (Expr) left_buffer.width() / 2);
                    which_pred = 1;
            }
            else if (right) {
                    scan_width = make_pair((Expr) left_buffer.width() / 2, (Expr) left_buffer.width() / 2);
                    which_pred = 2;
            }
            else {
                cout<<"Error"<<endl;
            }

            scan_range[0] = scan_width;
            scan_range[1] = scan_height;

            // so far 1e-05 per iteration

            scanner = (scan_range); //slow 0.006

            // these three kind of slow 0.002

            predicate[0] = scanner.x != xi && scanner.y != yi;
            predicate[1] = scanner.x != 0 && scanner.y != 0;
            predicate[2] = scanner.x != xi && scanner.y != yi;
            scanner.where(predicate[which_pred]);


            // this loop is 1e-05
            for (int i = 0; i < k; i++) { // say we want top 4 support points.
                top_k[3*i] = 10000.0f;
                top_k[3*i+1] = 0;
                top_k[3*i+2] = 0;
            }

            Func argmin("argmin");

            // 0.003
            argmin() = Tuple(top_k);


            next_val = abs(output_x(xi, yi) - output_x(scanner.x, scanner.y)) + abs(output_y(xi, yi) - output_y(scanner.x, scanner.y));

            next_x = scanner.x;
            next_y = scanner.y;

            top_k = Tuple(argmin()).as_vector();

            // Insert a single element into a sorted list without actually branching
            top_k.push_back(next_val);
            top_k.push_back(next_x);
            top_k.push_back(next_y);

            for (int i = k; i > 0; i--) {
                prev_val = top_k[(i-1)*3];
                prev_x = top_k[(i-1)*3 + 1];
                prev_y = top_k[(i-1)*3 + 2];
                should_swap = top_k[i*3] < prev_val;

                top_k[(i-1)*3] = select(should_swap, top_k[i*3], prev_val);
                top_k[(i-1)*3 + 1] = select(should_swap, top_k[i*3 + 1], prev_x);
                top_k[(i-1)*3 + 2] = select(should_swap, top_k[i*3 + 2], prev_y);
                top_k[i*3] = select(should_swap, prev_val, top_k[i*3]);
                top_k[i*3 + 1] = select(should_swap, prev_x, top_k[i*3 + 1]);
                top_k[i*3 + 2] = select(should_swap, prev_y, top_k[i*3 + 2]);
//                cout<<"(165)"<<endl;
            }

            // Discard the k+1th element
            top_k.pop_back(); top_k.pop_back(); top_k.pop_back();

            argmin() = Tuple(top_k);
            argmin.compute_root();
//            argmin.trace_stores();

            argmin.compile_to_lowered_stmt("argmin.html", {}, HTML);
//            cout<<"(184)"<<endl;
            Realization real = argmin.realize(); //1.5 ms
//            cout<<"(186)"<<endl;

            for (int i = 0; i < k; i++) {
                c.first = top_k[3*i+1];
                c.second = top_k[3*i+2];
                support_points.push_back(c);
            }

            count++;
        }
    }
    double t2 = current_time();
    cout<<count<<" "<<(t2-t1)/100<<" ms"<<endl;
    cout<<"executed"<<endl;
}

```
\$\endgroup\$
0
\$\begingroup\$

This part looks ripe for refactoring:

        top_k = Tuple(argmin()).as_vector();
        // Insert a single element into a sorted list without actually branching
        top_k.push_back(next_val);
        top_k.push_back(next_x);
        top_k.push_back(next_y);
        for (int i = k; i > 0; i--) {
            Expr prev_val = top_k[(i-1)*3];
            Expr prev_x = top_k[(i-1)*3 + 1];
            Expr prev_y = top_k[(i-1)*3 + 2];
            Expr should_swap = top_k[i*3] < prev_val;

            top_k[(i-1)*3] = select(should_swap, top_k[i*3], prev_val);
            top_k[(i-1)*3 + 1] = select(should_swap, top_k[i*3 + 1], prev_x);
            top_k[(i-1)*3 + 2] = select(should_swap, top_k[i*3 + 2], prev_y);
            top_k[i*3] = select(should_swap, prev_val, top_k[i*3]);
            top_k[i*3 + 1] = select(should_swap, prev_x, top_k[i*3 + 1]);
            top_k[i*3 + 2] = select(should_swap, prev_y, top_k[i*3 + 2]);
        }
        // Discard the k+1th element
        top_k.pop_back(); top_k.pop_back(); top_k.pop_back();

If I understand correctly, what you're doing here is bubble-sorting a triple (next_val, next_x, next_y) into the proper place in vector top_k, and then removing the now-smallest triple from the vector.

The first thing you could do is use proper data types. Instead of triples of ints (and the continual multiplications-by-3 that entails), you should make a struct type, like

struct Triple {
    int x;
    int y;
    int value;
    struct descending {
        bool operator()(const Triple& a, const Triple& b) const {
            return a.value > b.value;
        }
    };
};

And then your insertion can become

std::vector<Triple> top_k = [...]
Triple new_item = { next_x, next_y, next_val };
auto insertion_point = std::lower_bound(top_k.begin(), top_k.end(), new_item, Triple::descending);
top_k.insert(insertion_point, new_item);
top_k.pop_back();

(Or maybe my descending should be ascending. I didn't expend many brain cells on it.)


If your algorithm depends on this operation, then (A) maybe you should use a priority_queue, (B) maybe you should use a priority_queue that supports the replace_top operation, and (C) certainly you should use something that never allocates memory in your inner loop.

Right now it seems you're spending a lot of time resizing STL containers in the inner loop, which causes heap traffic. For performance, instead of

for (int yi = 0; yi < left_buffer.height(); yi++) {
    for (int xi = 0; xi < left_buffer.width() - 2; xi++) {
        [...]
        vector<pair<Expr, Expr>> scan_range = {scan_width, scan_height};

you should do

vector<pair<Expr, Expr>> scan_range(2);
for (int yi = 0; yi < left_buffer.height(); yi++) {
    for (int xi = 0; xi < left_buffer.width() - 2; xi++) {
        [...]
        scan_range[0] = scan_width;
        scan_range[1] = scan_height;

This puts the heap allocation outside the inner loop, which means now you're doing 1 allocation where in your original code you're doing left_buffer.height()*(left_buffer.width() - 2) allocations.

Applying this strategy to your other allocations (e.g. top_k itself) should give a huge performance boost.


EDIT: Looking for more wasted cycles in your updated code, I see the following four lines at the end of your loop:

        argmin.compile_to_lowered_stmt("argmin.html", {}, HTML);
//            cout<<"(184)"<<endl;
        Realization real = argmin.realize(); //1.5 ms
//            cout<<"(186)"<<endl;

Variable real is never used, so you can eliminate it. If I'm interpreting that comment correctly, eliminating real will save you 1.5ms per loop.

argmin.compile_to_lowered_stmt(...) seems to be Halide's equivalent of a "debugging printf." It will do a lot of work to compile the function, and then also open a file, write to it, and close the file again. Those are both horribly slow operations. So unless you need argmin.html — and I don't see how you possibly could, since you overwrite it for every pixel of the input — you should eliminate this line.

In general, when you're trying to make something run fast, you should go over each line with a critical eye and ask, "How expensive is this line? What benefit does this line get me?" If it's expensive and has no benefit, then remove or rewrite it.

\$\endgroup\$
  • \$\begingroup\$ Please check out my edit. \$\endgroup\$ – Prikshet Sharma Aug 3 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.