Point Cloud
(top) (pkg)

The PointCloud::Array is a simple container of array like data. It’s primary purpose:

• Be held by Dataset.
• Erase the C++ numeric type and shape information of the underlying array to enable Dataset to be heterogeneous.
• Allow option to share user array data or to maintain an internal copy.
• Support appending along the major axis.
• Supply array views in useful types (flat vector as boost::span<T> or shaed boost::multi_array<T,NDIM>).

Here are more creation examples:

// some user data
std::vector<int> user(10, 0);

// shape: (2 rows, 5 columns), user data is shared
Array arr1(user, {2,5}, true);

// special, shape: (10,) not shared
Array arr2(user);

// copy/move constructor/assignment all supported
arr1 = arr2;

The underlying data can be retrieved in a number of constant, read-only forms that are more useful for custom code operations. These forms include:

// Return a boost::span (same as C++20 std::span and that is
// std::vector like), of the flattened array.
auto flat = arr1.elements<int>();

// as properly shapped boost::multi_array
auto ma = arr1.indexed<int, 2>();

// fixme: is an Eigen3 array interface useful?

These formats do not impose a copy and they will throw a ValueError if the requested type or shape is not compatible with the underlying array content.

The contents of an Array may be fully replaced which drops and does not modify previous contents:

// Just as when we initially constructed arr1 above.
arr1.assign(user.data(), {2,5}, true);

// array of shape (10,) not shared
arr1.assign(user);

One of the few array-like operations supported is to append to an Array. If the underlying data is shared, the append() methods will first copy that data to an internal buffer and drop the original shared array (copy-on-write semantics).

// Reuse user data 
Array tail(user, {2,5}, true);

// append, causes internal copy so "user" not modified.
arr1.append(user);

The array begin appended (tail) must be of the same numeric type as the current array and the sizes of all non-major axes of tail must match those same axes of the current array. The following represents shapes which are legal to use in an append and their resulting shape:

(N,n2,n3,...) + (M,n2,n3,...) = (N+M,n2,n3,...)

And Array also carries but does not in any way use a metadata object implemented as WireCell::Configuration (a JsonCPP object).

Array a;
auto& md = a.metadata();
md["foo"] = "bar";  

A PointCloud::Dataset provides a collection of Array with minimal operations. It’s primary purposes:

• Hold instances of Array, each identified by a “name” of type std::string.
• Assure all accepted Array have common sized major axes.
• Implement the appending of another Dataset while keeping assurances.
• Accept user-provided callback hooks to be called on successful append.
• Provide various means of access to the collection of the arrays.

Here shows the successful adding of two arrays followed by a rejection of a misshapen third.

Dataset d;
d.add("one", Array({1,2,3,4,5}));
d.add("two", Array({1.1,2.2,3.3,4.4,5.5}));
d.add("broken", Array({1}));    // throws ValueError

One may add() higher dimension arrays as long as the number of their “elements” match. In general, for a point data of shape (v1, ...) the point data array must be of shape (nele, v1, ...) to be added to a dataset with num_elements() returning nele.

The names of the held arrays can be retrieved:

for (const auto& name : d.keys()) {
    // ...
}

One Dataset (“tail”) may be appended to another Dataset (“head”). This operation merely appends each Array in “tail” to the array of the same name in “head”. A ValueError is thrown if any names in “head” are not in “tail” or if any in the set of matching arrays from “tail” have differing number of elements.

Dataset d2 = ...;
d.append(d2);

A Dataset may be copy/move assigned/constructed.

Dataset d2 = d;
Dataset d3(d);

Like Array, Dataset carries but does not use a metadata object.

Dataset d;
auto& md = d.metadata();
md["foo"] = "baz";

A Dataset is mapped to an ITensorSet and each of its named Array items to an ITensor in the set. The two representations are close, but not exact.

• On creation, the ItensorSet metadata will be given a special attribute "_dataset_arrays" holding an array of strings representing the names of the Array items and in the same order as the ITensor vector which the set holds.
• When read, if this special attribute is missing then a "name" attribute from the metadata of each ITensor is used to provide the Array name in the Dataset. If this attribute is missing a name is generated as "array%d" which will be interpolated with the index of the corresponding ITensor.
• The ident value of the ITensorSet may be provided by an "ident" attribute of the Dataset metadata, else zero is stored. In the reverse direction, the Dataset will be given an "ident" metadata attribute holding the value from the ITensorSet.
• Otherwise, all Array and Dataet metadata objects are pass unchanged to/from those of ITensor and ITensorSet.

