У меня есть большое количество объектов, определенных координатами WGS84, которые я хотел бы преобразовать в пространство Меркатора. Есть другие операции, не показанные в MWE, которые мешают мне делать все однопоточным образом.
Мой текущий код выглядит следующим образом:
//Compile with: g++ -O3 temp.cpp -lproj -fopenmp
#include <proj_api.h>
#include <cstdlib>
#include <vector>
#include <utility>
#include <cmath>
typedef std::vector<double> dvec;
typedef std::pair<dvec,dvec> dvpair;
//Define these here so we're not making them every time we run `pj_transform()`
const char* wgs84_str = "+init=epsg:4326";
const char* merc_str = "+init=epsg:3857";
const projPJ pj_wgs84 = pj_init_plus(wgs84_str);
const projPJ pj_merc = pj_init_plus(merc_str);
//Use the global projections to do projecty stuff
void ToMercator(std::vector<double> &x, std::vector<double> &y){
pj_transform(pj_wgs84, pj_merc, x.size(), 1, x.data(), y.data(), NULL);
}
int main(int argc, char **argv){
std::vector<dvpair> dvs(1000);
for(int i=0;i<dvs.size();i++){
dvs[i].first.push_back(2*M_PI*rand()-M_PI); //Longitude in [-180,180]
dvs[i].second.push_back(M_PI*rand()-M_PI/2); //Latitude in [-90,90]
}
#pragma omp parallel for
for(int i=0;i<dvs.size();i++)
ToMercator(dvs[i].first,dvs[i].second);
}
Но это происходит с гейзенбаг. Получение дампов ядра с ulimit -c unlimited
показывает, что pj_transform
поднимает проблему. Это не слишком загадочно, так как некоторые раскопки показали, что proj4 не является потокобезопасным.
Способ исправить это — создать контексты для proj4 для каждого потока. за этот сайт. И это подводит нас к решению проблемы.
Каждый поток должен быть запущен pj_ctx_alloc()
который будет возвращать переменную типа потока projCtx
, который затем должен быть использован с pj_init_plus_ctx()
сделать projPJ
объекты, отвечающие за преобразование координат.
Как я могу создать глобальные потоки в OpenMP?
(Возможно, вы захотите ответить, что «глобальные перемены плохие». Я понимаю это, но учтите, что их удаление не меняет сути вопроса. У меня также есть отвращение к «бродячие данные«В моем реальном приложении есть как минимум шесть уровней промежуточных функций, прежде чем я достигну pj_transform()
У меня есть предпочтение не передавать параметры, которые глубоко в этом случае. Вы можете не согласиться с выбором дизайна, но это не подходящее место, чтобы спорить об этом.)
редактировать Это также не сработало:
#include <proj_api.h>
#include <cstdlib>
#include <vector>
#include <utility>
#include <cmath>
typedef std::vector<double> dvec;
typedef std::pair<dvec,dvec> dvpair;
//Define these here so we're not making them every time we run `pj_transform()`
const char* wgs84_str = "+init=epsg:4326";
const char* merc_str = "+init=epsg:3857";
projPJ pj_wgs84;
projPJ pj_merc;
projCtx pj_ctx;
//Use the global projections to do projecty stuff
void ToMercator(std::vector<double> &x, std::vector<double> &y){
pj_transform(pj_wgs84, pj_merc, x.size(), 1, x.data(), y.data(), NULL);
}
int main(int argc, char **argv){
std::vector<dvpair> dvs(1000);
for(int i=0;i<dvs.size();i++){
dvs[i].first.push_back(2*M_PI*rand()-M_PI); //Longitude in [-180,180]
dvs[i].second.push_back(M_PI*rand()-M_PI/2); //Latitude in [-90,90]
}
#pragma omp parallel private(pj_ctx, pj_wgs84, pj_merc)
{
pj_ctx = pj_ctx_alloc();
pj_wgs84 = pj_init_plus_ctx(pj_ctx,wgs84_str);
pj_merc = pj_init_plus_ctx(pj_ctx,merc_str);
#pragma omp for
for(int i=0;i<dvs.size();i++)
ToMercator(dvs[i].first,dvs[i].second);
}
}
редактировать
С помощью:
projPJ pj_wgs84;
projPJ pj_merc;
projCtx pj_ctx;
#pragma omp threadprivate(pj_wgs84)
#pragma omp threadprivate(pj_merc)
#pragma omp threadprivate(pj_ctx)
...
#pragma omp parallel private(pj_ctx, pj_wgs84, pj_merc)
{
pj_ctx = pj_ctx_alloc();
...
Выдает ошибку компоновщика:
/usr/bin/ld: pj_merc: TLS definition in /tmp/ccqZRl8N.o section .tbss mismatches non-TLS definition in /usr/lib/gcc/x86_64-linux-gnu/6/../../../x86_64-linux-gnu/libproj.so section .text
/usr/lib/gcc/x86_64-linux-gnu/6/../../../x86_64-linux-gnu/libproj.so: error adding symbols: Bad value
collect2: error: ld returned 1 exit status
Следующие работы:
#include <proj_api.h>
#include <cstdlib>
#include <vector>
#include <utility>
#include <cmath>
typedef std::vector<double> dvec;
typedef std::pair<dvec,dvec> dvpair;
//Define these here so we're not making them every time we run `pj_transform()`
const char* wgs84_str = "+init=epsg:4326";
const char* merc_str = "+init=epsg:3857";
//Using static means that these variables are file-scope. If we didn't use
//static then the variables might be accessed across many files and OpenMP would
//not like that
static projPJ pj_wgs84;
static projPJ pj_merc;
static projCtx pj_ctx;
//Set it so each thread has its own, private copy
#pragma omp threadprivate(pj_wgs84)
#pragma omp threadprivate(pj_merc)
#pragma omp threadprivate(pj_ctx)
//Use the global projections to do projecty stuff
void ToMercator(std::vector<double> &x, std::vector<double> &y){
pj_transform(pj_wgs84, pj_merc, x.size(), 1, x.data(), y.data(), NULL);
}
int main(int argc, char **argv){
std::vector<dvpair> dvs(10000);
for(int i=0;i<dvs.size();i++){
dvs[i].first.push_back(2*M_PI*rand()-M_PI); //Longitude in [-180,180]
dvs[i].second.push_back(M_PI*rand()-M_PI/2); //Latitude in [-90,90]
}
//Do this somewhere convenient ... like here
#pragma omp parallel
{
pj_ctx = pj_ctx_alloc();
pj_wgs84 = pj_init_plus_ctx(pj_ctx,wgs84_str);
pj_merc = pj_init_plus_ctx(pj_ctx,merc_str);
}
#pragma omp parallel for
for(int i=0;i<dvs.size();i++)
ToMercator(dvs[i].first,dvs[i].second);
}
Грустные новости: Тестирование этого на машине с большим количеством ядер (24 против 4) по-прежнему приводило к ошибкам сегмента.
Других решений пока нет …