基于Permutohedral Lattice 的Bilateral filter 源码及部分注释【C++】
【摘要】
基于Permutohedral Lattice 的Bilateral filter 源码及部分注释【来自于网络】
实现基于论文《Fast High-Dimensional Filtering Using the Permutohedral Lattice》 .
延伸阅读 saliency filte...
基于Permutohedral Lattice 的Bilateral filter 源码及部分注释【来自于网络】
实现基于论文《Fast High-Dimensional Filtering Using the Permutohedral Lattice》 .
延伸阅读 saliency filters精读之permutohedral lattice
1.bilateralPermutohedral 方法：
-
static Mat bilateralPermutohedral(Mat img, Mat edge, float sigma_s, float sigma_r) // img 和 edge 都必须是CV_32F类型
-
{
-
-
float invSpatialStdev = 1.0f / sigma_s;
-
float invColorStdev = 1.0f / sigma_r;
-
-
// Construct the position vectors out of x, y, r, g, and b.
-
int height = img.rows;
-
int width = img.cols;
-
-
int eCh = edge.channels(); // 1 或 3
-
int iCh = img.channels();
-
Image positions(1, width, height, 2 + eCh); // 只有一个子窗口
-
Image input(1, width, height, iCh);
-
-
//From Mat to Image
-
for (int y = 0; y < height; y++)
-
{
-
float *pimg = img.ptr<float>(y);
-
float *pedge = edge.ptr<float>(y);
-
for (int x = 0; x < width; x++)
-
{
-
// 参考论文 p4 3.1
-
// 5维的 positiion vector
-
positions(x, y)[0] = invSpatialStdev * x; // 0
-
positions(x, y)[1] = invSpatialStdev * y; // 1
-
-
for(int c = 0; c < eCh; c++)
-
positions(x, y)[2 + c] = invColorStdev * pedge[x * eCh + c]; // 2+
-
// 3维的 input vector
-
for(int c = 0; c < iCh; c++)
-
input(x, y)[c] = pimg[x * iCh + c];
-
}
-
}
-
-
// Filter the input with respect to the position vectors. (see permutohedral.h)
-
Image out = PermutohedralLattice::filter(input, positions);
-
-
// Save the result
-
Mat imgOut(img.size(), img.type());
-
for (int y = 0; y < height; y++)
-
{
-
float *pimgOut = imgOut.ptr<float>(y);
-
for (int x = 0; x < width; x++)
-
{
-
for(int c = 0; c < iCh; c++)
-
pimgOut[x * iCh + c] = out(x, y)[c];
-
}
-
}
-
-
return imgOut;
-
}
2. PermutohedralLattice 类：
-
/***************************************************************/
-
/* The algorithm class that performs the filter
-
*
-
* PermutohedralLattice::filter(...) does all the work.
-
*
-
*/
-
/***************************************************************/
-
class PermutohedralLattice
-
{
-
public:
-
-
/* Filters given image against a reference image.
-
* im : image to be bilateral-filtered. (input vector)
-
* ref : reference image whose edges are to be respected. (position vector)
-
*/
-
static Image filter(Image im, Image ref)
-
{
-
//timeval t[5];
-
-
// Create lattice
-
// gettimeofday(t+0, NULL);
-
// d = ref.channels (5)
-
// vd = im.channels + 1 (3+1)
-
PermutohedralLattice lattice(ref.channels, im.channels + 1, im.width * im.height * im.frames);
-
-
// Splat into the lattice
-
// gettimeofday(t+1, NULL);
-
// printf("Splatting...\n");
-
-
float *col = new float[im.channels + 1];
-
col[im.channels] = 1; // homogeneous coordinate
-
-
float *imPtr = im(0, 0, 0);
-
float *refPtr = ref(0, 0, 0); // position vector
-
for (int t = 0; t < im.frames; t++)
-
{
-
for (int y = 0; y < im.height; y++)
-
{
-
for (int x = 0; x < im.width; x++)
-
{
-
for (int c = 0; c < im.channels; c++)
-
{
-
col[c] = *imPtr++;
-
}
-
lattice.splat(refPtr, col);
-
refPtr += ref.channels;
-
}
-
}
-
}
-
-
// Blur the lattice
-
// gettimeofday(t+2, NULL);
-
// printf("Blurring...");
-
lattice.blur();
-
-
// Slice from the lattice
-
// gettimeofday(t+3, NULL);
-
// printf("Slicing...\n");
-
Image out(im.frames, im.width, im.height, im.channels);
-
-
lattice.beginSlice();
-
float *outPtr = out(0, 0, 0);
-
for (int t = 0; t < im.frames; t++)
-
{
-
for (int y = 0; y < im.height; y++)
-
{
-
for (int x = 0; x < im.width; x++)
-
{
-
lattice.slice(col);
-
float scale = 1.0f / col[im.channels];
-
for (int c = 0; c < im.channels; c++)
-
{
-
*outPtr++ = col[c] * scale;
-
}
-
}
-
}
-
}
-
-
// Print time elapsed for each step
-
// gettimeofday(t+4, NULL);
-
// const char *names[4] = {"Init ", "Splat ", "Blur ", "Slice "};
-
// for (int i = 1; i < 5; i++)
-
// printf("%s: %3.3f ms\n", names[i-1], (t[i].tv_sec - t[i-1].tv_sec) +
-
// (t[i].tv_usec - t[i-1].tv_usec)/1000000.0);
-
-
return out;
-
}
-
-
/* Constructor
-
* d_ : dimensionality of key vectors (ref.channels)
-
* vd_ : dimensionality of value vectors (im.channels + 1)
-
* nData_ : number of points in the input (im.size * im.frames)
-
*/
-
PermutohedralLattice(int d_, int vd_, int nData_) :
-
d(d_), vd(vd_), nData(nData_), hashTable(d_, vd_)
-
{
-
// Allocate storage for various arrays
-
elevated = new float[d + 1];
-
scaleFactor = new float[d];
-
-
greedy = new short[d + 1];
-
rank = new char[d + 1];
-
barycentric = new float[d + 2];
-
replay = new ReplayEntry[nData * (d + 1)];
-
nReplay = 0;
-
canonical = new short[(d + 1) * (d + 1)];
-
key = new short[d + 1];
-
-
// compute the coordinates of the canonical simplex, in which
-
// the difference between a contained point and the zero
-
// remainder vertex is always in ascending order. (See pg.4 of paper.)
-
// 论文第四页，d=4的矩阵例子（列主序）
-
for (int i = 0; i <= d; i++)
-
{
-
for (int j = 0; j <= d - i; j++)
-
canonical[i * (d + 1) + j] = i;
-
for (int j = d - i + 1; j <= d; j++)
-
canonical[i * (d + 1) + j] = i - (d + 1);
-
}
-
-
// Compute parts of the rotation matrix E. (See pg.4-5 of paper.)
-
for (int i = 0; i < d; i++)
-
{
-
// the diagonal entries for normalization
-
scaleFactor[i] = 1.0f / (sqrtf( (float)(i + 1) * (i + 2) ));
-
-
/* We presume that the user would like to do a Gaussian blur of standard deviation
-
* 1 in each dimension (or a total variance of d, summed over dimensions.)
-
* Because the total variance of the blur performed by this algorithm is not d,
-
* we must scale the space to offset this.
-
*
-
* The total variance of the algorithm is (See pg.6 and 10 of paper):
-
* [variance of splatting] + [variance of blurring] + [variance of splatting]
-
* = d(d+1)(d+1)/12 + d(d+1)(d+1)/2 + d(d+1)(d+1)/12
-
* = 2d(d+1)(d+1)/3.
-
*
-
* So we need to scale the space by (d+1)sqrt(2/3).
-
*/
-
// 论文 第四页 scale position vector
-
scaleFactor[i] *= (d + 1) * sqrtf(2.0 / 3);
-
}
-
}
-
-
-
/* Performs splatting with given position and value vectors */
-
// position: d-dimension position vector
-
// value: [r, g, b, 1]
-
void splat(float *position, float *value)
-
{
-
// first rotate position into the (d+1)-dimensional hyperplane
-
// 论文 第五页 Ex计算
-
elevated[d] = -d * position[d - 1] * scaleFactor[d - 1];
-
for (int i = d - 1; i > 0; i--)
-
elevated[i] = (elevated[i + 1] -
-
i * position[i - 1] * scaleFactor[i - 1] +
-
(i + 2) * position[i] * scaleFactor[i]);
-
elevated[0] = elevated[1] + 2 * position[0] * scaleFactor[0];
-
-
// prepare to find the closest lattice points
-
float scale = 1.0f / (d + 1);
-
char *myrank = rank;
-
short *mygreedy = greedy;
-
-
// greedily search for the closest zero-colored lattice point
-
// 论文 第三页
-
int sum = 0;
-
for (int i = 0; i <= d; i++)
-
{
-
float v = elevated[i] * scale;
-
float up = ceilf(v) * (d + 1); // 查找最近的整数点，up / down
-
float down = floorf(v) * (d + 1);
-
-
if (up - elevated[i] < elevated[i] - down)
-
mygreedy[i] = (short)up;
-
else
-
mygreedy[i] = (short)down;
-
-
sum += mygreedy[i];
-
}
-
sum /= d + 1; // consistent remainder (d+1)
-
-
// rank differential to find the permutation between this simplex and the canonical one.
-
// (See pg. 3-4 in paper.)
-
// 相对差值小的rank++
-
memset(myrank, 0, sizeof(char) * (d + 1));
-
for (int i = 0; i < d; i++)
-
for (int j = i + 1; j <= d; j++)
-
if (elevated[i] - mygreedy[i] < elevated[j] - mygreedy[j])
-
myrank[i]++;
-
else
-
myrank[j]++;
-
-
if (sum > 0)
-
{
-
// sum too large - the point is off the hyperplane.
-
// need to bring down the ones with the smallest differential
-
for (int i = 0; i <= d; i++)
-
{
-
if (myrank[i] >= d + 1 - sum)
-
{
-
mygreedy[i] -= d + 1;
-
myrank[i] += sum - (d + 1);
-
}
-
else
-
myrank[i] += sum;
-
}
-
}
-
else if (sum < 0)
-
{
-
// sum too small - the point is off the hyperplane
-
// need to bring up the ones with largest differential
-
for (int i = 0; i <= d; i++)
-
{
-
if (myrank[i] < -sum)
-
{
-
mygreedy[i] += d + 1;
-
myrank[i] += (d + 1) + sum;
-
}
-
else
-
myrank[i] += sum;
-
}
-
}
-
-
// Compute barycentric coordinates (See pg.10 of paper.)
-
memset(barycentric, 0, sizeof(float) * (d + 2));
-
for (int i = 0; i <= d; i++)
-
{
-
barycentric[d - myrank[i]] += (elevated[i] - mygreedy[i]) * scale;
-
barycentric[d + 1 - myrank[i]] -= (elevated[i] - mygreedy[i]) * scale;
-
}
-
barycentric[0] += 1.0f + barycentric[d + 1];
-
-
// Splat the value into each vertex of the simplex, with barycentric weights.
-
for (int remainder = 0; remainder <= d; remainder++)
-
{
-
// Compute the location of the lattice point explicitly (all but the last coordinate - it's redundant because they sum to zero)
-
for (int i = 0; i < d; i++)
-
key[i] = mygreedy[i] + canonical[remainder * (d + 1) + myrank[i]];
-
-
// Retrieve pointer to the value at this vertex.
-
float *val = hashTable.lookup(key, true);
-
-
// Accumulate values with barycentric weight.
-
for (int i = 0; i < vd; i++)
-
val[i] += barycentric[remainder] * value[i];
-
-
// Record this interaction to use later when slicing
-
replay[nReplay].offset = val - hashTable.getValues();
-
replay[nReplay].weight = barycentric[remainder];
-
nReplay++;
-
-
}
-
}
-
-
// Prepare for slicing
-
void beginSlice()
-
{
-
nReplay = 0;
-
}
-
-
/* Performs slicing out of position vectors. Note that the barycentric weights and the simplex
-
* containing each position vector were calculated and stored in the splatting step.
-
* We may reuse this to accelerate the algorithm. (See pg. 6 in paper.)
-
*/
-
void slice(float *col)
-
{
-
float *base = hashTable.getValues();
-
for (int j = 0; j < vd; j++)
-
col[j] = 0;
-
for (int i = 0; i <= d; i++)
-
{
-
ReplayEntry r = replay[nReplay++];
-
for (int j = 0; j < vd; j++)
-
{
-
col[j] += r.weight * base[r.offset + j];
-
}
-
}
-
}
-
-
/* Performs a Gaussian blur along each projected axis in the hyperplane. */
-
void blur()
-
{
-
// Prepare arrays
-
short *neighbor1 = new short[d + 1];
-
short *neighbor2 = new short[d + 1];
-
float *newValue = new float[vd * hashTable.size()];
-
float *oldValue = hashTable.getValues();
-
float *hashTableBase = oldValue;
-
-
float *zero = new float[vd];
-
for (int k = 0; k < vd; k++)
-
zero[k] = 0;
-
-
// For each of d+1 axes,
-
for (int j = 0; j <= d; j++)
-
{
-
printf("blur %d\t", j);
-
fflush(stdout);
-
-
// For each vertex in the lattice,
-
for (int i = 0; i < hashTable.size(); i++) // blur point i in dimension j
-
{
-
short *key = hashTable.getKeys() + i * (d); // keys to current vertex
-
for (int k = 0; k < d; k++)
-
{
-
neighbor1[k] = key[k] + 1;
-
neighbor2[k] = key[k] - 1;
-
}
-
neighbor1[j] = key[j] - d;
-
neighbor2[j] = key[j] + d; // keys to the neighbors along the given axis.
-
-
float *oldVal = oldValue + i * vd;
-
float *newVal = newValue + i * vd;
-
-
float *vm1, *vp1;
-
//printf("first neighbor\n");
-
vm1 = hashTable.lookup(neighbor1, false); // look up first neighbor
-
if (vm1)
-
vm1 = vm1 - hashTableBase + oldValue;
-
else
-
vm1 = zero;
-
//printf("second neighbor\n");
-
vp1 = hashTable.lookup(neighbor2, false); // look up second neighbor
-
if (vp1)
-
vp1 = vp1 - hashTableBase + oldValue;
-
else
-
vp1 = zero;
-
-
// Mix values of the three vertices
-
for (int k = 0; k < vd; k++)
-
newVal[k] = (0.25f * vm1[k] + 0.5f * oldVal[k] + 0.25f * vp1[k]);
-
}
-
float *tmp = newValue;
-
newValue = oldValue;
-
oldValue = tmp;
-
// the freshest data is now in oldValue, and newValue is ready to be written over
-
}
-
-
// depending where we ended up, we may have to copy data
-
if (oldValue != hashTableBase)
-
{
-
memcpy(hashTableBase, oldValue, hashTable.size()*vd * sizeof(float));
-
delete oldValue;
-
}
-
else
-
{
-
delete newValue;
-
}
-
printf("\n");
-
-
delete zero;
-
delete neighbor1;
-
delete neighbor2;
-
}
-
-
private:
-
-
int d, vd, nData;
-
float *elevated, *scaleFactor, *barycentric;
-
short *canonical;
-
short *key;
-
-
// slicing is done by replaying splatting (ie storing the sparse matrix)
-
struct ReplayEntry
-
{
-
int offset;
-
float weight;
-
} *replay;
-
int nReplay, nReplaySub;
-
-
public:
-
char *rank;
-
short *greedy;
-
HashTablePermutohedral hashTable;
-
};
3. 用于permutohedral lattice的哈希表：
-
/***************************************************************/
-
/* Hash table implementation for permutohedral lattice
-
*
-
* The lattice points are stored sparsely using a hash table.
-
* The key for each point is its spatial location in the (d+1)-
-
* dimensional space.
-
*/
-
/***************************************************************/
-
-
class HashTablePermutohedral
-
{
-
public:
-
/* Constructor
-
* kd_: the dimensionality of the position vectors on the hyperplane.
-
* vd_: the dimensionality of the value vectors
-
*/
-
HashTablePermutohedral(int kd_, int vd_) : kd(kd_), vd(vd_)
-
{
-
capacity = 1 << 15;
-
filled = 0;
-
entries = new Entry[capacity];
-
keys = new short[kd * capacity / 2]; // 多维 键－值对
-
values = new float[vd * capacity / 2];
-
memset(values, 0, sizeof(float)*vd * capacity / 2);
-
}
-
-
// Returns the number of vectors stored.
-
int size()
-
{
-
return filled;
-
}
-
-
// Returns a pointer to the keys array.
-
short *getKeys()
-
{
-
return keys;
-
}
-
-
// Returns a pointer to the values array.
-
float *getValues()
-
{
-
return values;
-
}
-
-
/* Returns the index into the hash table for a given key.
-
* key: a pointer to the position vector.
-
* h: hash of the position vector.
-
* create: a flag specifying whether an entry should be created,
-
* should an entry with the given key not found.
-
*/
-
// 返回 value 指针的偏移量
-
int lookupOffset(short *key, size_t h, bool create = true)
-
{
-
-
// Double hash table size if necessary
-
// 如果存储的数据达到或超过容量的一半
-
if (filled >= (capacity / 2) - 1)
-
{
-
grow();
-
}
-
-
// Find the entry with the given key
-
// 根据给定的 hash 索引 entry
-
while (1)
-
{
-
Entry e = entries[h];
-
// check if the cell is empty
-
// 检查该 entry 的 key 是否存在
-
if (e.keyIdx == -1)
-
{
-
if (!create)
-
return -1; // Return not found.
-
-
// need to create an entry. Store the given key.
-
for (int i = 0; i < kd; i++)
-
keys[filled * kd + i] = key[i];
-
-
e.keyIdx = filled * kd;
-
e.valueIdx = filled * vd;
-
entries[h] = e;
-
filled++;
-
-
return e.valueIdx;
-
}
-
-
// check if the cell has a matching key
-
bool match = true;
-
for (int i = 0; i < kd && match; i++)
-
match = keys[e.keyIdx + i] == key[i];
-
if (match)
-
return e.valueIdx;
-
-
// increment the bucket with wraparound
-
// 顺序查找下一个 entry 【计算出的hash值相同的情况】
-
h++;
-
// 如果到达最后一个 entry， 则从第一个 entry 开始找
-
if (h == capacity)
-
h = 0;
-
}
-
}
-
-
/* Looks up the value vector associated with a given key vector.
-
* k : pointer to the key vector to be looked up.
-
* create : true if a non-existing key should be created.
-
*/
-
float *lookup(short *k, bool create = true)
-
{
-
size_t h = hash(k) % capacity;
-
int offset = lookupOffset(k, h, create);
-
if (offset < 0)
-
return NULL;
-
else
-
return values + offset;
-
};
-
-
/* Hash function used in this implementation. A simple base conversion. */
-
size_t hash(const short *key)
-
{
-
size_t k = 0;
-
for (int i = 0; i < kd; i++)
-
{
-
k += key[i];
-
k *= 2531011;
-
}
-
return k;
-
}
-
-
private:
-
/* Grows the size of the hash table */
-
void grow()
-
{
-
printf("Resizing hash table\n");
-
-
size_t oldCapacity = capacity;
-
capacity *= 2; // 变为2倍容量
-
-
// Migrate the value vectors.
-
float *newValues = new float[vd * capacity / 2];
-
memset(newValues, 0, sizeof(float)*vd * capacity / 2);
-
memcpy(newValues, values, sizeof(float)*vd * filled);
-
delete[] values;
-
values = newValues;
-
-
// Migrate the key vectors.
-
short *newKeys = new short[kd * capacity / 2];
-
memcpy(newKeys, keys, sizeof(short)*kd * filled);
-
delete[] keys;
-
keys = newKeys;
-
-
Entry *newEntries = new Entry[capacity];
-
-
// Migrate the table of indices.
-
for (size_t i = 0; i < oldCapacity; i++)
-
{
-
if (entries[i].keyIdx == -1)
-
continue;
-
// 根据键值计算hash
-
size_t h = hash(keys + entries[i].keyIdx) % capacity;
-
// 如果hash对应entry的keyidx已经被占用，则顺序往后找 entry，直到发现该 entry 的 keyidx 未被占用
-
while (newEntries[h].keyIdx != -1)
-
{
-
h++;
-
if (h == capacity)
-
h = 0;
-
}
-
newEntries[h] = entries[i];
-
}
-
delete[] entries;
-
entries = newEntries;
-
}
-
-
// Private struct for the hash table entries.
-
struct Entry
-
{
-
Entry() : keyIdx(-1), valueIdx(-1) {}
-
int keyIdx; // keys 的索引
-
int valueIdx; // values 的索引
-
};
-
-
short *keys;
-
float *values;
-
Entry *entries;
-
size_t capacity, filled; // 分别表示 entry 的容量 和 已填充的 entry 数
-
int kd, vd; // keys 和 values 数组的维度（PermutohedraLattice 会将数据 splat 到高维空间）
-
};
效果图：
文章来源: panda1234lee.blog.csdn.net，作者：panda1234lee，版权归原作者所有，如需转载，请联系作者。
原文链接：panda1234lee.blog.csdn.net/article/details/52892359
【版权声明】本文为华为云社区用户转载文章，如果您发现本社区中有涉嫌抄袭的内容，欢迎发送邮件进行举报，并提供相关证据，一经查实，本社区将立刻删除涉嫌侵权内容，举报邮箱：
cloudbbs@huaweicloud.com
- 点赞
- 收藏
- 关注作者
评论（0）