FHWS/museumLoc
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

version for plotting the figures of kde paper

Browse Source
This commit is contained in:
toni
2018-02-19 18:18:45 +01:00
parent 46abeae148
commit 97ce18f5f1
12 changed files with 1374 additions and 1576 deletions
Download Patch File Download Diff File Expand all files Collapse all files

72
navMesh/filter.h
View File

@@ -3,6 +3,7 @@
#include "mesh.h"
#include "../Settings.h"
#include <omp.h>
#include <Indoor/geo/Heading.h>
#include <Indoor/math/Distributions.h>
@@ -12,15 +13,15 @@
#include <Indoor/smc/filtering/ParticleFilterInitializer.h>
#include <Indoor/smc/filtering/resampling/ParticleFilterResamplingSimple.h>
#include <Indoor/smc/filtering/resampling/ParticleFilterResamplingKLD.h>
#include <Indoor/smc/filtering/estimation/ParticleFilterEstimationBoxKDE.h>
#include <Indoor/smc/filtering/estimation/ParticleFilterEstimationWeightedAverage.h>
#include <Indoor/smc/filtering/estimation/ParticleFilterEstimationMax.h>
#include <Indoor/navMesh/walk/NavMeshWalkSimple.h>
#include <Indoor/navMesh/walk/NavMeshWalkEval.h>
#include <Indoor/sensors/radio/WiFiMeasurements.h>
#include <Indoor/data/Timestamp.h>
#include <Indoor/sensors/radio/WiFiProbabilityFree.h>
#include <Indoor/sensors/activity/ActivityDetector.h>
struct MyState {
@@ -32,6 +33,8 @@ struct MyState {
MyState() : pos(), heading(0) {;}
MyState(Point3 p) : pos(p, nullptr), heading(0){;}
MyState& operator += (const MyState& o) {
pos.tria = nullptr; // impossible
pos.pos += o.pos.pos;
@@ -50,6 +53,14 @@ struct MyState {
return res;
}
float getX(){
return pos.pos.x;
}
float getY() {
return pos.pos.y;
}
float getBinValue(const int dim) const {
switch (dim) {
case 0: return this->pos.pos.x;
@@ -87,6 +98,9 @@ struct MyObservation {
//time
Timestamp currentTime;
//activity
Activity activity;
};
class MyPFInitUniform : public SMC::ParticleFilterInitializer<MyState> {
@@ -156,21 +170,28 @@ public:
void transition(std::vector<SMC::Particle<MyState>>& particles, const MyControl* control) override {
Distribution::Normal<float> dStepSizeFloor(0.70, 0.1);
Distribution::Normal<float> dStepSizeFloor(0.70, 0.1);
Distribution::Normal<float> dStepSizeStair(0.35, 0.1);
Distribution::Normal<float> dHeading(0.0, 0.10);
Distribution::Normal<float> dHeading(0.0, 0.1);
for (SMC::Particle<MyState>& p : particles) {
#pragma omp parallel for num_threads(3)
for (int i = 0; i < particles.size(); ++i) {
SMC::Particle<MyState>& p = particles[i];
// how to walk
MyNavMeshWalkParams params;
params.heading = p.state.heading + control->headingChangeSinceLastEval + dHeading.draw();
params.numSteps = control->numStepsSinceLastEval;
params.start = p.state.pos;
params.stepSizes.stepSizeFloor_m = dStepSizeFloor.draw();
params.stepSizes.stepSizeStair_m = dStepSizeStair.draw();
if(params.stepSizes.stepSizeFloor_m < 0.1 || params.stepSizes.stepSizeStair_m < 0.1){
params.stepSizes.stepSizeFloor_m = 0.1;
params.stepSizes.stepSizeStair_m = 0.1;
}
// walk
MyNavMeshWalk::ResultEntry res = walker.getOne(params);
@@ -194,6 +215,29 @@ class MyPFEval : public SMC::ParticleFilterEvaluation<MyState, MyObservation> {
WiFiModel& wifiModel;
WiFiObserverFree wifiProbability;
//TODO: add this to transition probability
double getStairProb(const SMC::Particle<MyState>& p, const Activity act) {
const float kappa = 0.75;
switch (act) {
case Activity::WALKING:
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::FLOOR_INDOOR) {return kappa;}
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::DOOR) {return kappa;}
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::STAIR_LEVELED) {return kappa;}
{return 1-kappa;}
case Activity::WALKING_UP:
case Activity::WALKING_DOWN:
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::STAIR_SKEWED) {return kappa;}
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::STAIR_LEVELED) {return kappa;}
if (p.state.pos.tria->getType() == (int) NM::NavMeshType::ELEVATOR) {return kappa;}
{return 1-kappa;}
}
return 1.0;
}
public:
MyPFEval(WiFiModel& wifiModel) : wifiModel(wifiModel), wifiProbability(Settings::WiFiModel::sigma, wifiModel){}
@@ -203,13 +247,23 @@ public:
double sum = 0;
const WiFiMeasurements wifiObs = Settings::WiFiModel::vg_eval.group(observation.wifi);
for (SMC::Particle<MyState>& p : particles) {
#pragma omp parallel for num_threads(3)
for (int i = 0; i < particles.size(); ++i) {
SMC::Particle<MyState>& p = particles[i];
double pWifi = wifiProbability.getProbability(p.state.pos.pos, observation.currentTime, wifiObs);
double pStair = getStairProb(p, observation.activity);
const double prob = pWifi;
//HACK HACK HACK HACK
double prob = 1.0;
if(observation.currentTime.ms() > 20801){
prob = pWifi * pStair;
}
p.weight *= prob;
#pragma omp atomic
sum += prob;
}

183
navMesh/main.h
View File

@@ -1,183 +0,0 @@
#ifndef NAV_MESH_MAIN_H
#define NAV_MESH_MAIN_H
#include "mesh.h"
#include "filter.h"
#include "../Settings.h"
#include <memory>
#include <thread>
#include <Indoor/floorplan/v2/FloorplanReader.h>
#include <Indoor/sensors/radio/model/WiFiModelLogDistCeiling.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/MotionDetection.h>
#include <Indoor/sensors/pressure/RelativePressure.h>
#include <Indoor/data/Timestamp.h>
#include <Indoor/sensors/radio/setup/WiFiOptimizerLogDistCeiling.h>
void navMeshMain() {
//std::string mapFile = "/apps/paper/diss/data/maps/museum31.xml";
std::string mapFile = "../map/map42_ap.xml";
// reading file
Floorplan::IndoorMap* map = Floorplan::Reader::readFromFile(mapFile);
Offline::FileReader fr("../measurements/museum/Pixel/Path1_2468.csv");
WiFiFingerprints fingerprints("../measurements/museum/Nexus/fingerprints/wifi_fp.dat");
const std::string wifiModelFile = "../measurements/museum/wifimodel.dat";
std::ifstream inp(wifiModelFile, std::ifstream::binary);
// wifi
WiFiModelLogDistCeiling WiFiModel(map);
// with optimization
if(Settings::WiFiModel::optimize){
if (!inp.good() || (inp.peek()&&0) || inp.eof()) {
Assert::isFalse(fingerprints.getFingerprints().empty(), "no fingerprints available!");
WiFiOptimizer::LogDistCeiling opt(map, Settings::WiFiModel::vg_calib);
for (const WiFiFingerprint& fp : fingerprints.getFingerprints()) {
opt.addFingerprint(fp);
}
const WiFiOptimizer::LogDistCeiling::APParamsList res = opt.optimizeAll(opt.NONE);
for (const WiFiOptimizer::LogDistCeiling::APParamsMAC& ap : res.get()) {
const WiFiModelLogDistCeiling::APEntry entry(ap.params.getPos(), ap.params.txp, ap.params.exp, ap.params.waf);
WiFiModel.addAP(ap.mac, entry);
}
WiFiModel.saveXML(wifiModelFile);
} else {
WiFiModel.loadXML(wifiModelFile);
}
} else {
// without optimization
WiFiModel.loadAPs(map, Settings::WiFiModel::TXP, Settings::WiFiModel::EXP, Settings::WiFiModel::WAF);
Assert::isFalse(WiFiModel.getAllAPs().empty(), "no AccessPoints stored within the map.xml");
}
// mesh
NM::NavMeshSettings set;
MyNavMesh mesh;
MyNavMeshFactory fac(&mesh, set);
fac.build(map);
const Point3 srcPath0(26, 43, 7.5);
const Point3 srcPath1(62, 38, 1.8);
// add shortest-path to destination
//const Point3 dst(51, 45, 1.7);
//const Point3 dst(25, 45, 0);
//NM::NavMeshDijkstra::stamp<MyNavMeshTriangle>(mesh, dst);
// debug show
NM::NavMeshDebug dbg;
dbg.addMesh(mesh);
//dbg.addDijkstra(mesh);
dbg.draw();
// particle-filter
const int numParticles = 1000;
//auto init = std::make_unique<MyPFInitFixed>(&mesh, srcPath1); // known position
auto init = std::make_unique<MyPFInitUniform>(&mesh); // uniform distribution
auto eval = std::make_unique<MyPFEval>(WiFiModel);
auto trans = std::make_unique<MyPFTrans>(mesh);
auto resample = std::make_unique<SMC::ParticleFilterResamplingSimple<MyState>>();
auto estimate = std::make_unique<SMC::ParticleFilterEstimationWeightedAverage<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(1);
// sensors
MyControl ctrl;
MyObservation obs;
StepDetection sd;
PoseDetection pd;
TurnDetection td(&pd);
RelativePressure relBaro;
relBaro.setCalibrationTimeframe( Timestamp::fromMS(5000) );
Timestamp lastTimestamp = Timestamp::fromMS(0);
// parse each sensor-value within the offline data
for (const Offline::Entry& e : fr.getEntries()) {
const Timestamp ts = Timestamp::fromMS(e.ts);
if (e.type == Offline::Sensor::WIFI) {
obs.wifi = fr.getWiFiGroupedByTime()[e.idx].data;
} else if (e.type == Offline::Sensor::ACC) {
if (sd.add(ts, fr.getAccelerometer()[e.idx].data)) {
++ctrl.numStepsSinceLastEval;
}
const Offline::TS<AccelerometerData>& _acc = fr.getAccelerometer()[e.idx];
pd.addAccelerometer(ts, _acc.data);
} else if (e.type == Offline::Sensor::GYRO) {
const Offline::TS<GyroscopeData>& _gyr = fr.getGyroscope()[e.idx];
const float delta_gyro = td.addGyroscope(ts, _gyr.data);
ctrl.headingChangeSinceLastEval += delta_gyro;
} else if (e.type == Offline::Sensor::BARO) {
relBaro.add(ts, fr.getBarometer()[e.idx].data);
obs.relativePressure = relBaro.getPressureRealtiveToStart();
obs.sigmaPressure = relBaro.getSigma();
}
if (ts.ms() - lastTimestamp.ms() > 500 && ctrl.numStepsSinceLastEval > 0) {
obs.currentTime = ts;
// if(ctrl.numStepsSinceLastEval > 0){
// pf.updateTransitionOnly(&ctrl);
// ctrl.afterEval();
// }
// MyState est = pf.updateEvaluationOnly(obs);
// lastTimestamp = ts;
MyState est = pf.update(&ctrl, obs);
ctrl.afterEval();
lastTimestamp = ts;
// try {
// MyNavMeshLocation loc = mesh.getLocationNearestTo(est.pos.pos);
// auto path = loc.tria->getPathToDestination<MyNavMeshTriangle>(loc.pos);
// dbg.addDijkstra(path);
// } catch (...) {;}
// const int d = (i * 1) % 360;
// dbg.plot.getView().setCamera(60, d);
dbg.showParticles(pf.getParticles());
dbg.setCurPos(est.pos.pos);
//dbg.gp.setOutput("/tmp/123/" + std::to_string(i) + ".png");
//dbg.gp.setTerminal("pngcairo", K::GnuplotSize(60, 30));
// std::cout << i << std::endl;
dbg.draw();
std::this_thread::sleep_for(std::chrono::milliseconds(1));
}
}
}
#endif

242
navMesh/meshPlotter.h Normal file
View File

@@ -0,0 +1,242 @@
#ifndef MESHPLOTTER_H
#define MESHPLOTTER_H
#include <KLib/misc/gnuplot/Gnuplot.h>
#include <KLib/misc/gnuplot/GnuplotSplot.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementLines.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementPoints.h>
#include <KLib/misc/gnuplot/GnuplotSplotElementColorPoints.h>
#include <KLib/misc/gnuplot/objects/GnuplotObjectPolygon.h>
#include <Indoor/math/Distributions.h>
#include <Indoor/navMesh/NavMesh.h>
#include <Indoor/floorplan/v2/Floorplan.h>
class NavMeshTriangleDijkstra;
/**
* debug plot NavMeshes
*/
class MeshPlotter {
public:
K::Gnuplot gp;
K::GnuplotSplot plot;
K::GnuplotSplotElementLines pFloor;
K::GnuplotSplotElementLines pOutline;
K::GnuplotSplotElementLines lines;
K::GnuplotSplotElementPoints border;
K::GnuplotSplotElementColorPoints particles;
K::GnuplotSplotElementColorPoints distances;
K::GnuplotSplotElementLines pathEstimated;
K::GnuplotSplotElementLines shortestPath;
K::GnuplotSplotElementLines groundtruthPath;
private:
K::GnuplotFill gFill[6] = {
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#0000ff"), 1), // unknown
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#999999"), 1), // indoor
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#44ffee"), 1), // outdoor
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#666699"), 1), // door
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#444444"), 1), // stairs_level
K::GnuplotFill(K::GnuplotFillStyle::SOLID, K::GnuplotColor::fromHexStr("#666666"), 1) // stairs_skewed
};
public:
MeshPlotter() {
gp << "set view equal xy\n";
plot.add(&lines); lines.setShowPoints(true);
plot.add(&border);
plot.add(&particles); particles.setPointType(7); particles.setPointSize(0.2);
plot.add(&pathEstimated); pathEstimated.getStroke().setWidth(2); pathEstimated.setShowPoints(false); pathEstimated.getStroke().getColor().setHexStr("#00ff00");
plot.add(&groundtruthPath); groundtruthPath.getStroke().setWidth(2); groundtruthPath.getStroke().getColor().setHexStr("#000000");
plot.add(&distances); distances.setPointSize(0.75); distances.setPointType(7);
plot.add(&shortestPath); shortestPath.getStroke().setWidth(3);
plot.add(&pFloor);
plot.add(&pOutline); pOutline.getStroke().getColor().setHexStr("#999999");
}
void draw() {
gp.draw(plot);
gp.flush();
}
template <typename T> void showParticles(const std::vector<T>& particles) {
this->particles.clear();
double min = +999;
double max = -999;
for (const T& p : particles) {
const K::GnuplotPoint3 p3(p.state.pos.pos.x, p.state.pos.pos.y, p.state.pos.pos.z);
const double prob = std::pow(p.weight, 0.25);
this->particles.add(p3, prob);
if (prob > max) {max = prob;}
if (prob < min) {min = prob;}
}
plot.getAxisCB().setRange(min, max + 0.000001);
}
template <typename Tria> void addMesh(NM::NavMesh<Tria>& nm) {
K::GnuplotStroke gStroke = K::GnuplotStroke(K::GnuplotDashtype::SOLID, 1, K::GnuplotColor::fromHexStr("#666600"));
const BBox3 bbox = nm.getBBox();
border.add(K::GnuplotPoint3(bbox.getMin().x,bbox.getMin().y,bbox.getMin().z));
border.add(K::GnuplotPoint3(bbox.getMax().x,bbox.getMax().y,bbox.getMax().z));
// lines.add(K::GnuplotPoint3(bbox.getMin().x,bbox.getMin().y,bbox.getMin().z), K::GnuplotPoint3(bbox.getMax().x, 0, 0));
// lines.add(K::GnuplotPoint3(bbox.getMin().x,bbox.getMin().y,bbox.getMin().z), K::GnuplotPoint3(0,bbox.getMax().y,0));
// lines.addSegment(K::GnuplotPoint3(bbox.getMin().x,bbox.getMin().y,bbox.getMin().z), K::GnuplotPoint3(0,0,bbox.getMax().z));
//stairs in eigene group? vlt gehen dann auch die dellen weg?
for (const Tria* tria : nm) {
const uint8_t type = tria->getType();
if (type < 0 || type > 5) {
throw std::runtime_error("out of type-bounds");
}
K::GnuplotObjectPolygon* pol = new K::GnuplotObjectPolygon(gFill[type], gStroke);
pol->add(K::GnuplotCoordinate3(tria->getP1().x, tria->getP1().y, tria->getP1().z, K::GnuplotCoordinateSystem::FIRST));
pol->add(K::GnuplotCoordinate3(tria->getP2().x, tria->getP2().y, tria->getP2().z, K::GnuplotCoordinateSystem::FIRST));
pol->add(K::GnuplotCoordinate3(tria->getP3().x, tria->getP3().y, tria->getP3().z, K::GnuplotCoordinateSystem::FIRST));
pol->close();
pol->setZIndex(tria->getP3().z);
plot.getObjects().add(pol);
//for (int i = 0; i < nm.getNumNeighbors(tria); ++i) {
// const Tria* o = nm.getNeighbor(tria, i);
// const Point3 p1 = tria->getCenter();
// const Point3 p2 = o.getCenter();
// //lines.addSegment(K::GnuplotPoint3(p1.x,p1.y,p1.z+0.1), K::GnuplotPoint3(p2.x,p2.y,p2.z+0.1));
//}
for (const NM::NavMeshTriangle* o : *tria) {
const Point3 p1 = tria->getCenter();
const Point3 p2 = o->getCenter();
// lines.addSegment(K::GnuplotPoint3(p1.x,p1.y,p1.z+0.1), K::GnuplotPoint3(p2.x,p2.y,p2.z+0.1));
}
}
plot.getObjects().reOrderByZIndex();
}
template <typename Tria> void addDijkstra(NM::NavMesh<Tria>& mesh) {
distances.clear();
// ensure Tria extends NavMeshTriangleDijkstra
StaticAssert::AinheritsB<Tria, NavMeshTriangleDijkstra>();
NM::NavMeshRandom<Tria> rnd = mesh.getRandom();
for (int i = 0; i < 5000; ++i) {
NM::NavMeshLocation<Tria> loc = rnd.draw();
float v = loc.tria->interpolate(loc.pos, loc.tria->spFromP1.distance, loc.tria->spFromP2.distance, loc.tria->spFromP3.distance);
distances.add(K::GnuplotPoint3(loc.pos.x, loc.pos.y, loc.pos.z), v);
}
// Distribution::Uniform<float> dist (-0.5, +0.5);
// for (const Tria* t : mesh) {
// const Point3 posC = t->getCenter();
// distances.add(K::GnuplotPoint3(posC.x+dist.draw(), posC.y+dist.draw(), posC.z), t->distAtCenter);
// const Point3 pos1 = t->getP1();
// distances.add(K::GnuplotPoint3(pos1.x+dist.draw(), pos1.y+dist.draw(), pos1.z), t->distAtP1);
// const Point3 pos2 = t->getP2();
// distances.add(K::GnuplotPoint3(pos2.x+dist.draw(), pos2.y+dist.draw(), pos2.z), t->distAtP2);
// const Point3 pos3 = t->getP3();
// distances.add(K::GnuplotPoint3(pos3.x+dist.draw(), pos3.y+dist.draw(), pos3.z), t->distAtP3);
// }
}
template <typename Tria> void addDijkstra(std::vector<NM::NavMeshLocation<Tria>>& path) {
shortestPath.clear();
for (auto& e : path) {
K::GnuplotPoint3 gp(e.pos.x, e.pos.y, e.pos.z);
shortestPath.add(gp);
}
}
void addGroundTruthNode(const Point3 pos) {
K::GnuplotPoint3 gp(pos.x, pos.y, std::round(pos.z * 10) / 10);
groundtruthPath.add(gp);
}
void addEstimationNode(const Point3 pos){
K::GnuplotPoint3 est(pos.x, pos.y, std::round(pos.z * 10) / 10);
pathEstimated.add(est);
}
void setTimeInMinute(const int minutes, const int seconds) {
gp << "set label 1002 at screen 0.02, 0.94 'Time: " << minutes << ":" << seconds << "'\n";
}
void setGT(const Point3 pt) {
gp << "set arrow 31337 from " << pt.x << "," << pt.y << "," << (pt.z+1.4) << " to " << pt.x << "," << pt.y << "," << pt.z << " front \n";
}
void setCurPos(const Point3 pt) {
gp << "set arrow 31338 from " << pt.x << "," << pt.y << "," << (pt.z+0.9) << " to " << pt.x << "," << pt.y << "," << pt.z << " lw 2 lc 'green' front \n";
}
void saveToFile(std::ofstream& stream){
gp.draw(plot);
stream << "set terminal x11 size 2000,1500\n";
stream << gp.getBuffer();
stream << "pause -1\n";
gp.flush();
}
void printOverview(const std::string& path) {
gp << "set terminal png size 1280,720\n";
gp << "set output '" << path << "_overview" << ".png'\n";
gp << "set view 75,60\n";
gp << "set autoscale xy\n";
gp << "set autoscale z\n";
draw();
}
//meshless drawing
void addFloors(Floorplan::IndoorMap* map) {
for (Floorplan::Floor* f : map->floors) {
for (Floorplan::FloorObstacle* obs : f->obstacles) {
Floorplan::FloorObstacleLine* line = dynamic_cast<Floorplan::FloorObstacleLine*>(obs);
if (line) {
K::GnuplotPoint3 p1(line->from.x, line->from.y, f->atHeight);
K::GnuplotPoint3 p2(line->to.x, line->to.y, f->atHeight);
pFloor.addSegment(p1, p2);
}
}
}
}
void addOutline(Floorplan::IndoorMap* map) {
for (Floorplan::Floor* f : map->floors) {
for (Floorplan::FloorOutlinePolygon* poly : f->outline) {
const int cnt = poly->poly.points.size();
for (int i = 0; i < cnt; ++i) {
Point2 p1 = poly->poly.points[(i+0)];
Point2 p2 = poly->poly.points[(i+1)%cnt];
K::GnuplotPoint3 gp1(p1.x, p1.y, f->atHeight);
K::GnuplotPoint3 gp2(p2.x, p2.y, f->atHeight);
pOutline.addSegment(gp1, gp2);
}
}
}
}
};
#endif // MESHPLOTTER_H