10 float nom[WIDTH * HEIGHT], denom[WIDTH * HEIGHT];
11 unsigned char pix[WIDTH * HEIGHT * 3];
20 bool operator< (const point &other) const
22 if (fabs(lat - other.lat) > EPS)
23 return (lat < other.lat);
24 if (fabs(lon - other.lon) > EPS)
25 return (lon < other.lon);
26 return (strength < other.strength);
29 bool same_place_as(const point &other) const
31 return !((fabs(lat - other.lat) > EPS) || (fabs(lon - other.lon) > EPS));
35 void collapse_points(std::vector<point> &points, std::vector<point> &upoints)
37 // average non-unique points
38 point cp = { 0.0, 0.0, 0 };
41 for (std::vector<point>::const_iterator i = points.begin(); i != points.end(); ++i) {
48 if (i->same_place_as(cp)) {
49 cp.strength += i->strength;
53 upoints.push_back(cp);
61 upoints.push_back(cp);
66 std::vector<point> points, upoints;
67 double max_lat = -HUGE_VAL, min_lat = HUGE_VAL;
68 double max_lon = -HUGE_VAL, min_lon = HUGE_VAL;
72 if (scanf("%lf %lf %lf", &p.lat, &p.lon, &p.strength) != 3)
75 max_lon = std::max(max_lon, p.lon);
76 min_lon = std::min(min_lon, p.lon);
78 max_lat = std::max(max_lat, p.lat);
79 min_lat = std::min(min_lat, p.lat);
84 double x_res = WIDTH / (max_lon - min_lon);
85 double y_res = HEIGHT / (max_lat - min_lat);
88 double lat_dist = max_lat - min_lat;
89 double lat_wanted = HEIGHT / x_res;
90 double extra_lat = lat_wanted - lat_dist;
92 min_lat -= extra_lat * 0.5;
93 max_lat += extra_lat * 0.5;
95 y_res = HEIGHT / (max_lat - min_lat);
97 fprintf(stderr, "FIXME\n");
100 std::sort(points.begin(), points.end());
101 collapse_points(points, upoints);
104 for (std::vector<point>::iterator i = upoints.begin(); i != upoints.end(); ++i) {
105 i->x = (int)((i->lon - min_lon) * WIDTH / (max_lon - min_lon));
106 i->y = (int)((i->lat - min_lat) * HEIGHT / (max_lat - min_lat));
107 i->strength = 0.02 * (i->strength + 80);
108 if (i->strength < 0.0)
113 for (std::vector<point>::iterator i = upoints.begin(); i != upoints.end(); ++i) {
115 for (std::vector<point>::iterator j = upoints.begin(); j != upoints.end(); ++j) {
117 (i->x - j->x) * (i->x - j->x) +
118 (i->y - j->y) * (i->y - j->y);
121 min_rad = std::min(min_rad, rad);
124 i->radius = std::min(sqrt(min_rad), 15.0);
128 std::fill(nom, nom + WIDTH*HEIGHT, 0.0);
129 std::fill(denom, denom + WIDTH*HEIGHT, 0.001);
131 for (std::vector<point>::iterator i = upoints.begin(); i != upoints.end(); ++i) {
132 double r = i->radius * 3.0;
133 double sigma = i->radius;
135 for (int y = -r; y <= r; ++y) {
136 for (int x = -r; x <= r; ++x) {
140 if (rx < 0 || ry < 0 || rx >= WIDTH || ry >= HEIGHT)
143 double z = - (x*x + y*y) / (2.0 * sigma * sigma);
144 double fac = exp(z) / sigma;
146 nom [ry * WIDTH + rx] += i->strength * fac;
147 denom[ry * WIDTH + rx] += fac;
153 for (unsigned y = 0; y < HEIGHT; ++y) {
154 for (unsigned x = 0; x < WIDTH; ++x) {
156 if (denom[y * WIDTH + x] < EPS)
159 ch = (unsigned)(255.0 * nom[y * WIDTH + x] / denom[y * WIDTH + x]);
164 pix[(y * WIDTH + x) * 3] = ch;
165 pix[(y * WIDTH + x) * 3 + 1] = ch;
166 pix[(y * WIDTH + x) * 3 + 2] = ch;
171 FILE *track = fopen("track.txt", "r");
174 if (fscanf(track, "%lf %lf", &lat, &lon) != 2)
177 int x = (int)((lon - min_lon) * WIDTH / (max_lon - min_lon));
178 int y = (int)((lat - min_lat) * HEIGHT / (max_lat - min_lat));
180 if (x < 0 || y < 0 || x >= WIDTH || y >= HEIGHT)
182 pix[(y * WIDTH + x) * 3] = 255;
185 // and the sampling points
186 for (std::vector<point>::iterator i = upoints.begin(); i != upoints.end(); ++i) {
187 int x = (int)((i->lon - min_lon) * WIDTH / (max_lon - min_lon));
188 int y = (int)((i->lat - min_lat) * HEIGHT / (max_lat - min_lat));
190 if (x < 0 || y < 0 || x >= WIDTH || y >= HEIGHT)
192 pix[(y * WIDTH + x) * 3 + 2] = 255;
196 printf("P6\n%u %u\n255\n", WIDTH, HEIGHT);
197 fwrite(pix, 3, WIDTH*HEIGHT, stdout);