static void __lmrcImageSmoothingMeanFilter(mrcImage *dst, mrcImage *src, lmrcImageSmoothingInfo *info, long mode);
static void __lmrcImageSmoothingSurfaceFitFilter(mrcImage *dst, mrcImage *src, lmrcImageSmoothingInfo *info, long mode);
static void __lmrcImageSmoothingLeeSigmaFilter(mrcImage *dst, mrcImage *src, lmrcImageSmoothingInfo *info, long mode);
+static void __lmrcImageSmoothingLocalWeightedMeanFilter(mrcImage *dst, mrcImage *src, lmrcImageSmoothingInfo *info, long mode);
void
lmrcImageSmoothingModeInfo(FILE* fpt)
fprintf(fpt, " 2: Mean filter \n");
fprintf(fpt, " 3: Surface Fit filter\n");
fprintf(fpt, " 4: Lee-Sigma filter\n");
+ fprintf(fpt, " 5: Local weighted mean filter\n");
fprintf(fpt, ">>> Mode2 <<<< \n");
fprintf(fpt, " 0: Kernel is cube. : required argumnent: -sx, -sy, [-sz]\n");
fprintf(fpt, " 1: Kernel is sphere.(Only mode 1 supported) : required argumnent: -r \n");
DEBUGPRINT("__lmrcImageSmoothingLeeSigmaFilter\n");
__lmrcImageSmoothingLeeSigmaFilter(dst, src, info, mode);
break;
+ case 5: /* Local Weighted Mean Filtering */
+ DEBUGPRINT("__lmrcImageSmoothingLocalWeightedMeanFilter\n");
+ __lmrcImageSmoothingLocalWeightedMeanFilter(dst, src, info, mode);
+ break;
default: {
fprintf(stderr, "Not supported mode in lmrcImageSmoothing: %ld\n", info->mode);
lmrcImageSmoothingModeInfo(stderr);
long n, i, index, center ,lsx ,lsy ,lsz;
float x, y, z ,r2, rr;
float tmpx, tmpy, tmpz;
- long offsetx, offsety, offsetz;
+ int offsetx, offsety, offsetz;
double data;
float* kernel;
float* kernel2;
break;
}
default: {
- fprintf(stderr, "Not supported mode: %d\n", mode);
+ fprintf(stderr, "Not supported mode: %ld\n", mode);
exit(EXIT_FAILURE);
break;
}
}
}
+void
+__lmrcImageSmoothingLocalWeightedMeanFilter(mrcImage* dst, mrcImage* src, lmrcImageSmoothingInfo* info, long mode)
+{
+ long n, index;
+ float x, y, z;
+ float tmpx, tmpy, tmpz;
+ long offsetx, offsety, offsetz ,lsx ,lsy ,lsz;
+ double data, sum, sumW, W;
+
+ if (mode==0){
+ lsx=floor(info->sx+0.5);
+ lsy=floor(info->sy+0.5);
+ lsz=floor(info->sz+0.5);
+ if(IsOdd(lsx)) {
+ offsetx = (lsx - 1)/2;
+ } else {
+ fprintf(stderr, "Not supported kernel size: sx = %ld\n", lsx);
+ exit(EXIT_FAILURE);
+ }
+ if(IsOdd(lsy)) {
+ offsety = (lsy - 1)/2;
+ } else {
+ fprintf(stderr, "Not supported kernel size: sy = %ld\n", lsy);
+ exit(EXIT_FAILURE);
+ }
+ if(IsOdd(lsz)) {
+ offsetz = (lsz - 1)/2;
+ } else {
+ fprintf(stderr, "Not supported kernel size: sz = %ld\n", lsz);
+ }
+ n = lsx*lsy*lsz;
+ }
+ if(src->HeaderN.z!=1) {
+ if(IsOdd(lsz)) {
+ offsetz = (lsz - 1)/2;
+ } else {
+ fprintf(stderr, "Not supported kernel size: sz = %ld\n", lsz);
+ }
+ n = n*lsz;
+ }
+ DEBUGPRINT5("x: %ld %ld y: %ld %ld datapoints: %ld\n", offsetx, lsx, offsety, lsy, n);
+ dst->Header = src->Header;
+ mrcInit(dst, NULL);
+ for(z = 0; z < src->HeaderN.z; z++) {
+ for(y = 0; y < src->HeaderN.y; y++) {
+ for(x = 0; x < src->HeaderN.x; x++) {
+ if(x < offsetx
+ || y < offsety
+ || z < offsetz
+ || src->HeaderN.x-1-offsetx<x
+ || src->HeaderN.y-1-offsety<y
+ || src->HeaderN.z-1-offsetz<z ) {
+ mrcPixelDataGet(src, x, y, 0.0, &data, mrcPixelRePart, mrcPixelHowNearest);
+ mrcPixelDataSet(dst, x, y, 0.0, data, mrcPixelRePart);
+ } else {
+ index = 0;
+ sum = 0.0;
+ sumW = 0.0;
+ for(tmpz=z-offsetz; tmpz<=z+offsetz; tmpz++) {
+ for(tmpy=y-offsety; tmpy<=y+offsety; tmpy++) {
+ for(tmpx=x-offsetx; tmpx<=x+offsetx; tmpx++) {
+ W = (offsetx - fabs(tmpx - x) + 1)
+ *(offsety - fabs(tmpy - y) + 1)
+ *(offsetz - fabs(tmpz - z) + 1);
+ mrcPixelDataGet(src, tmpx, tmpy, tmpz, &data, mrcPixelRePart, mrcPixelHowNearest);
+ sum += data*W;
+ sumW += W;
+ }
+ }
+ }
+ data = sum / sumW;
+ mrcPixelDataSet(dst, x, y, tmpz, data, mrcPixelRePart);
+ }
+ }
+ }
+ }
+}
+
static long int __SumOfSquare(long int);
static long int __SumOfSquareSquare(long int n, long int m);
static long int __SumOf4(long int n);
OSDN Git Service
