Classify Example Source Files¶
Table of Contents
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