[FEA] Implement SSPD (Symmetric Segment-Path Distance)

[mpearmain]

Is your feature request related to a problem? Please describe.
I wish cuSpatial could find the SSPD (Symmetric Segment-Path Distance) of trajectories
https://arxiv.org/abs/1508.04904

Describe the solution you'd like
As this is a distance measurement similar to Hausdorff (that is already implemented) having a new method which encapsulates distance metrics with the same signature would be ideal so that a user can easily switch between trajectory distance metrics i.e

import numpy as np
from cudf import Series
import cuspatial

in_trajs = []
in_trajs.append(np.array([[1, 0], [2, 1], [3, 2], [5, 3], [7, 1]]))
in_trajs.append(np.array([[0, 3], [2, 5], [3, 6], [6, 5]]))
in_trajs.append(np.array([[1, 4], [3, 7], [6, 4]]))
out_trajs = np.concatenate([np.asarray(traj) for traj in in_trajs], 0)
py_x = np.array(out_trajs[:, 0])
py_y = np.array(out_trajs[:, 1])
py_cnt = []
for traj in in_trajs:
    py_cnt.append(len(traj))
pnt_x = Series(py_x)
pnt_y = Series(py_y)
cnt = Series(py_cnt)

for distance in ["sspd",  "hausdorff"]:
    cuspatial.trajectory_distance(pnt_x, pnt_y, cnt, metric=distance)

Describe alternatives you've considered
I currently use the traj_dist package (https://github.com/bguillouet/traj-dist)

Additional context
My specific use case is for trajectory clustering of maritime vessels and the routes (trajectories) that they undertake from specific ports. (please find attached a small example of just two of the n clusters i have form this port)

In this case i have 11,000 trajectories i want to compute the SSPD distance for, i have around 8,000 ports in total and need to run for every port.

In the current setup i have a geopandas dataframe with a col called 'trajectories` and in each row i have an numpy array representing the points in a linestring

this leads to a line such as

import traj_dist.distance as tdist
p_dist = tdist.pdist(np.array([np.array(xi) for xi in df['clipped_trajectory']]), metric="sspd")

This currently takes around 1800 secs

if i were to try and find similar routes without subsetting on ports, it would lead to around 200,000 trajectory comparisons, on a single threaded machine this would likely take longer than the heat death of the universe :)

[thomcom]

Thanks @mpearmain ! This is a great and well structured feature request. We'll talk about it at our next meeting. We don't have new features planned for the next few weeks, but I think it would be an interesting target for many team members.

[github-actions]

This issue has been marked stale due to no recent activity in the past 30d. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed. This issue will be marked rotten if there is no activity in the next 60d.

[github-actions]

This issue has been marked rotten due to no recent activity in the past 90d. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed.

[mpearmain]

Coming back to this , as im now working on a problem in this area again.

[github-actions]

This issue has been labeled inactive-30d due to no recent activity in the past 30 days. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed. This issue will be labeled inactive-90d if there is no activity in the next 60 days.

[github-actions]

This issue has been labeled inactive-90d due to no recent activity in the past 90 days. Please close this issue if no further response or action is needed. Otherwise, please respond with a comment indicating any updates or changes to the original issue and/or confirm this issue still needs to be addressed.

[mpearmain]

A few years later - does something like this solve?

#include <cuda_runtime.h>
#include <cuspatial/cuspatial.h>
#include <iostream>

__device__ float point_distance(float2 p1, float2 p2) {
    float dx = p1.x - p2.x;
    float dy = p1.y - p2.y;
    return sqrtf(dx * dx + dy * dy);
}

__device__ float segment_point_distance(float2 a, float2 b, float2 p) {
    float2 ab = make_float2(b.x - a.x, b.y - a.y);
    float2 ap = make_float2(p.x - a.x, p.y - a.y);
    float t = fmaxf(0, fminf(1, (ap.x * ab.x + ap.y * ab.y) / (ab.x * ab.x + ab.y * ab.y)));
    float2 projection = make_float2(a.x + t * ab.x, a.y + t * ab.y);
    return point_distance(projection, p);
}

__global__ void compute_sspd(float2* segments1, int len1, float2* segments2, int len2, float* result) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < len1 * len2) {
        int i = idx / len2;
        int j = idx % len2;

        float2 a1 = segments1[2 * i];
        float2 b1 = segments1[2 * i + 1];
        float2 a2 = segments2[2 * j];
        float2 b2 = segments2[2 * j + 1];

        float d1 = segment_point_distance(a1, b1, a2);
        float d2 = segment_point_distance(a1, b1, b2);
        float d3 = segment_point_distance(a2, b2, a1);
        float d4 = segment_point_distance(a2, b2, b1);

        float min_dist = fminf(fminf(d1, d2), fminf(d3, d4));

        atomicMin(result, min_dist);
    }
}

#include <cuspatial/cuspatial.h>
#include <rmm/device_vector.hpp>
#include <iostream>

