Artifact 868cf8d0a5e9c4efcb655fa440980b09e8fec0f1:
- File
operator/bilateral-grey8.crimp
— part of check-in
[1da2e834b8]
at
2011-11-18 02:24:25
on branch infinite-plane
— New branch.
Embed images into the infinite 2d plane, i.e. make them bounded rectangles with a location in that plane. We had a start on that with the operators for projective warp, storing the origin point in the user-meta data. This information is now in the C level structures.
This commit provides the basic change, plus new accessors for the geometry data, and updates all users of said geometry data in C to use these new accessors.
The semantics of the various operations which could make use of the location information has not changed yet. This is for future commits in this branch. (user: andreask size: 6898) [more...]
bilateral_grey8 Tcl_Obj* imageObj double sigma_space double sigma_range /* * Bilateral filter. Uses a bilateral grid downsampled by sigma-s and * sigma-r for higher performance, and lesser memory use. A sigma of 1 * implies 'no downsampling', filtering is on the full grid: high * memory use, slow speed. */ crimp_image* result; crimp_image* image; crimp_volume* wi; /* Bilateral grid, accumulated pixel intensities */ crimp_volume* w; /* Bilateral grid, accumulated pixel counts, = weight factor */ int x, y, z; int bgrid_width, bgrid_height, bgrid_range, bgrid_maxdim; double* nw; double* nwi; /* * Process and validate the arguments. */ crimp_input (imageObj, image, grey8); CRIMP_ASSERT (sigma_space >= 1, "Cannot use sigma/s < 1"); CRIMP_ASSERT (sigma_range >= 1, "Cannot use sigma/r < 1"); result = crimp_new_like (image); /* * Determine the size of the bilateral grid. * +1 = One more, in case the float->int of the ceil result rounded down. * +4 = Borders for the convolution of the grid. * * TODO NOTE: The SParis BF code obtains the min and max grey levels from the * TODO NOTE: image and uses that for the range, instead of a fixed 256 (Also * TODO NOTE: assumes that intensity is in [0,1]). */ bgrid_width = 4 + 1 + (int) ceil (crimp_w (image)/sigma_space); bgrid_height = 4 + 1 + (int) ceil (crimp_w (image)/sigma_space); bgrid_range = 4 + 1 + (int) ceil (256/sigma_range); bgrid_maxdim = MAX (bgrid_width, MAX (bgrid_height, bgrid_range)); /* * Phase I. Allocate and initialize the bilateral grid (2 volumes). */ wi = crimp_vnew_float (bgrid_width, bgrid_height, bgrid_range); w = crimp_vnew_float (bgrid_width, bgrid_height, bgrid_range); for (z = 0; z < bgrid_range; z++) { for (y = 0; y < bgrid_height; y++) { for (x = 0; x < bgrid_width; x++) { VFLOATP (wi, x, y, z) = 0.0; VFLOATP (w, x, y, z) = 0.0; } } } /* * Phase II. Update the bilateral grid with the downsampled image data. */ for (y = 0; y < crimp_h (image); y++) { for (x = 0; x < crimp_w (image); x++) { double p = GREY8 (image, x, y); /* +2 is the offset to keep the borders empty. */ int xr = 2 + lrint (((double) x) / sigma_space); int yr = 2 + lrint (((double) y) / sigma_space); int pr = 2 + lrint (p / sigma_range); VFLOATP (wi, xr, yr, pr) += p; VFLOATP (w, xr, yr, pr) += 1; } } /* * Phase III. Convolve the grid using gaussian [1/16 (1 4 6 4 1)] along each * of the three axes. The convolution is hard-wired. Note that the grid was * allocated with the necessary borders, and the previous phases made sure to * keep the borders empty. * * NOTE: I sort of tried to optimize here, using a buffer just for a single * slice through the grid. The SParis code creates a whole duplicate of the * grid. It doesn't have to access the grid in strided fashion either, except * for the neighbour calculations. * * It also uses a simpler gaussian (1/4 (1 2 1)), and applies it twice. */ nw = CRIMP_ALLOC_ARRAY (bgrid_maxdim, double); /* Helper arrays to buffer the convolution */ nwi = CRIMP_ALLOC_ARRAY (bgrid_maxdim, double); /* result per scan line. */ /* gauss(a,b,c,d,e) = 1a+4b+6c+4d+1e = a+e+4(b+d)+6c = a+e+4(b+d+c)+2c */ #define GAUSS(a, b, c, d, e) ((((a)+(e)) + 4*((b)+(d)) + 6*(c))/16.) #define GX(f, x, y, z) \ GAUSS (VFLOATP (f, x-2, y, z), \ VFLOATP (f, x-1, y, z), \ VFLOATP (f, x , y, z), \ VFLOATP (f, x+1, y, z), \ VFLOATP (f, x+2, y, z)) #define GY(f, x, y, z) \ GAUSS (VFLOATP (f, x, y-2, z), \ VFLOATP (f, x, y-1, z), \ VFLOATP (f, x, y , z), \ VFLOATP (f, x, y+1, z), \ VFLOATP (f, x, y+2, z)) #define GZ(f, x, y, z) \ GAUSS (VFLOATP (f, x, y, z-2), \ VFLOATP (f, x, y, z-1), \ VFLOATP (f, x, y, z ), \ VFLOATP (f, x, y, z+1), \ VFLOATP (f, x, y, z+2)) /* Gaussian @ X */ for (z = 2; z < bgrid_range-2; z++) { for (y = 2; y < bgrid_height-2; y++) { for (x = 2; x < bgrid_width-2; x++) { nw [x-2] = GX(w, x, y, z); nwi [x-2] = GX(wi, x, y, z); } for (x = 2; x < bgrid_width-2; x++) { VFLOATP (w, x, y, z) = nw [x-2]; } for (x = 2; x < bgrid_width-2; x++) { VFLOATP (wi, x, y, z) = nwi[x-2]; } } } /* Gaussian @ Y */ for (z = 2; z < bgrid_range-2; z++) { for (x = 2; x < bgrid_width-2; x++) { for (y = 2; y < bgrid_height-2; y++) { nw [y-2] = GY(w, x, y, z); nwi [y-2] = GY(wi, x, y, z); } for (y = 2; y < bgrid_height-2; y++) { VFLOATP (w, x, y, z) = nw [y-2]; } for (y = 2; y < bgrid_height-2; y++) { VFLOATP (wi, x, y, z) = nwi[y-2]; } } } /* Gaussian @ Z */ for (y = 2; y < bgrid_height-2; y++) { for (x = 2; x < bgrid_width-2; x++) { for (z = 2; z < bgrid_range-2; z++) { nw [z-2] = GZ(w, x, y, z); nwi [z-2] = GZ(wi, x, y, z); } for (z = 2; z < bgrid_range-2; z++) { VFLOATP (w, x, y, z) = nw [z-2]; } for (z = 2; z < bgrid_range-2; z++) { VFLOATP (wi, x, y, z) = nwi[z-2]; } } } #undef GX #undef GY #undef GZ #undef GAUSS /* * Phase IV. Resample the image using the updated bilateral grid and trilinear * interpolation. * * #define I(a,b,s) ((b) + ((a)-(b))*(s)) */ #define BETWEEN(a,b,s) ((a)*(s) + (b)*(1-(s))) for (y = 0; y < crimp_h (image); y++) { for (x = 0; x < crimp_w (image); x++) { double winew, wnew, p = GREY8 (image, x, y); /* Continuous grid location */ double xf = 2 + ((double) x) / sigma_space; double yf = 2 + ((double) y) / sigma_space; double pf = 2 + p / sigma_range; /* Integral grid location */ int xr = lrint (xf); int yr = lrint (yf); int pr = lrint (pf); /* Fractional grid location from the integral */ if (xr > xf) { xr -- ; } ; xf = xf - xr; if (yr > yf) { yr -- ; } ; yf = yf - yr; if (pr > pf) { pr -- ; } ; pf = pf - pr; /* Trilinear interpolate over the grid */ winew = BETWEEN (BETWEEN (BETWEEN (VFLOATP (wi, xr, yr, pr), VFLOATP (wi, xr+1, yr, pr), xf), BETWEEN (VFLOATP (wi, xr, yr+1, pr), VFLOATP (wi, xr+1, yr+1, pr), xf), yf), BETWEEN (BETWEEN (VFLOATP (wi, xr, yr, pr+1), VFLOATP (wi, xr+1, yr, pr+1), xf), BETWEEN (VFLOATP (wi, xr, yr+1, pr+1), VFLOATP (wi, xr+1, yr+1, pr+1), xf), yf), pf); wnew = BETWEEN (BETWEEN (BETWEEN (VFLOATP (w, xr, yr, pr), VFLOATP (w, xr+1, yr, pr), xf), BETWEEN (VFLOATP (w, xr, yr+1, pr), VFLOATP (w, xr+1, yr+1, pr), xf), yf), BETWEEN (BETWEEN (VFLOATP (w, xr, yr, pr+1), VFLOATP (w, xr+1, yr, pr+1), xf), BETWEEN (VFLOATP (w, xr, yr+1, pr+1), VFLOATP (w, xr+1, yr+1, pr+1), xf), yf), pf); GREY8 (result, x, y) = CLAMP (0, (winew / wnew), 255); } } #undef BETWEEN Tcl_SetObjResult(interp, crimp_new_image_obj (result)); return TCL_OK; /* vim: set sts=4 sw=4 tw=80 et ft=c: */ /* * Local Variables: * mode: c * c-basic-offset: 4 * fill-column: 78 * End: */