Each of the four possible conversions has a corresponding function. Their interface is very simple and the WireCellAux/TensorTools.h header file gives the essentials:

Note that share may be passed as true to enable the optimization of zero-copy with optional copy-on-write. It is false by default as the user must assure that the ITensorSet or ITensor remains alive.

The aux/test/test_tensor_tools.cxx unit test contains more examples.

A nanoflann dataset adaptor must be provided for each unique type of PC data. One that makes use of Eigen3 arrays is provided by nanoflann. WCT provides one for boost::multi_array that comes from an ITensor as described above. A dataset adaptor must provide two methods and may provide a third, optional method. As copied from comments in the source and embelished these are:

template<...>
class MyDatasetAdaptor {
public:
    using point_type = ...;

    // Must return the number of data poins
    inline size_t kdtree_get_point_count() const
    {
        // return npoints
    }

    // Must return the dim'th component of the idx'th point in the
    // class:
    inline point_type kdtree_get_pt(const size_t idx, const size_t dim) const
    {
        // return a point
    }

    // Optional bounding-box computation: return false to default to a
    // standard bbox computation loop.  Return true if the BBOX was
    // already computed by the class and returned in "bb" so it can be
    // avoided to redo it again.  Look at bb.size() to find out the
    // expected dimensionality (e.g. 2 or 3
    template <class BBOX>
    bool kdtree_get_bbox(BBOX& bb) const
    {
        bb[0].low = ...; bb[0].high = ...;  // 0th dimension limits
        bb[1].low = ...; bb[1].high = ...;  // 1st dimension limits
        // ...
        return true;
    }
};

A distance metric must be provided to measure separation of points in the space. L1, L2 and some other metrics are provided. Note, the units of “distance” depend on the metric. For example, when the L2 metric is used the distance is in units of length-squared.

A nanoflann index adaptor provides the k-d tree query interface. KDTreeSingleIndexAdaptor provides one implementation. It is templated in terms of a distance type, a dataset adaptor, the dimensionality and the index type. It operates assuming the adapted dataset remains unchanged.

The KDTreeSingleIndexDynamicAdaptor is another which allows for points to be added or removed from the k-d tree. Points may be removed from the k-d tree without removing them from the dataset. When adding points, they are first added to the dataset and the adaptor is notified about the range of indices which are added.

The nanoflann library also defines an KDTreeEigenMatrixAdaptor. This class shears across both dataset adaptor and index adaptor design layers. It provides a dataset adaptor for Eigen::Matrix and then internally creates an maintains a KDTreeSingleIndexAdaptor using itself as the dataset adaptor.

There are two primary query methods provided by KDTree*Index*Adaptor classes (and via the exposed KDTreeEigenMatrixAdaptor::index pointer to a KDTreeSingleIndexAdaptor).

Size knnSearch(
    const ElementType* query_point, const Size num_closest,
    AccessorType* out_indices, DistanceType* out_distances_sq) const;

Size radiusSearch(
    const ElementType* query_point, const DistanceType& radius,
    std::vector<std::pair<AccessorType, DistanceType>>& IndicesDists,
    const SearchParams&                                 searchParams) const;

The knnSearch() finds the k-nearest neighbors (num_closest) of the query point and radiusSearch() finds neighbors withing the given distance.

The ElementType* query_point provides a point in the space as a C array. The type ElementType is the type of the point coordinates. The DistanceType is the type in which distances between points is expressed. These are set through the distance metric type and by default they are the same type.

The AccessorType is that of the index of a point in the point cloud. The SearchParams are a small bundle, presumably to tune the k-d tree.

Both of these rely on findNeighbors() method.

template <typename RESULTSET>
bool findNeighbors(
    RESULTSET& result, const ElementType* query_point,
    const SearchParams& searchParams) const;

The RESULTSET passed to findNeighbors() may be a KNNResultSet<DistanceType> for knnSearch() or a RadiusResultSet<DistanceType> for radiusSearch(). The only reason to call this directly instead of through the high level methods is if the result set provides additional required information that they do not.