Classify Example Source Files

Main.cpp

  1/*
  2 * Copyright (c) 2014-2023 National Technology and Engineering
  3 * Solutions of Sandia, LLC. Under the terms of Contract DE-NA0003525
  4 * with National Technology and Engineering Solutions of Sandia, LLC,
  5 * the U.S. Government retains certain rights in this software.
  6 *
  7 * Redistribution and use in source and binary forms, with or without
  8 * modification, are permitted provided that the following conditions
  9 * are met:
 10 *
 11 * 1. Redistributions of source code must retain the above copyright
 12 * notice, this list of conditions and the following disclaimer.
 13 *
 14 * 2. Redistributions in binary form must reproduce the above copyright
 15 * notice, this list of conditions and the following disclaimer in the
 16 * documentation and/or other materials provided with the distribution.
 17 *
 18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 22 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 29 */
 30
 31#include "AssignHeadings.h"
 32#include "Mapping.h"
 33#include "TrackFilter.h"
 34
 35#include <tracktable/CommandLineFactories/AssemblerFromCommandLine.h>
 36#include <tracktable/CommandLineFactories/PointReaderFromCommandLine.h>
 37#include <tracktable/Core/Geometry.h>
 38#include <tracktable/Domain/Terrestrial.h>
 39#include <tracktable/RW/KmlOut.h>
 40
 41#include <boost/timer/timer.hpp>
 42
 43#include <algorithm>
 44
 45using TrajectoryT = tracktable::domain::terrestrial::trajectory_type;
 46using PointT = typename TrajectoryT::point_type;
 47using PointReaderT = tracktable::PointReader<PointT>;
 48using PointReaderIteratorT = typename PointReaderT::iterator;
 49using AssemblerT = tracktable::AssembleTrajectories<TrajectoryT, PointReaderIteratorT>;
 50constexpr auto length = tracktable::length<TrajectoryT>;
 51
 52using tracktable::kml;
 53
 54static constexpr auto helpmsg = R"(
 55--------------------------------------------------------------------------------
 56The classify example demonstrates:
 57    - Using command line factories to read points and assemble trajectories
 58    - Using boost program options to take parameters from command lines(in addition to the factories)
 59    - Filtering trajectories on any combination of the following:
 60      - length
 61      - curvature
 62      - hull-gyration ratio
 63      - length ratio
 64      - hull-aspect ratio
 65      - straightness
 66      - number of turn arounds
 67    - A few different methods of applying these filters are demonstrated int he code
 68    - Writing trajectories as KML
 69
 70Typical use:
 71    ./ classify-- input=/data/flights.tsv --min-turn-arounds=10 --output =/results/mappingflights.kml
 72
 73Defaults assume a tab separated file formatted as :
 74
 75OBJECTID TIMESTAMP LON LAT
 76
 77Default output is standard out
 78--------------------------------------------------------------------------------)";
 79
 80int main(int _argc, char* _argv[]) {
 81    tracktable::set_log_level(tracktable::log::info);
 82    boost::program_options::options_description commandLineOptions;
 83    commandLineOptions.add_options()("help", "Print help");
 84    boost::program_options::options_description classifyOptions("Classify");
 85    // clang-format off
 86    /** Here we manually add some filter options
 87    *     This is a syntax created by boost program options using some clever tricks.
 88    *     This is legal C++ but it should look weird to most users. */
 89    classifyOptions.add_options()
 90      ("assign-headings",
 91      "Assign headings to points")
 92      ("min-length",
 93        bpo::value<double>(), //We intentionally leave off default, this allows us to skip filtering if no value
 94      "minimum length of trajectory")
 95      ("max-length",
 96        bpo::value<double>(),
 97      "maximum length of trajectory")
 98      ("min-curvature",
 99        bpo::value<double>(),
100      "minimum curvature of trajectory")
101      ("max-curvature",
102        bpo::value<double>(),
103      "maximum curvature of trajectory")
104    ;
105    // clang-format on
106    /** We can use the Trackfilter class to help us with the repetitive task of filtering
107     *   We provide the TrackFilter with a lambda that corresponds to the 'ShouldKeep'
108     *   scenario the filter. */
109    MinMaxTrackFilter<double> hullGyrationFilter("hull-gyration-ratio", [](const TrajectoryT& _t) {
110        return tracktable::convex_hull_area(_t) / tracktable::radius_of_gyration(_t);
111    });
112    hullGyrationFilter.addOptions(classifyOptions);
113
114    MinMaxTrackFilter<double> lengthRatioFilter("length-ratio", [](const TrajectoryT& _t) {
115        return tracktable::end_to_end_distance(_t) / tracktable::length(_t);
116    });
117    lengthRatioFilter.addOptions(classifyOptions);
118
119    MinMaxTrackFilter<double> hullAspectRatioFilter("hull-aspect-ratio",
120                                                    tracktable::convex_hull_aspect_ratio<TrajectoryT>);
121    hullAspectRatioFilter.addOptions(classifyOptions);
122
123    MinMaxTrackFilter<double> straightnessFilter("straightness", StraightFraction);
124    straightnessFilter.addOptions(classifyOptions);
125
126    MinMaxTrackFilter<unsigned int> turnAroundsFilter("turn-arounds", TurnArounds);
127    turnAroundsFilter.addOptions(classifyOptions);
128
129    commandLineOptions.add(classifyOptions);
130
131    // Reader factories, which add lots of command line options
132    tracktable::PointReaderFromCommandLine<PointT> readerFactory;
133    tracktable::AssemblerFromCommandLine<TrajectoryT> assemblerFactory;
134    readerFactory.addOptions(commandLineOptions);
135    assemblerFactory.addOptions(commandLineOptions);
136
137    // By creating separate program options objects, we are creating groups for the
138    //   help menu.
139    boost::program_options::options_description outputOptions("Output Options");
140    // clang-format off
141    // We manually add options for output
142    outputOptions.add_options()
143      ("no-output",
144      "specifies no output is wanted")
145      ("separate-kmls",
146      "indicate whether to separate output to different kmls in [result-dir]")
147      ("result-dir",
148        bpo::value<std::string>()->default_value("result/"),
149      "directory to store separate kmls")
150      ("output",
151        bpo::value<std::string>()->default_value("-"),
152        "file to write to (use '-' for stdout), overridden by 'separate-kmls'"
153      )
154    ;
155    // clang-format on
156    commandLineOptions.add(outputOptions);
157
158    /** Boost program options using a variable map to tie everything together.
159     * one parse will have a single variable map. We need to let the factories know
160     * about this variable map so they can pull information out of it */
161    auto vm = std::make_shared<boost::program_options::variables_map>();
162    readerFactory.setVariables(vm);
163    assemblerFactory.setVariables(vm);
164
165    // Parse the command lines, don't forget the 'notify' after
166    try {
167        // We use this try/catch to automatically display help when an unknown option is used
168        boost::program_options::store(
169            boost::program_options::command_line_parser(_argc, _argv).options(commandLineOptions).run(), *vm);
170        boost::program_options::notify(*vm);
171    } catch (boost::program_options::error e) {
172        std::cerr << e.what();
173        std::cerr << helpmsg << "\n\n";
174        std::cerr << commandLineOptions << std::endl;
175        return 1;
176    }
177    /** Parsing will give an error of an incorrect option is used, but it won't
178     * display the help unless we tell it too */
179    if (vm->count("help") != 0) {
180        std::cerr << helpmsg << "\n\n";
181        std::cerr << commandLineOptions << std::endl;
182        return 1;
183    }
184
185    // Create Point Reader and assembler
186    auto pointReader = readerFactory.createPointReader();
187    auto assembler = assemblerFactory.createAssembler(pointReader);
188
189    std::vector<TrajectoryT> trajectories = {};
190    // This block exists for easy timing of trajectory assembling using the boost auto timer
191    // Note that all feedback to the user is done on std::cerr, this allows us to only
192    // put desired results into std::cout, this make processing output easier.
193    {
194        std::cerr << "Assemble Trajectories" << std::endl;
195        boost::timer::auto_cpu_timer timer3(std::cerr);
196        auto count = 0u;
197        std::cerr << std::right;
198        for (auto tIter = assembler->begin(); tIter != assembler->end(); ++tIter) {
199            std::cerr << "\b\b\b\b\b\b\b\b\b\b" << std::setw(10)  // Using backspaces for in place counter
200                      << count++;
201            trajectories.push_back(*tIter);
202        }
203        std::cerr << std::left << "\nStarting with " << trajectories.size() << " trajectories" << std::endl;
204    }
205
206    // Filter tracks based on length
207    // This is the typical method, fairly compact and reasonably efficient
208    // we can check if an option was used by checking its 'count' in the variable map
209    if (vm->count("min-length") != 0 || vm->count("max-length") != 0) {
210        std::cerr << "Filtering based on length" << std::endl;
211        boost::timer::auto_cpu_timer t(std::cerr);
212        auto lower = vm->count("min-length") != 0 ? (*vm)["min-length"].as<double>() : 0.0;
213        auto upper = vm->count("max-length") != 0 ? (*vm)["max-length"].as<double>()
214                                                  : std::numeric_limits<double>::max();
215        // We use remove_if to move all trajectories where the lambda returns true to the end
216        // remove_if returns the end of the list of kept items.
217        // we then erase everything from that point to the 'end()'
218        trajectories.erase(std::remove_if(trajectories.begin(), trajectories.end(),
219                                          [&](const TrajectoryT& _t) {
220                                              auto l = length(_t);
221                                              return (upper < l) || (l < lower);
222                                          }),
223                           trajectories.end());
224        std::cerr << trajectories.size() << " trajectories after filtering for length" << std::endl;
225    }
226
227    // Filter tracks based on curvature
228    // This can be used to looks for loops.
229    // This method is identical to that above but the steps are seperated out for demonstration
230    if (vm->count("min-curvature") != 0 || vm->count("max-curvature") != 0) {
231        std::cerr << "Filtering based on curvature" << std::endl;
232        boost::timer::auto_cpu_timer t(std::cerr);
233        auto lower = vm->count("min-curvature") != 0 ? (*vm)["min-curvature"].as<double>() : 0.0;
234        auto upper = vm->count("max-curvature") != 0 ? (*vm)["max-curvature"].as<double>()
235                                                     : std::numeric_limits<double>::max();
236        // Create and save a lambda to use as filter predicate
237        auto filterFunc = [&](const TrajectoryT& _t) {
238            auto c = std::abs(TotalCurvature(_t));
239            return (upper < c) || (c < lower);
240        };
241        // we save the result of remove if to use in the erase
242        auto lastGood = std::remove_if(trajectories.begin(), trajectories.end(), filterFunc);
243        // we then erase everything from that point to the 'end()'
244        trajectories.erase(lastGood, trajectories.end());
245        std::cerr << trajectories.size() << " trajectories after filtering for curvature" << std::endl;
246    }
247
248    // use track filter object, we still manually check options ot skip if not asked for.
249    // The track filter object uses the same method as above, but centralizes the code
250    // for improved reliability and maintenance.
251    if (vm->count("min-hull-gyration-ratio") != 0 || vm->count("max-hull-gyration-ratio") != 0) {
252        auto& filter = hullGyrationFilter;
253        std::cerr << "Filtering based on " << filter.getName() << std::endl;
254        boost::timer::auto_cpu_timer t(std::cerr);
255        filter(trajectories);
256        std::cerr << trajectories.size() << " trajectories after filtering for " << filter.getName()
257                  << std::endl;
258    }
259
260    // This ratio represents the directness of a flight
261    if (vm->count("min-length-ratio") != 0 || vm->count("max-length-ratio") != 0) {
262        auto& filter = lengthRatioFilter;
263        std::cerr << "Filtering based on " << filter.getName() << std::endl;
264        boost::timer::auto_cpu_timer t(std::cerr);
265        filter(trajectories);
266        std::cerr << trajectories.size() << " trajectories after filtering for " << filter.getName()
267                  << std::endl;
268    }
269
270    if (vm->count("min-hull-aspect-ratio") != 0 || vm->count("max-hull-aspect-ratio") != 0) {
271        auto& filter = hullAspectRatioFilter;
272        std::cerr << "Filtering based on " << filter.getName() << std::endl;
273        boost::timer::auto_cpu_timer t(std::cerr);
274        filter(trajectories);
275        std::cerr << trajectories.size() << " trajectories after filtering for " << filter.getName()
276                  << std::endl;
277    }
278    // Apply headings if needed
279    if (0 != (vm->count("assign-headings") + vm->count("min-straightness") + vm->count("max-straightness") +
280              vm->count("min-turn-arounds") + vm->count("max-turn-arounds"))) {
281        AssignTrajectoryHeadings(trajectories);
282    }
283    // Mapping flights will have a high straightness with a high number of turn arounds
284    if (vm->count("min-straightness") != 0 || vm->count("max-straightness") != 0) {
285        auto& filter = straightnessFilter;
286        std::cerr << "Filtering based on " << filter.getName() << std::endl;
287        boost::timer::auto_cpu_timer t(std::cerr);
288        filter(trajectories);
289        std::cerr << trajectories.size() << " trajectories after filtering for " << filter.getName()
290                  << std::endl;
291    }
292    if (vm->count("min-turn-arounds") != 0 || vm->count("max-turn-arounds") != 0) {
293        auto& filter = turnAroundsFilter;
294        std::cerr << "Filtering based on " << filter.getName() << std::endl;
295        boost::timer::auto_cpu_timer t(std::cerr);
296        filter(trajectories);
297        std::cerr << trajectories.size() << " trajectories after filtering for " << filter.getName()
298                  << std::endl;
299    }
300
301    if (0 != vm->count("no-output")) {
302        std::cerr << "No Output" << std::endl;
303        return 0;
304    }
305    if (0 != vm->count("separate-kmls")) {
306        auto resultsDirectory = (*vm)["result-dir"].as<std::string>();
307        std::cerr << "Writing separate kml files to " << resultsDirectory << std::endl;
308        kml::writeToSeparateKmls(trajectories, resultsDirectory);
309        return 0;
310    }
311    auto filename = (*vm)["output"].as<std::string>();
312    if ("-" != filename) {
313        std::ofstream outfile(filename);
314        std::cerr << "Writing to " << filename << std::endl;
315        outfile << kml(trajectories);
316        outfile.close();
317        return 0;
318    }
319    std::cerr << "Writing to stdout" << std::endl;
320    std::cout << kml::header;
321    std::cout << kml(trajectories);
322    std::cout << kml::footer;
323    return 0;
324}

