This repository has been archived on 2020-04-08. You can view files and clone it, but cannot push or open issues or pull requests.
Files
FtmPrologic/code/main.cpp

524 lines
18 KiB
C++

#include "main.h"
#include "mesh.h"
#include "filter.h"
#include "Settings.h"
#include "meshPlotter.h"
#include "Plotty.h"
#include "Plotta.h"
#include <array>
#include <memory>
#include <thread>
#include <filesystem>
#include <chrono>
#include <Indoor/floorplan/v2/FloorplanReader.h>
#include <Indoor/sensors/offline/FileReader.h>
#include <Indoor/geo/Heading.h>
#include <Indoor/geo/Point2.h>
#include <Indoor/sensors/imu/TurnDetection.h>
#include <Indoor/sensors/imu/StepDetection.h>
#include <Indoor/sensors/imu/PoseDetection.h>
#include <Indoor/sensors/imu/MotionDetection.h>
#include <Indoor/sensors/pressure/RelativePressure.h>
#include <Indoor/data/Timestamp.h>
#include <Indoor/math/stats/Statistics.h>
#include "FtmKalman.h"
#include "mainFtm.h"
#include "misc.h"
#include <sys/stat.h>
using namespace std::chrono_literals;
std::vector<std::tuple<float, float, float>> getFtmValues(Offline::FileReader& fr, Interpolator<uint64_t, Point3>& gtInterpolator, const MACAddress nuc)
{
std::vector<std::tuple<float, float, float>> result;
for (const Offline::Entry& e : fr.getEntries())
{
if (e.type == Offline::Sensor::WIFI_FTM)
{
const Timestamp ts = Timestamp::fromMS(e.ts);
Point3 gtPos = gtInterpolator.get(static_cast<uint64_t>(ts.ms())) + Point3(0, 0, 1.3);
auto wifi = fr.getWifiFtm()[e.idx].data;
if (wifi.getAP().getMAC() == nuc)
{
Point3 apPos = Settings::data.CurrentPath.NUCs.find(wifi.getAP().getMAC())->second.position;
float apDist = gtPos.getDistance(apPos);
float ftmDist = wifi.getFtmDist();
float rssi = wifi.getRSSI();
result.push_back({ apDist, ftmDist, rssi });
}
}
}
return result;
}
std::pair<float, float> optimizeFtm(std::vector<std::tuple<float, float, float>>& values)
{
std::vector<std::pair<float, float>> error;
for (float offset = 0; offset < 10.0f; offset += 0.25)
{
Stats::Statistics<float> diffs;
for (const auto& data : values)
{
float apDist = std::get<0>(data);
float ftmDist = std::get<1>(data);
ftmDist += offset;
float diff = (apDist - ftmDist);
diffs.add(diff);
}
error.push_back({ offset, diffs.getSquaredSumAvg() });
}
auto minElement = std::min_element(error.begin(), error.end(), [](std::pair<float, float> a, std::pair<float, float> b) {
return a.second < b.second;
});
std::cout << "Min ftm offset \t" << minElement->first << "\t" << minElement->second << "\n";
return *minElement;
}
std::pair<float, float> optimizeRssi(std::vector<std::tuple<float, float, float>>& values)
{
std::vector<std::pair<float, float>> error;
for (float pathLoss = 2.0f; pathLoss < 4.0f; pathLoss += 0.125)
{
Stats::Statistics<float> diffs;
for (const auto& data : values)
{
float apDist = std::get<0>(data);
float rssi = std::get<2>(data);
float rssiDist = LogDistanceModel::rssiToDistance(-40, pathLoss, rssi);
float diff = (apDist - rssiDist);
diffs.add(diff);
}
error.push_back({ pathLoss, diffs.getSquaredSumAvg() });
}
auto minElement = std::min_element(error.begin(), error.end(), [](std::pair<float, float> a, std::pair<float, float> b) {
return a.second < b.second;
});
std::cout << "Min path loss \t" << minElement->first << "\t" << minElement->second << "\n";
return *minElement;
}
void optimizeWifiParameters(Offline::FileReader& fr, Interpolator<uint64_t, Point3>& gtInterpolator)
{
int i = 1;
for (auto nuc : { Settings::NUC1, Settings::NUC2, Settings::NUC3, Settings::NUC4 })
{
auto values = getFtmValues(fr, gtInterpolator, nuc);
std::cout << "NUC" << i++ << "\n";
optimizeFtm(values);
optimizeRssi(values);
}
}
void exportFtmValues(Offline::FileReader& fr, Interpolator<uint64_t, Point3>& gtInterpolator)
{
std::fstream fs;
fs.open("test.txt", std::fstream::out);
fs << "timestamp;nucid;dist;rssiDist;ftmDist;ftmStdDev;numMeas;numSuccesMeas" << "\n";
for (const Offline::Entry& e : fr.getEntries())
{
if (e.type == Offline::Sensor::WIFI_FTM)
{
const Timestamp ts = Timestamp::fromMS(e.ts);
Point3 gtPos = gtInterpolator.get(static_cast<uint64_t>(ts.ms())) + Point3(0, 0, 1.3);
auto wifi = fr.getWifiFtm()[e.idx].data;
int nucid = Settings::data.CurrentPath.NUCs.at(wifi.getAP().getMAC()).ID;
float ftm_offset = Settings::data.CurrentPath.NUCs.at(wifi.getAP().getMAC()).ftm_offset;
float rssi_pathloss = Settings::data.CurrentPath.NUCs.at(wifi.getAP().getMAC()).rssi_pathloss;
float rssiDist = LogDistanceModel::rssiToDistance(-40, rssi_pathloss, wifi.getRSSI());
float ftmDist = wifi.getFtmDist() + ftm_offset; //in m; plus static offset
float ftmStdDev = wifi.getFtmDistStd();
int numMeas = wifi.getNumAttemptedMeasurements();
int numSuccessMeas = wifi.getNumSuccessfulMeasurements();
Point3 apPos = Settings::data.CurrentPath.NUCs.find(wifi.getAP().getMAC())->second.position;
float apDist = gtPos.getDistance(apPos);
fs << ts.ms() << ";" << nucid << ";" << apDist << ";" << rssiDist << ";" << ftmDist << ";" << ftmStdDev << ";" << numMeas << ";" << numSuccessMeas << "\n";
}
}
fs.close();
}
static float kalman_procNoiseDistStdDev = 1.2f; // standard deviation of distance for process noise
static float kalman_procNoiseVelStdDev = 0.1f; // standard deviation of velocity for process noise
static CombinedStats<float> run(Settings::DataSetup setup, int walkIdx, std::string folder) {
// reading file
std::string currDir = std::filesystem::current_path().string();
Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(setup.map);
Offline::FileReader fr(setup.training[walkIdx], setup.HasNanoSecondTimestamps);
// ground truth
std::vector<int> gtPath = setup.gtPath;
Interpolator<uint64_t, Point3> gtInterpolator = fr.getGroundTruthPath(map, gtPath);
CombinedStats<float> errorStats;
//calculate distance of path
std::vector<Interpolator<uint64_t, Point3>::InterpolatorEntry> gtEntries = gtInterpolator.getEntries();
double distance = 0;
for(int i = 1; i < gtEntries.size(); ++i){
distance += gtEntries[i].value.getDistance(gtEntries[i-1].value);
}
std::cout << "Distance of Path: " << distance << std::endl;
// error file
const long int t = static_cast<long int>(time(NULL));
auto evalDir = std::filesystem::path(Settings::errorDir);
evalDir.append(folder);
if (!std::filesystem::exists(evalDir)) {
std::filesystem::create_directory(evalDir);
}
std::ofstream errorFile;
errorFile.open (evalDir.string() + "/" + std::to_string(walkIdx) + "_" + std::to_string(t) + ".csv");
// wifi
auto kalmanMap = std::make_shared<std::unordered_map<MACAddress, Kalman>>();
kalmanMap->insert({ Settings::NUC1, Kalman(1, setup.NUCs.at(Settings::NUC1).kalman_measStdDev, kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev) });
kalmanMap->insert({ Settings::NUC2, Kalman(2, setup.NUCs.at(Settings::NUC2).kalman_measStdDev, kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev) });
kalmanMap->insert({ Settings::NUC3, Kalman(3, setup.NUCs.at(Settings::NUC3).kalman_measStdDev, kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev) });
kalmanMap->insert({ Settings::NUC4, Kalman(4, setup.NUCs.at(Settings::NUC4).kalman_measStdDev, kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev) });
std::cout << "Optimal wifi parameters for " << setup.training[walkIdx] << "\n";
optimizeWifiParameters(fr, gtInterpolator);
// mesh
NM::NavMeshSettings set;
set.maxQuality_m = 0.10; // because of narrow hallways and small rooms reduce min. triangle size (default is 0.2)
MyNavMesh mesh;
MyNavMeshFactory fac(&mesh, set);
fac.build(map);
// debug show
//MeshPlotter dbg;
//dbg.addFloors(map);
//dbg.addOutline(map);
//dbg.addMesh(mesh);
////dbg.addDijkstra(mesh);
//dbg.draw();
Plotty plot(map);
plot.buildFloorplan();
plot.setGroundTruth(gtPath);
plot.setView(30, 0);
// APs Positions
plot.addCircle(100000 + 0, Settings::data.CurrentPath.nucInfo(0).position.xy(), 0.05);
plot.addCircle(100000 + 1, Settings::data.CurrentPath.nucInfo(1).position.xy(), 0.05);
plot.addCircle(100000 + 2, Settings::data.CurrentPath.nucInfo(2).position.xy(), 0.05);
plot.addCircle(100000 + 3, Settings::data.CurrentPath.nucInfo(3).position.xy(), 0.05);
plot.plot();
// particle-filter
const int numParticles = 5000;
//auto init = std::make_unique<MyPFInitFixed>(&mesh, srcPath0); // known position
auto init = std::make_unique<MyPFInitUniform>(&mesh); // uniform distribution
auto eval = std::make_unique<MyPFEval>();
eval->ftmKalmanFilters = kalmanMap;
auto trans = std::make_unique<MyPFTransRandom>();
//auto trans = std::make_unique<MyPFTransStatic>();
auto resample = std::make_unique<SMC::ParticleFilterResamplingSimple<MyState>>();
auto estimate = std::make_unique<SMC::ParticleFilterEstimationWeightedAverage<MyState>>();
//auto estimate = std::make_unique<SMC::ParticleFilterEstimationMax<MyState>>();
// setup
MyFilter pf(numParticles, std::move(init));
pf.setEvaluation(std::move(eval));
pf.setTransition(std::move(trans));
pf.setResampling(std::move(resample));
pf.setEstimation(std::move(estimate));
pf.setNEffThreshold(0.85);
// sensors
MyControl ctrl;
MyObservation obs;
Timestamp lastTimestamp = Timestamp::fromMS(0);
std::vector<float> errorValuesFtm, errorValuesRssi;
std::vector<int> timestamps;
std::vector<std::array<float, 4>> gtDistances, ftmDistances, rssiDistances; // distance per AP
Plotta::Plotta errorPlot("errorPlot", Settings::plotDataDir + "errorData.py");
Plotta::Plotta distsPlot("distsPlot", Settings::plotDataDir + "distances.py");
for (const Offline::Entry& e : fr.getEntries())
{
if (e.type != Offline::Sensor::WIFI_FTM) {
continue;
}
// TIME
const Timestamp ts = Timestamp::fromMS(e.ts);
auto wifiFtm = fr.getWifiFtm()[e.idx].data;
obs.ftm.push_back(wifiFtm);
if (ts - lastTimestamp >= Timestamp::fromMS(500))
{
// Do update step
Point2 gtPos = gtInterpolator.get(static_cast<uint64_t>(ts.ms())).xy();
plot.setGroundTruth(Point3(gtPos.x, gtPos.y, 0.1));
gtDistances.push_back({ gtPos.getDistance(Settings::data.CurrentPath.nucInfo(0).position.xy()),
gtPos.getDistance(Settings::data.CurrentPath.nucInfo(1).position.xy()),
gtPos.getDistance(Settings::data.CurrentPath.nucInfo(2).position.xy()),
gtPos.getDistance(Settings::data.CurrentPath.nucInfo(3).position.xy()) });
Point3 estPos;
float distErrorFtm = 0;
float distErrorRssi = 0;
// Run PF
obs.currentTime = ts;
ctrl.currentTime = ts;
MyState est = pf.update(&ctrl, obs);
ctrl.afterEval();
lastTimestamp = ts;
estPos = est.pos.pos;
ctrl.lastEstimate = estPos;
plot.setCurEst(Point3(estPos.x, estPos.y, 0.1));
plot.addEstimationNode(Point3(estPos.x, estPos.y, 0.1));
// Error
distErrorFtm = gtPos.getDistance(estPos.xy());
errorStats.ftm.add(distErrorFtm);
// draw wifi ranges
plot.clearDistanceCircles();
for (size_t i = 0; i < obs.ftm.size(); i++)
{
WiFiMeasurement wifi2 = obs.ftm[i];
Point3 apPos = Settings::data.CurrentPath.nuc(wifi2.getAP().getMAC()).position;
K::GnuplotColor color;
switch (Settings::data.CurrentPath.nuc(wifi2.getAP().getMAC()).ID)
{
case 1: color = K::GnuplotColor::fromRGB(0, 255, 0); break;
case 2: color = K::GnuplotColor::fromRGB(0, 0, 255); break;
case 3: color = K::GnuplotColor::fromRGB(255, 255, 0); break;
default: color = K::GnuplotColor::fromRGB(255, 0, 0); break;
}
plot.addDistanceCircle(apPos.xy(), wifi2.getFtmDist(), color);
}
obs.wifi.clear();
obs.ftm.clear();
errorValuesFtm.push_back(distErrorFtm);
errorValuesRssi.push_back(distErrorRssi);
timestamps.push_back(ts.ms());
// Error plot
errorPlot.add("t", timestamps);
errorPlot.add("errorFtm", errorValuesFtm);
errorPlot.add("errorRssi", errorValuesRssi);
errorPlot.frame();
// Distances plot
//distsPlot.add("t", timestamps);
//distsPlot.add("gtDists", gtDistances);
//distsPlot.add("ftmDists", ftmDistances);
//distsPlot.frame();
// Plotting
plot.showParticles(pf.getParticles());
plot.setCurEst(estPos);
plot.setGroundTruth(Point3(gtPos.x, gtPos.y, 0.1));
plot.addEstimationNode(estPos);
//plot.setActivity((int)act.get());
//plot.splot.getView().setEnabled(false);
//plot.splot.getView().setCamera(0, 0);
//plot.splot.getView().setEqualXY(true);
plot.plot();
}
}
printErrorStats(errorStats);
//system("pause");
return errorStats;
}
int main(int argc, char** argv)
{
CmdArguments args(argc, argv);
if (args.hasFlag("prob"))
{
std::cout << "Probabilistic" << "\n";
return mainProp(argc, argv);
}
else if (args.hasFlag("trilat"))
{
std::cout << "Trilateration" << "\n";
return mainTrilat(argc, argv);
}
CombinedStats<float> statsAVG;
CombinedStats<float> statsMedian;
CombinedStats<float> statsSTD;
CombinedStats<float> statsQuantil;
CombinedStats<float> tmp;
std::string evaluationName = "prologic/tmp";
for (size_t walkIdx = 0; walkIdx < 1 /*Settings::data.CurrentPath.training.size()*/; walkIdx++)
{
std::cout << "Executing walk " << walkIdx << "\n";
for (int i = 0; i < 5; ++i)
{
std::cout << "Start of iteration " << i << "\n";
tmp = run(Settings::data.CurrentPath, walkIdx, evaluationName);
statsAVG.ftm.add(tmp.ftm.getAvg());
statsMedian.ftm.add(tmp.ftm.getMedian());
statsSTD.ftm.add(tmp.ftm.getStdDev());
statsQuantil.ftm.add(tmp.ftm.getQuantile(0.75));
statsAVG.rssi.add(tmp.rssi.getAvg());
statsMedian.rssi.add(tmp.rssi.getMedian());
statsSTD.rssi.add(tmp.rssi.getStdDev());
statsQuantil.rssi.add(tmp.rssi.getQuantile(0.75));
std::cout << "Iteration " << i << " completed" << std::endl;
}
}
std::cout << "==========================================================" << std::endl;
std::cout << "Average of all statistical data FTM: " << std::endl;
std::cout << "Median: " << statsMedian.ftm.getAvg() << std::endl;
std::cout << "Average: " << statsAVG.ftm.getAvg() << std::endl;
std::cout << "Standard Deviation: " << statsSTD.ftm.getAvg() << std::endl;
std::cout << "75 Quantil: " << statsQuantil.ftm.getAvg() << std::endl;
std::cout << "==========================================================" << std::endl;
std::cout << "==========================================================" << std::endl;
std::cout << "Average of all statistical data RSSI: " << std::endl;
std::cout << "Median: " << statsMedian.rssi.getAvg() << std::endl;
std::cout << "Average: " << statsAVG.rssi.getAvg() << std::endl;
std::cout << "Standard Deviation: " << statsSTD.rssi.getAvg() << std::endl;
std::cout << "75 Quantil: " << statsQuantil.rssi.getAvg() << std::endl;
std::cout << "==========================================================" << std::endl;
//std::vector<std::array<float, 3>> error;
//std::ofstream error_out;
//error_out.open(Settings::errorDir + evaluationName + "_error_path1" + ".csv", std::ios_base::app);
//for (kalman_procNoiseDistStdDev = 0.8f; kalman_procNoiseDistStdDev < 1.5f; kalman_procNoiseDistStdDev += 0.1f)
//{
// for (kalman_procNoiseVelStdDev = 0.1f; kalman_procNoiseVelStdDev < 0.5f; kalman_procNoiseVelStdDev += 0.1f)
// {
//for (size_t walkIdx = 0; walkIdx < Settings::data.CurrentPath.training.size(); walkIdx++)
//{
// std::cout << "Executing walk " << walkIdx << "\n";
// for (int i = 0; i < 1; ++i)
// {
// std::cout << "Start of iteration " << i << "\n";
// tmp = run(Settings::data.CurrentPath, walkIdx, evaluationName);
// statsMedian.add(tmp.getMedian());
// statsAVG.add(tmp.getAvg());
// statsSTD.add(tmp.getStdDev());
// statsQuantil.add(tmp.getQuantile(0.75));
// std::cout << kalman_procNoiseDistStdDev << " " << kalman_procNoiseVelStdDev << std::endl;
// std::cout << "Iteration " << i << " completed" << std::endl;
// }
//}
// error.push_back({{ kalman_procNoiseDistStdDev, kalman_procNoiseVelStdDev, statsAVG.getAvg() }});
// auto minElement = std::min_element(error.begin(), error.end(), [](std::array<float, 3> a, std::array<float, 3> b) {
// return a[2] < b[2];
// });
// std::cout << "Current min error " << (*minElement)[2] << "\t Q(0)=\t" << (*minElement)[0] << "\t Q(1)=" << (*minElement)[1] << "\n";
// error_out << kalman_procNoiseDistStdDev << ";" << kalman_procNoiseVelStdDev << ";" << statsAVG.getAvg() << std::endl;
// // reset stats
// statsAVG.reset();
// statsMedian.reset();
// statsSTD.reset();
// statsQuantil.reset();
// }
//}
//auto minElement = std::min_element(error.begin(), error.end(), [](std::array<float, 3> a, std::array<float, 3> b) {
// return a[2] < b[2];
//});
//std::cout << "Global Min error " << (*minElement)[2] << "\t Q(0)=\t" << (*minElement)[0] << "\t Q(1)=" << (*minElement)[1] << "\n";
//error_out.close();
//std::this_thread::sleep_for(std::chrono::seconds(60));
}