Harmozing the API of the usual and gpu octree to faciliate swapping implementations easily

When working with the GPU and CPU versions of the octree class, I had to make several modifications to my codebase before running the same data on both host and device. In this ticket, I would first like to document the necessary changes based on the (abridged) octree tutorial. Furthermore, I thought of a few ways to harmonize APIs to facilitate swapping different implementations easily.

The (abridged) tutorial I run was:

#include <pcl/octree/octree_search.h>
#include <pcl/point_cloud.h>

#include <ctime>
#include <iostream>
#include <vector>

int main() {
    srand((unsigned int)time(NULL));

    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(
        new pcl::PointCloud<pcl::PointXYZ>);

    // Generate pointcloud data
    cloud->width = 1000;
    cloud->height = 1;
    cloud->points.resize(cloud->width * cloud->height);

    for (std::size_t i = 0; i < cloud->size(); ++i) {
        (*cloud)[i].x = 1024.0f * rand() / (RAND_MAX + 1.0f);
        (*cloud)[i].y = 1024.0f * rand() / (RAND_MAX + 1.0f);
        (*cloud)[i].z = 1024.0f * rand() / (RAND_MAX + 1.0f);
    }

    float resolution = 128.0f;

    pcl::octree::OctreePointCloudSearch<pcl::PointXYZ> octree(resolution);

    octree.setInputCloud(cloud);
    octree.addPointsFromInputCloud();

    pcl::PointXYZ searchPoint;

    searchPoint.x = 1024.0f * rand() / (RAND_MAX + 1.0f);
    searchPoint.y = 1024.0f * rand() / (RAND_MAX + 1.0f);
    searchPoint.z = 1024.0f * rand() / (RAND_MAX + 1.0f);

    std::vector<int> pointIdxRadiusSearch;
    std::vector<float> pointRadiusSquaredDistance;

    float radius = 256.0f * rand() / (RAND_MAX + 1.0f);

    std::cout << "Neighbors within radius search at (" << searchPoint.x << " "
              << searchPoint.y << " " << searchPoint.z
              << ") with radius=" << radius << std::endl;

    if (octree.radiusSearch(searchPoint, radius, pointIdxRadiusSearch,
                            pointRadiusSquaredDistance) > 0) {
        for (std::size_t i = 0; i < pointIdxRadiusSearch.size(); ++i)
            std::cout << "    " << (*cloud)[pointIdxRadiusSearch[i]].x << " "
                      << (*cloud)[pointIdxRadiusSearch[i]].y << " "
                      << (*cloud)[pointIdxRadiusSearch[i]].z
                      << " (squared distance: " << pointRadiusSquaredDistance[i]
                      << ")" << std::endl;
    }
}

The GPU version including the necessary revisions is:

#include <pcl/octree/octree_search.h>
#include <pcl/point_cloud.h>
#include <pcl/gpu/octree/octree.hpp>

#include <ctime>
#include <iostream>
#include <vector>

int main() {
    srand((unsigned int)time(NULL));

    //pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(
        //new pcl::PointCloud<pcl::PointXYZ>);

    pcl::gpu::DeviceArray<pcl::PointXYZ> cloud(1000);
    std::vector<pcl::PointXYZ> tmp(1000);

    // Generate pointcloud data
    //cloud->width = 1000;
    //cloud->height = 1;
    //cloud->points.resize(cloud->width * cloud->height);

    for (std::size_t i = 0; i < cloud.size(); ++i) {
        (tmp)[i].x = 1024.0f * rand() / (RAND_MAX + 1.0f);
        (tmp)[i].y = 1024.0f * rand() / (RAND_MAX + 1.0f);
        (tmp)[i].z = 1024.0f * rand() / (RAND_MAX + 1.0f);
    }
    cloud.upload(&tmp[0], 1000);

    //float resolution = 128.0f;
    //pcl::octree::OctreePointCloudSearch<pcl::PointXYZ> octree(resolution);
    pcl::gpu::Octree octree{};

    //octree.setInputCloud(cloud);
    octree.setCloud(cloud);
    //octree.addPointsFromInputCloud();
    octree.build();

    pcl::PointXYZ searchPoint;
    pcl::gpu::DeviceArray<pcl::PointXYZ> gpuSearchPointVec(1);

    searchPoint.x = 1024.0f * rand() / (RAND_MAX + 1.0f);
    searchPoint.y = 1024.0f * rand() / (RAND_MAX + 1.0f);
    searchPoint.z = 1024.0f * rand() / (RAND_MAX + 1.0f);
    gpuSearchPointVec.upload(&searchPoint, 1);

    // Neighbors within radius search

    //std::vector<int> pointIdxRadiusSearch;
    //std::vector<float> pointRadiusSquaredDistance;
    int max_answers = 1000;
    pcl::gpu::NeighborIndices results(gpuSearchPointVec.size(), max_answers);

    float radius = 256.0f * rand() / (RAND_MAX + 1.0f);

    std::cout << "Neighbors within radius search at (" << searchPoint.x << " "
              << searchPoint.y << " " << searchPoint.z
              << ") with radius=" << radius << std::endl;

    octree.radiusSearch(gpuSearchPointVec, radius, max_answers,
                            results);
}

Modification for GPU Octree
I could think of modifying the gpu::Octree module as follows:

• Add the member functions setInputCloud and addPointsFromInputCloud to the GPU Octree module
• Add another overload to the radiusSearch function returning the sequaredDistances too. A member function such as:
  void radiusSearch(const Queries& centers, float radius, int max_results, NeighborIndices& result, ResultSqrDists& sqr_distances) const might do that.
• Do we need the max_results argument? When I used Octree at first, I found it difficult to specify the number, and I am still unsure why we need it.
• Should we provide another overload that accepts one search point as a query instead of requiring a DeviceArray?

Modification for DeviceArray

• For the DeviceArray, I would like to add a subscript operator to load a host point to the device. At the moment, we have to load points into a std::vector and upload these points afterward.

What do others think about the example and the modifications?

Start by comparing pcl/octree/octree_search.h with pcl/gpu/octree/octree.hpp and the abridged octree tutorial in the issue. Review the GPU Octree, NeighborIndices, and DeviceArray APIs, then determine which proposed overloads and input-cloud changes have an agreed scope. Done should include a documented CPU/GPU usage comparison and tests for any accepted API changes.