AssignHeadings.h

 1/*
 2 * Copyright (c) 2014-2023 National Technology and Engineering
 3 * Solutions of Sandia, LLC. Under the terms of Contract DE-NA0003525
 4 * with National Technology and Engineering Solutions of Sandia, LLC,
 5 * the U.S. Government retains certain rights in this software.
 6 *
 7 * Redistribution and use in source and binary forms, with or without
 8 * modification, are permitted provided that the following conditions
 9 * are met:
10 *
11 * 1. Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 *
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
22 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31#ifndef AssignHeadings_h
32#define AssignHeadings_h
33
34#include <tracktable/Core/detail/algorithm_signatures/Bearing.h>
35#include <tracktable/Core/detail/algorithm_signatures/TurnAngle.h>
36
37#include <algorithm>
38#include <vector>
39
40template <typename TrajectoryT>
41void AssignTrajectoryHeadings(TrajectoryT &trajectory) {
42    if (trajectory.size() == 0) return;
43    if (trajectory.size() == 1) {
44        trajectory[0].set_property("heading", 0.0);
45        return;
46    }
47    for (unsigned int i = 0; i < trajectory.size() - 1; ++i) {
48        trajectory[i].set_property("heading", tracktable::bearing(trajectory[i], trajectory[i + 1]));
49    }
50    trajectory[trajectory.size() - 1].set_property(
51        "heading", trajectory[trajectory.size() - 2].real_property("heading"));
52}
53// overload for vector type
54template <typename TrajectoryT>
55void AssignTrajectoryHeadings(std::vector<TrajectoryT> &trajectories) {
56    // use static cast to specify which overload
57    std::for_each(trajectories.begin(), trajectories.end(),
58                  static_cast<void (*)(TrajectoryT &)>(AssignTrajectoryHeadings<TrajectoryT>));
59}
60
61template <typename TrajectoryT>
62double TotalCurvature(TrajectoryT const &trajectory) {
63    if (trajectory.size() < 3) {
64        return 0.0;
65    }
66    auto curvature = 0.0;
67    for (auto i = 1u; i < trajectory.size() - 1; ++i)
68        curvature += signed_turn_angle(trajectory[i - 1], trajectory[i], trajectory[i + 1]);
69    return curvature;
70}
71
72#endif

Mapping.h

 1/*
 2 * Copyright (c) 2014-2023 National Technology and Engineering
 3 * Solutions of Sandia, LLC. Under the terms of Contract DE-NA0003525
 4 * with National Technology and Engineering Solutions of Sandia, LLC,
 5 * the U.S. Government retains certain rights in this software.
 6 *
 7 * Redistribution and use in source and binary forms, with or without
 8 * modification, are permitted provided that the following conditions
 9 * are met:
10 *
11 * 1. Redistributions of source code must retain the above copyright
12 * notice, this list of conditions and the following disclaimer.
13 *
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
22 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31#ifndef Mapping_h
32#define Mapping_h
33
34#include <tracktable/Domain/Terrestrial.h>
35
36// TODO: Use template to get rid of this
37using TrajectoryT = tracktable::domain::terrestrial::trajectory_type;
38using PointT = typename TrajectoryT::point_type;
39
40unsigned int TurnArounds(TrajectoryT const& trajectory);
41double StraightFraction(TrajectoryT const& trajectory);
42double HeadingDifference(const double h2, const double h1);
43double HeadingDifference(const PointT& t2, const PointT& t1);
44
45#endif

Mapping.cpp

  1/*
  2 * Copyright (c) 2014-2023 National Technology and Engineering
  3 * Solutions of Sandia, LLC. Under the terms of Contract DE-NA0003525
  4 * with National Technology and Engineering Solutions of Sandia, LLC,
  5 * the U.S. Government retains certain rights in this software.
  6 *
  7 * Redistribution and use in source and binary forms, with or without
  8 * modification, are permitted provided that the following conditions
  9 * are met:
 10 *
 11 * 1. Redistributions of source code must retain the above copyright
 12 * notice, this list of conditions and the following disclaimer.
 13 *
 14 * 2. Redistributions in binary form must reproduce the above copyright
 15 * notice, this list of conditions and the following disclaimer in the
 16 * documentation and/or other materials provided with the distribution.
 17 *
 18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 22 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 29 *
 30 * Mapping: A simple example for finding mapping flights with lots of
 31 * back-and-forth motion.
 32 *
 33 * Created by Danny Rintoul.
 34 */
 35
 36#include "Mapping.h"
 37
 38#include <algorithm>
 39
 40// TODO: Move this into tracktable, if there isn't already an equivalent
 41double HeadingDifference(const double h2, const double h1) {
 42    return (h2 - h1) - 360.0 * int((h2 - h1) / 180.0);
 43}
 44
 45double HeadingDifference(const PointT& t2, const PointT& t1) {
 46    auto h2 = t2.real_property("heading");
 47    auto h1 = t1.real_property("heading");
 48    return HeadingDifference(h2, h1);
 49}
 50
 51// TODO: investigate whether this is prudent without assigning headings.
 52// TODO: Move into tracktable.
 53unsigned int TurnArounds(TrajectoryT const& trajectory) {
 54    constexpr size_t window = 5;
 55    unsigned int ctr = 0;
 56
 57    if (trajectory.size() <= window) return 0;
 58
 59    auto itr1 = trajectory.begin();
 60    auto itr2 = itr1 + window;
 61    double diff = 0;
 62    bool found = false;
 63
 64    do {
 65        diff = std::abs(itr1->real_property("heading") - itr2->real_property("heading"));
 66
 67        if (std::abs(diff - 180.0) < 2.0) {  // 2.0 is the fudge factor for comparing against 180 degrees
 68            if (!found) {
 69                ++ctr;
 70            }
 71            found = true;
 72        } else {
 73            found = false;
 74        }
 75
 76        if (found) {
 77            int leap = std::min(static_cast<int>(std::distance(itr2, trajectory.end())), 5);
 78            // The 5 above is related to how far away you want to go before you
 79            // define another turnaround
 80            itr1 += leap;
 81            itr2 += leap;
 82        } else {
 83            ++itr1;
 84            ++itr2;
 85        }
 86    } while (itr2 != trajectory.end());
 87
 88    return ctr;
 89}
 90
 91// ----------------------------------------------------------------------
 92// TODO: possible to do this without assigning headings using unsigned_angle
 93double StraightFraction(TrajectoryT const& trajectory) {
 94    int sum = 0;
 95    constexpr size_t min_straight_size = 5;
 96
 97    auto itr1 = trajectory.begin();
 98    auto itr2 = trajectory.end();
 99
100    do {
101        // The 2.0 below is related to how different the straightness can be from 180
102        itr2 = std::adjacent_find(itr1, trajectory.end(), [](const PointT& _lhs, const PointT& _rhs) {
103            return std::abs(HeadingDifference(_lhs, _rhs)) < 2.0;
104        });
105        if (itr2 != trajectory.end()) {
106            ++itr2;
107        }
108
109        if (size_t(itr2 - itr1) >= min_straight_size) {
110            sum += (itr2 - itr1);
111        }
112
113        itr1 = itr2;
114    } while (itr2 != trajectory.end());
115
116    return static_cast<double>(sum) / static_cast<double>(trajectory.size());
117}

