264 lines
7.3 KiB
C++
264 lines
7.3 KiB
C++
#pragma once
|
|
|
|
#include <cmath>
|
|
#include <vector>
|
|
|
|
#include "BoxSizes.h"
|
|
#include "Image3D.h"
|
|
|
|
template <class T>
|
|
struct BoxGaus3D
|
|
{
|
|
void boxFilter(std::vector<T>& input, unsigned w, unsigned h, unsigned d, unsigned filterSize)
|
|
{
|
|
std::vector<T> buffer(input.size());
|
|
|
|
boxFilter(input.data(), buffer.data(), w, h, d, filterSize);
|
|
|
|
input.assign(buffer.begin(), buffer.end());
|
|
}
|
|
|
|
void boxFilter(T* input, T* output, unsigned w, unsigned h, unsigned d, unsigned filterSize)
|
|
{
|
|
assertMsg((filterSize % 2 == 1), "filterSize must be odd");
|
|
unsigned radius = filterSize / 2;
|
|
|
|
ImageView3D<T> in (w, h, d, input);
|
|
ImageView3D<T> out(w, h, d, output);
|
|
|
|
boxBlur(in, out, radius, radius, radius);
|
|
}
|
|
|
|
void approxGaus(Image3D<T>& input, _Point3<T> sigma, unsigned nFilt)
|
|
{
|
|
approxGaus(input.data(), input.getWidth(), input.getHeight(), input.getDepth(), sigma, nFilt);
|
|
}
|
|
|
|
void approxGaus(std::vector<T>& input, unsigned w, unsigned h, unsigned d, _Point3<T> sigma, unsigned nFilt)
|
|
{
|
|
std::vector<T> buffer(input.size());
|
|
|
|
approxGaus(input.data(), buffer.data(), w, h, d, sigma, nFilt);
|
|
|
|
input.assign(buffer.begin(), buffer.end());
|
|
}
|
|
|
|
void approxGaus(T* input, T* output, unsigned w, unsigned h, unsigned d, _Point3<T> sigma, unsigned nFilt)
|
|
{
|
|
ImageView3D<T> in (w, h, d, input);
|
|
ImageView3D<T> out(w, h, d, output);
|
|
|
|
approxGaus(in, out, sigma, nFilt);
|
|
}
|
|
|
|
void approxGaus(ImageView3D<T>& in, ImageView3D<T>& out, _Point3<T> sigma, unsigned nFilt)
|
|
{
|
|
BoxSizes<T> bsX(sigma.x, nFilt);
|
|
BoxSizes<T> bsY(sigma.y, nFilt);
|
|
BoxSizes<T> bsZ(sigma.z, nFilt);
|
|
|
|
assertMsg((nFilt % 2 == 1), "nFilt must be odd");
|
|
|
|
// if equal, we can save some cond's inside the loop
|
|
if (bsX.m == bsY.m && bsX.m == bsZ.m)
|
|
{
|
|
const size_t m = bsX.m;
|
|
|
|
for (size_t i = 1; i <= m; i++)
|
|
{
|
|
if (i % 2) {
|
|
boxBlur(in, out, bsX.wl, bsY.wl, bsZ.wl);
|
|
} else {
|
|
boxBlur(out, in, bsX.wl, bsY.wl, bsZ.wl);
|
|
}
|
|
}
|
|
|
|
for (size_t i = 1; i <= nFilt - m; i++)
|
|
{
|
|
if (i % 2) {
|
|
boxBlur(in, out, bsX.wu, bsY.wu, bsZ.wu);
|
|
} else {
|
|
boxBlur(out, in, bsX.wu, bsY.wu, bsZ.wu);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (size_t i = 1; i <= nFilt; i++)
|
|
{
|
|
const size_t rX = (i < bsX.m ? bsX.wl : bsX.wu);
|
|
const size_t rY = (i < bsY.m ? bsY.wl : bsY.wu);
|
|
const size_t rZ = (i < bsZ.m ? bsZ.wl : bsZ.wu);
|
|
|
|
if (i % 2) {
|
|
boxBlur(in, out, rX, rY, rZ);
|
|
} else{
|
|
boxBlur(out, in, rX, rY, rZ);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
private:
|
|
void boxBlur(ImageView3D<T>& in, ImageView3D<T>& out, size_t rX, size_t rY, size_t rZ)
|
|
{
|
|
//box(in, out, r);
|
|
boxBlurX(in, out, rX);
|
|
boxBlurY(out, in, rY);
|
|
boxBlurZ(in, out, rZ);
|
|
}
|
|
|
|
void boxBlurX(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
|
|
{
|
|
T iarr = (T)1.0 / (r + r + 1);
|
|
|
|
for (size_t z = 0; z < src.getDepth(); z++)
|
|
{
|
|
for (size_t y = 0; y < src.getHeight(); y++)
|
|
{
|
|
auto srcView = src.getViewYZ(y, z);
|
|
auto dstView = dst.getViewYZ(y, z);
|
|
|
|
boxBlur1D<SlicePlane::YZ>(srcView, dstView, src.getWidth(), r, iarr);
|
|
}
|
|
}
|
|
}
|
|
|
|
void boxBlurY(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
|
|
{
|
|
T iarr = (T)1.0 / (r + r + 1);
|
|
|
|
for (size_t x = 0; x < src.getWidth(); x++)
|
|
{
|
|
for (size_t z = 0; z < src.getDepth(); z++)
|
|
{
|
|
auto srcView = src.getViewXZ(x, z);
|
|
auto dstView = dst.getViewXZ(x, z);
|
|
|
|
boxBlur1D<SlicePlane::XZ>(srcView, dstView, src.getHeight(), r, iarr);
|
|
}
|
|
}
|
|
}
|
|
|
|
void boxBlurZ(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
|
|
{
|
|
T iarr = (T)1.0 / (r + r + 1);
|
|
|
|
for (size_t x = 0; x < src.getWidth(); x++)
|
|
{
|
|
for (size_t y = 0; y < src.getHeight(); y++)
|
|
{
|
|
auto srcView = src.getViewXY(x, y);
|
|
auto dstView = dst.getViewXY(x, y);
|
|
|
|
boxBlur1D<SlicePlane::XY>(srcView, dstView, src.getDepth(), r, iarr);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template <SlicePlane S>
|
|
void boxBlur1D(const typename ImageView3D<T>::template ConstLineView<S>& src, typename ImageView3D<T>::template LineView<S>& dst, size_t len, size_t r, T iarr)
|
|
{
|
|
size_t li = 0; // left index
|
|
size_t ri = r; // right index
|
|
|
|
T fv = src(0);
|
|
T lv = src(len - 1);
|
|
T val = (r + 1)*fv;
|
|
|
|
for (size_t i = 0; i < r; i++)
|
|
val += src(i);
|
|
|
|
// Überhangbereich links vom Bild
|
|
for (size_t i = 0; i <= r; i++)
|
|
{
|
|
val += src(ri++) - fv;
|
|
dst(i) = val*iarr;
|
|
}
|
|
|
|
// Bildbereich
|
|
for (size_t i = r + 1; i < len - r; i++)
|
|
{
|
|
val += src(ri++) - src(li++);
|
|
dst(i) = val*iarr;
|
|
}
|
|
|
|
// Überhangbereich rechts vom Bild
|
|
for (size_t i = len - r; i < len; i++)
|
|
{
|
|
val += lv - src(li++);
|
|
dst(i) = val*iarr;
|
|
}
|
|
}
|
|
|
|
void box(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
|
|
{
|
|
T iarr = (T)1.0 / (r + r + 1);
|
|
|
|
ImageView3D<T> buffer(src);
|
|
ImageView3D<T> buffer2(dst.getHeight(), dst.getDepth(), dst.getWidth(), dst.val_begin());
|
|
|
|
for (size_t i = 0; i < 3; i++)
|
|
{
|
|
for (size_t idx2 = 0; idx2 < buffer.getHeight(); idx2++)
|
|
{
|
|
for (size_t idx3 = 0; idx3 < buffer.getDepth(); idx3++)
|
|
{
|
|
box1D(buffer, buffer2, buffer.getWidth(), idx2, idx3, r, iarr);
|
|
}
|
|
}
|
|
|
|
// reinteprete buffer sizes after rotation and swap data pointers
|
|
buffer = ImageView3D<T>(buffer.getHeight(), buffer.getDepth(), buffer.getWidth(), buffer.val_begin());
|
|
buffer2 = ImageView3D<T>(buffer2.getHeight(), buffer2.getDepth(), buffer2.getWidth(), buffer2.val_begin());
|
|
buffer.swap(buffer2);
|
|
}
|
|
|
|
buffer.swap(dst);
|
|
}
|
|
|
|
void box1D(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t len, size_t idx2, size_t idx3, size_t r, T iarr)
|
|
{
|
|
for (size_t i = 0; i < len; i++)
|
|
{
|
|
// dst(idx2, idx3, i) = src(i, idx2, idx3);
|
|
|
|
size_t li = 0; // left index
|
|
size_t ri = r; // right index
|
|
|
|
T fv = src(0, idx2, idx3);
|
|
T lv = src(len - 1, idx2, idx3);
|
|
T val = (r + 1)*fv;
|
|
|
|
for (size_t i = 0; i < r; i++)
|
|
val += src(i, idx2, idx3);
|
|
|
|
// Überhangbereich links vom Bild
|
|
for (size_t i = 0; i <= r; i++)
|
|
{
|
|
val += src(ri++, idx2, idx3) - fv;
|
|
dst(idx2, idx3, i) = val*iarr;
|
|
}
|
|
|
|
// Bildbereich
|
|
for (size_t i = r + 1; i < len - r; i++)
|
|
{
|
|
val += src(ri++, idx2, idx3) - src(li++, idx2, idx3);
|
|
dst(idx2, idx3, i) = val*iarr;
|
|
}
|
|
|
|
// Überhangbereich rechts vom Bild
|
|
for (size_t i = len - r; i < len; i++)
|
|
{
|
|
val += lv - src(li++, idx2, idx3);
|
|
dst(idx2, idx3, i) = val*iarr;
|
|
}
|
|
}
|
|
}
|
|
};
|
|
|
|
|
|
template struct BoxGaus3D<float>;
|
|
template struct BoxGaus3D<double>;
|