void compute_sspd_cuspatial(std::vector<float2> segments1, std::vector<float2> segments2) {
    // Error checking omitted for brevity

    int len1 = segments1.size() / 2;
    int len2 = segments2.size() / 2;
    float initial_distance = FLT_MAX;
    
    rmm::device_vector<float2> d_segments1(segments1);
    rmm::device_vector<float2> d_segments2(segments2);
    rmm::device_vector<float> d_result(1, initial_distance);

    int threadsPerBlock = 256;
    int blocksPerGrid = (len1 * len2 + threadsPerBlock - 1) / threadsPerBlock;

    compute_sspd<<<blocksPerGrid, threadsPerBlock>>>(d_segments1.data().get(), len1, d_segments2.data().get(), len2, d_result.data().get());

    cudaDeviceSynchronize();

    float result;
    cudaMemcpy(&result, d_result.data().get(), sizeof(float), cudaMemcpyDeviceToHost);

    std::cout << "SSPD: " << result << std::endl;
}

int main() {
    std::vector<float2> segments1 = {make_float2(0.0f, 0.0f), make_float2(1.0f, 1.0f), make_float2(1.0f, 0.0f), make_float2(2.0f, 1.0f)};
    std::vector<float2> segments2 = {make_float2(0.0f, 1.0f), make_float2(1.0f, 2.0f), make_float2(2.0f, 0.0f), make_float2(3.0f, 1.0f)};

    compute_sspd_cuspatial(segments1, segments2);

    return 0;
}

[harrism]

Conceptually, yes. I wouldn't do a global atomicMin like that, because all threads will be bottlenecking on a single memory location. Better to use a thrust::transform_reduce and trust in Thrust to do it efficiently (in the long run it will be more consistently faster than code you can write across a variety of hardware).

[mpearmain]

Thanks for the pro tip @harrism

im guessing something like this (thanks CudaGPT!):

#include <thrust/transform_reduce.h>
#include <thrust/device_vector.h>
#include <rmm/device_vector.hpp>
#include <iostream>
#include <limits>

__device__ float point_distance(float2 p1, float2 p2) {
    float dx = p1.x - p2.x;
    float dy = p1.y - p2.y;
    return sqrtf(dx * dx + dy * dy);
}

__device__ float segment_point_distance(float2 a, float2 b, float2 p) {
    float2 ab = make_float2(b.x - a.x, b.y - a.y);
    float2 ap = make_float2(p.x - a.x, p.y - a.y);
    float t = fmaxf(0, fminf(1, (ap.x * ab.x + ap.y * ab.y) / (ab.x * ab.x + ab.y * ab.y)));
    float2 projection = make_float2(a.x + t * ab.x, a.y + t * ab.y);
    return point_distance(projection, p);
}

struct SSPD_Functor {
    float2 *segments1, *segments2;
    int len2;

    SSPD_Functor(float2 *s1, float2 *s2, int l2) : segments1(s1), segments2(s2), len2(l2) {}

    __device__ float operator()(int idx) const {
        int i = idx / len2;
        int j = idx % len2;

        float2 a1 = segments1[2 * i];
        float2 b1 = segments1[2 * i + 1];
        float2 a2 = segments2[2 * j];
        float2 b2 = segments2[2 * j + 1];

        float d1 = segment_point_distance(a1, b1, a2);
        float d2 = segment_point_distance(a1, b1, b2);
        float d3 = segment_point_distance(a2, b2, a1);
        float d4 = segment_point_distance(a2, b2, b1);

        return fminf(fminf(d1, d2), fminf(d3, d4));
    }
};

void compute_sspd_cuspatial(std::vector<float2> &segments1, std::vector<float2> &segments2) {
    int len1 = segments1.size() / 2;
    int len2 = segments2.size() / 2;
    int num_pairs = len1 * len2;

    rmm::device_vector<float2> d_segments1(segments1);
    rmm::device_vector<float2> d_segments2(segments2);

    thrust::device_vector<int> indices(num_pairs);
    thrust::sequence(indices.begin(), indices.end());

    float max_distance = std::numeric_limits<float>::max();

    float result = thrust::transform_reduce(
        indices.begin(), indices.end(),
        SSPD_Functor(thrust::raw_pointer_cast(d_segments1.data()), thrust::raw_pointer_cast(d_segments2.data()), len2),
        max_distance, thrust::minimum<float>());

    return result;
}

What do you think about having an encapsulating function for trajectory distance in cuSpatial as per original idea? if this isn't of interest, I will close this and implement the above for my purposes

for distance in ["sspd",  "hausdorff"]:
    cuspatial.trajectory_distance(pnt_x, pnt_y, cnt, metric=distance)

[harrism]

We have definitely discussed abstracting distance metrics. Feel free to open an issue (and if you like, a PR).

Currently development on cuSpatial is not very active.

[ikaryss]

@mpearmain thank you so much for the code you provided!! It is really helpful. I'm currently working on a similar trajectory clustering problem and have been studying your implementation to better understand it

In the paper you referenced, the SSPD (Symmetric Segment Path Distance) is calculated in several steps:

1. calculate D_spd(T1, T2), which is the mean of all minimum distances between points from trajectory T1 to segments from trajectory T2.
2. calculate D_spd(T2, T1) in a similar manner.
3. Take the mean of these two D_spd values to obtain the SSPD.

In your implementation, however, it appears that you calculate the SSPD as the minimum value of the minimum distances between all segments from the two trajectories (if I understand correctly).

Could you please explain why you chose to calculate the SSPD in this manner? I’m curious about the reasoning behind this approach.