TrackFilter.h

  1/*
  2 * Copyright (c) 2014-2023 National Technology and Engineering
  3 * Solutions of Sandia, LLC. Under the terms of Contract DE-NA0003525
  4 * with National Technology and Engineering Solutions of Sandia, LLC,
  5 * the U.S. Government retains certain rights in this software.
  6 *
  7 * Redistribution and use in source and binary forms, with or without
  8 * modification, are permitted provided that the following conditions
  9 * are met:
 10 *
 11 * 1. Redistributions of source code must retain the above copyright
 12 * notice, this list of conditions and the following disclaimer.
 13 *
 14 * 2. Redistributions in binary form must reproduce the above copyright
 15 * notice, this list of conditions and the following disclaimer in the
 16 * documentation and/or other materials provided with the distribution.
 17 *
 18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 19 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 20 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 21 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 22 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 23 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 24 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 25 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 26 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 28 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 29 */
 30
 31#ifndef TrackFilter_h
 32#define TrackFilter_h
 33
 34#include <tracktable/Domain/Terrestrial.h>
 35
 36#include <boost/program_options.hpp>
 37
 38#include <limits>
 39
 40using TrajectoryT = tracktable::domain::terrestrial::trajectory_type;
 41using PointT = typename TrajectoryT::point_type;
 42// TODO: Template these
 43
 44/** @brief Utility interface that provides framework for filtering tracks
 45 * In it's simplest form the only advantage it provides is each of maintenance
 46 * of the erase loops.
 47 *
 48 * The real advantage comes when the command line automation is added in
 49 * the MinMaxTrackFilter implementation
 50 * @sa MinMaxTrackFilter
 51 */
 52class TrackFilter {
 53   public:
 54    TrackFilter(const std::string& _name) : name(_name) {}
 55    void operator()(std::vector<TrajectoryT>& _tracks) { filterTracks(_tracks); }
 56
 57    void filterTracks(std::vector<TrajectoryT>& _tracks) {
 58        // use lambda to wrap non-static function
 59        _tracks.erase(std::remove_if(_tracks.begin(), _tracks.end(),
 60                                     [&](const TrajectoryT& _t) { return this->shouldNotKeepTrack(_t); }),
 61                      _tracks.end());
 62    }
 63    void inverseFilterTracks(std::vector<TrajectoryT>& _tracks) {
 64        _tracks.erase(std::remove_if(_tracks.begin(), _tracks.end(),
 65                                     [&](const TrajectoryT& _t) { return this->shouldKeepTrack(_t); }),
 66                      _tracks.end());
 67    }
 68
 69    const std::string& getName() const { return name; }
 70
 71   protected:
 72    /** This is filter criteria */
 73    virtual bool shouldKeepTrack(const TrajectoryT& _t) = 0;
 74    bool shouldNotKeepTrack(const TrajectoryT& _t) { return !shouldKeepTrack(_t); }
 75
 76
 77   protected:
 78    std::string name;
 79};
 80
 81namespace bpo = boost::program_options;
 82
 83/** @brief Provides a basic min max filter and facilitates command line options to set the parameters
 84 *
 85 * The following filters on length, adding command line options "--min-length" & "--max-length"
 86 * @code
 87 * boost::program_options::options_description commandLineOptions;
 88 * MinMaxTrackFilter<double> lengthFilter("length",tracktable::length<TrajectoryT>);
 89 * lengthFilter.addOptions(commandLineOptions;)
 90 * //parse and notify
 91 * lengthFilter(trajectories);
 92 * @endcode
 93 *
 94 * @tparam MeasurementT The type of value we are expecting from the measurement function
 95 */
 96template <typename MeasurementT>
 97class MinMaxTrackFilter : public TrackFilter {
 98   protected:
 99    using ThisType = MinMaxTrackFilter<MeasurementT>;
100    using BaseType = TrackFilter;
101
102   public:
103    /**
104     * @brief Construct a new Min Max Track Filter object
105     *
106     * @param _name Used to create command line arguments for min/max
107     * @param _measureFunc Function that returns a value to be compared against min/max
108     *
109     * @note Can us inverseFilterTracks to exclude the range
110     */
111    MinMaxTrackFilter(const std::string& _name, MeasurementT (*_measureFunc)(const TrajectoryT&))
112        : BaseType(_name), measureFunc(_measureFunc) {}
113    bool shouldKeepTrack(const TrajectoryT& _t) {
114        auto r = measureFunc(_t);
115        return max >= r && r >= min;
116    }
117    void addOptions(bpo::options_description& _desc) {
118        auto minString = "min-" + name;
119        auto maxString = "max-" + name;
120        bpo::options_description filterOptions(name);
121        // clang-format off
122        filterOptions.add_options()
123            (minString.c_str(), bpo::value<MeasurementT>(&min),
124             ("minimum value for " + name).c_str())
125            (maxString.c_str(), bpo::value<MeasurementT>(&max),
126             ("maximum value for " + name).c_str())
127        ;
128        // clang-format on
129        _desc.add(filterOptions);
130    }
131
132    MeasurementT getMin() const { return min; }
133    void setMin(MeasurementT _v) { min = _v; }
134    MeasurementT getMax() const { return max; }
135    void setMax(MeasurementT _v) { max = _v; }
136
137   private:
138    MeasurementT min = std::numeric_limits<MeasurementT>::min();
139    MeasurementT max = std::numeric_limits<MeasurementT>::max();
140    MeasurementT (*measureFunc)(const TrajectoryT&) = nullptr;
141};
142
143#endif  // TrackFilter_h