几个月前我有类似的问题,但我只需要计算一维边缘。这是对我有用的解决方案的大纲,我想它也可以适应多维边缘:
我基本上是将pdf存储在一维数组/矢量中(使用任何你喜欢的):
double* pdf = new double[a*b*c*d*e];
然后我用您可以在二维阵列a[width][height]
存储为一维阵列b[widht*height]
和访问任何元素作为a[x][y]
b[width*x + y]
。您可以将此公式推广到任意维度,并且正确使用模/整数除法,您还可以计算反比。
使用模板,从一维索引到N维索引以及反之亦然的计算都非常简单。这使您可以将取决于维度的符号PDF[a][b][c][d][e]
转换为易于扩展到任意维度的符号,如PDF(std::vector<size_t>{a,b,c,d,e})
,因为您可以预先在循环中填充该向量。
如果您认为这种方法可能对您有所帮助,我可以尝试获取我的实施的一些关键功能并将其添加到此处。
编辑:
template <size_t DIM>
inline void posToPosN(const size_t& pos,
const size_t* const size,
size_t* const posN){
size_t r = pos;
for (size_t i = DIM; i > 0; --i){
posN[i - 1] = r % size[i - 1];
r /= size[i - 1];
}
}
template <size_t DIM>
inline void posNToPos(size_t& pos,
const size_t* const size,
const size_t* const posN){
pos = 0;
size_t mult = 1;
for (size_t i = DIM; i > 0; --i){
pos += mult * posN[i - 1];
mult *= size[i - 1];
}
}
template<typename type, size_t DIM>
class Iterator{
private:
type* const _data; //pointer to start of Array
size_t _pos; //1-dimensional position
size_t _posN[DIM]; //n-dimensional position
size_t const * const _size; //pointer to the _size-Member of Array
size_t _total;
private:
public:
Iterator(type* const data, const size_t* const size, size_t total, size_t pos)
: _data(data), _pos(pos), _size(size), _total(total)
{
if (_pos > _total || _pos < 0) _pos = _total;
posToPosN<DIM>(_pos, _size, _posN);
}
bool operator!= (const Iterator& other) const
{
return _pos != other._pos;
}
type& operator*() const{
if (_pos >= _total)
std::cout << "ERROR, dereferencing too high operator";
return *(_data + _pos);
}
const Iterator& operator++()
{
++_pos;
if (_pos > _total) _pos = _total;
posToPosN<DIM>(_pos, _size, _posN);
return *this;
}
Iterator& operator +=(const size_t& b)
{
_pos += b;
if (_pos > _total) _pos = _total;
posToPosN<DIM>(_pos, _size, _posN);
return *this;
}
const Iterator& operator--()
{
if (_pos == 0)
_pos = _total;
else
--_pos;
posToPosN<DIM>(_pos, _size, _posN);
return *this;
}
//returns position in n-th dimension
size_t operator[](size_t n){
return _posN[n];
}
//returns a new iterator, advanced by n steps in the dim Dimension
Iterator advance(size_t dim, int steps = 1){
if (_posN[dim] + steps < 0 || _posN[dim] + steps >= _size[dim]){
return Iterator(_data, _size, _total, _total);
}
size_t stride = 1;
for (size_t i = DIM - 1; i > dim; --i){
stride *= _size[i];
}
return Iterator(_data, _size, _total, _pos + steps*stride);
}
};
template <typename type, size_t DIM>
class Array{
type* _data;
size_t _size[DIM];
size_t _total;
void init(const size_t* const dimensions){
_total = 1;
for (int i = 0; i < DIM; i++){
_size[i] = dimensions[i];
_total *= _size[i];
}
_data = new type[_total];
}
public:
Array(const size_t* const dimensions){
init(dimensions);
}
Array(const std::array<size_t, DIM>& dimensions){
init(&dimensions[0]);
}
~Array(){
delete _data;
}
Iterator<type, DIM> begin(){
return Iterator<type, DIM>(_data, _size, _total, 0);
}
Iterator<type, DIM> end(){
return Iterator<type, DIM>(_data, _size, _total, _total);
}
const size_t* const size(){
return _size;
}
};
//for projections of the PDF
void calc_marginals(size_t dir, double* p_xPos, double* p_yPos){
assert(dir < N_THETA);
std::lock_guard<std::mutex> lock(calcInProgress);
//reset to 0
for (size_t i = 0; i < _size[dir]; ++i){
p_yPos[i] = 0;
}
//calc projection
double sum = 0;
for (auto it = _p_theta.begin(); it != _p_theta.end(); ++it){
p_yPos[it[dir]] += (*it);
sum += (*it);
}
if (abs(sum - 1) > 0.001){ cout << "Warning: marginal[" << dir << "] not normalized" << endl; }
//calc x-Axis
for (size_t i = 0; i < _size[dir]; ++i){
p_xPos[i] = _p[dir].start + double(i)/double(_size[dir] - 1)*(_p[dir].stop - _p[dir].start);
}
}
该代码由几部分组成:
- 两个函数
posToPosN()
和posNToPos()
其中做一维和DIM维坐标之间的变换提到。尺寸在这里作为模板参数DIM给出。 pos
只是一维位置,posN
一个指向DIM30维坐标的尺寸为的数组的指针,而size
是一个尺寸为的数组,其宽度在不同的方向上(在您的情况下,{a,b, c,d,e})
Iterator
是一个Iterator类,允许在DIM-Dimensional Array上使用基于范围或基于迭代器的for-loops。注意它返回DIM维坐标的第n个分量和advance()
函数返回一个迭代器具有坐标的元素{posN[0], posN[1], ...,posN[dim] + steps , ... posN[DIM]}
Array
应该是相当简单明了
calcMarginals
的operator[](size_t n)
为i用于计算边缘人的函数。 dir
这里是我想要计算边际的方向(记住:一维边际),并写入p_xPos
和p_yPos
,_p_theta
是Array
。注意基于迭代器的for循环,(*it)
这里指的是存储在数组中的pdf的double值,就像通常的迭代器一样。另外,it[dim]
返回方向dim中实际值的坐标。写入p_xPos的最后一个循环仅仅是因为我不想在这个数组中使用索引而是实际值。
我想如果你重新定义Iterator::operator[]
采取维度索引的矢量/阵列并返回相应坐标的矢量/阵列并添加Array::operator[]
随机访问,这需要矢量/阵列,以及你应该几乎完成了。
谢谢您的回复!我认为这种方法可能会有所帮助。你能提供这些关键功能吗? – Rene
@Rene将我的代码添加到它,希望它有帮助 – Anedar