Panda3D
|
00001 // Filename: pnm-image-filter.cxx 00002 // Created by: 00003 // 00004 //////////////////////////////////////////////////////////////////// 00005 // 00006 // PANDA 3D SOFTWARE 00007 // Copyright (c) Carnegie Mellon University. All rights reserved. 00008 // 00009 // All use of this software is subject to the terms of the revised BSD 00010 // license. You should have received a copy of this license along 00011 // with this source code in a file named "LICENSE." 00012 // 00013 //////////////////////////////////////////////////////////////////// 00014 00015 // The functions in this module support spatial filtering of an image by 00016 // convolution with an (almost) arbitrary kernel. There are a broad class of 00017 // image filters to which this principle applies, including box filtering, 00018 // Bartlett, and Gaussian. 00019 00020 // This particular approach breaks the 2-d kernel into two 1-d kernels, which 00021 // improves performance at the expense of limiting the types of filters that 00022 // may be applied. In general, any of the square-based kernels are still 00023 // applicable; the only circle-based kernel that still works with this 00024 // approach is Gaussian. Furthermore, the implementation assume that the 00025 // kernel is symmetric about zero. 00026 00027 // The image is filtered first along one axis, then along the other. This 00028 // decreases the complexity of the convolution operation: it is faster to 00029 // convolve twice with a one-dimensional kernel than once with a two- 00030 // dimensional kernel. In the interim, a temporary matrix of type StoreType 00031 // (a numeric type, described below) is built which contains the results 00032 // from the first convolution. The entire process is then repeated for 00033 // each channel in the image. 00034 00035 #include "pandabase.h" 00036 #include <math.h> 00037 #include "cmath.h" 00038 00039 #include "pnmImage.h" 00040 00041 // WorkType is an abstraction that allows the filtering process to be 00042 // recompiled to use either floating-point or integer arithmetic. On SGI 00043 // machines, there doesn't seem to be much of a performance difference-- 00044 // if anything, the integer arithmetic is slower--though certainly other 00045 // architectures may differ. 00046 00047 // A StoreType is a numeric type that is used to store the intermediate values 00048 // of the filtering computation; temporary calculations are performed using 00049 // WorkTypes. source_max represents the largest value that will be stored in 00050 // a StoreType, while filter_max represents the largest value that will 00051 // multiplied by it as a weighted filter (source_max * filter_max must fit 00052 // comfortably within a WorkType, with room to spare). 00053 00054 // Floating-point arithmetic is slightly faster and slightly more precise, 00055 // but the main reason to use it is that it is conceptually easy to understand. 00056 // All values are scaled in the range 0..1, as they should be. 00057 00058 // The biggest reason to use integer arithmetic is space. A table of 00059 // StoreTypes must be allocated to match the size of the image. On an SGI, 00060 // sizeof(float) is 4, while sizeof(short) is 2 and sizeof(char) is, of 00061 // course, 1. Since the source precision is probably 8 bits anyway (there 00062 // are usually 8 bits per channel), it doesn't cost any precision at all to 00063 // use shorts, and not very much to use chars. 00064 00065 /* 00066 // To use double-precision floating point, 8 bytes: (strictly for the neurotic) 00067 typedef double WorkType; 00068 typedef double StoreType; 00069 static const WorkType source_max = 1.0; 00070 static const WorkType filter_max = 1.0; 00071 */ 00072 00073 // To use single-precision floating point, 4 bytes: 00074 typedef double WorkType; 00075 typedef float StoreType; 00076 static const WorkType source_max = 1.0; 00077 static const WorkType filter_max = 1.0; 00078 00079 /* 00080 // To use 16-bit integer arithmetic, 2 bytes: 00081 typedef unsigned long WorkType; 00082 typedef unsigned short StoreType; 00083 static const WorkType source_max = 65535; 00084 static const WorkType filter_max = 255; 00085 */ 00086 00087 /* 00088 // To use 8-bit integer arithmetic, 1 byte: 00089 typedef unsigned long WorkType; 00090 typedef unsigned char StoreType; 00091 static const WorkType source_max = 255; 00092 static const WorkType filter_max = 255; 00093 */ 00094 00095 00096 00097 // filter_row() filters a single row by convolving with a one-dimensional 00098 // kernel filter. The kernel is defined by an array of weights in filter[], 00099 // where the ith element of filter corresponds to abs(d * scale), if scale>1.0, 00100 // and abs(d), if scale<=1.0, where d is the offset from the center and varies 00101 // from -filter_width to filter_width. 00102 00103 // Note that filter_width is not necessarily the length of the array; it is 00104 // the radius of interest of the filter function. The array may need to be 00105 // larger (by a factor of scale), to adequately cover all the values. 00106 00107 static void 00108 filter_row(StoreType dest[], int dest_len, 00109 const StoreType source[], int source_len, 00110 double scale, // == dest_len / source_len 00111 const WorkType filter[], 00112 double filter_width) { 00113 // If we are expanding the row (scale>1.0), we need to look at a fractional 00114 // granularity. Hence, we scale our filter index by scale. If we are 00115 // compressing (scale<1.0), we don't need to fiddle with the filter index, so 00116 // we leave it at one. 00117 double iscale = max(scale, 1.0); 00118 00119 // Similarly, if we are expanding the row, we want to start the new row at 00120 // the far left edge of the original pixel, not in the center. So we will 00121 // have a non-zero offset. 00122 int offset = (int)cfloor(iscale*0.5); 00123 00124 for (int dest_x=0; dest_x<dest_len; dest_x++) { 00125 double center = (dest_x-offset)/scale; 00126 00127 // left and right are the starting and ending ranges of the radius of 00128 // interest of the filter function. We need to apply the filter to each 00129 // value in this range. 00130 int left = max((int)cfloor(center - filter_width), 0); 00131 int right = min((int)cceil(center + filter_width), source_len-1); 00132 00133 // right_center is the point just to the right of the center. This 00134 // allows us to flip the sign of the offset when we cross the center point. 00135 int right_center = (int)cceil(center); 00136 00137 WorkType net_weight = 0; 00138 WorkType net_value = 0; 00139 00140 int index, source_x; 00141 00142 // This loop is broken into two pieces--the left of center and the right 00143 // of center--so we don't have to incur the overhead of calling fabs() 00144 // each time through the loop. 00145 for (source_x=left; source_x<right_center; source_x++) { 00146 index = (int)(iscale*(center-source_x)); 00147 net_value += filter[index] * source[source_x]; 00148 net_weight += filter[index]; 00149 } 00150 00151 for (; source_x<=right; source_x++) { 00152 index = (int)(iscale*(source_x-center)); 00153 net_value += filter[index] * source[source_x]; 00154 net_weight += filter[index]; 00155 } 00156 00157 if (net_weight>0) { 00158 dest[dest_x] = (StoreType)(net_value / net_weight); 00159 } else { 00160 dest[dest_x] = 0; 00161 } 00162 } 00163 Thread::consider_yield(); 00164 } 00165 00166 00167 // The various filter functions are called before each axis scaling to build 00168 // an kernel array suitable for the given scaling factor. Given a scaling 00169 // ratio of the axis (dest_len / source_len), and a width parameter supplied 00170 // by the user, they must build an array of filter values (described above) 00171 // and also set the radius of interest of the filter function. 00172 00173 // The values of the elements of filter must completely cover the range 00174 // 0..filter_max; the array must have enough elements to include all indices 00175 // corresponding to values in the range -filter_width to filter_width. 00176 00177 typedef void FilterFunction(double scale, double width, 00178 WorkType *&filter, double &filter_width); 00179 00180 static void 00181 box_filter_impl(double scale, double width, 00182 WorkType *&filter, double &filter_width) { 00183 double fscale; 00184 if (scale < 1.0) { 00185 // If we are compressing the image, we want to expand the range of 00186 // the filter function to prevent dropping below the Nyquist rate. 00187 // Hence, we divide by scale. 00188 fscale = 1.0 / scale; 00189 } else { 00190 00191 // If we are expanding the image, we want to increase the granularity 00192 // of the filter function since we will need to access fractional cel 00193 // values. Hence, we multiply by scale. 00194 fscale = scale; 00195 } 00196 filter_width = width; 00197 int actual_width = (int)cceil((filter_width+1) * fscale); 00198 00199 filter = (WorkType *)PANDA_MALLOC_ARRAY(actual_width * sizeof(WorkType)); 00200 00201 for (int i=0; i<actual_width; i++) { 00202 filter[i] = (i<=filter_width*fscale) ? filter_max : 0; 00203 } 00204 } 00205 00206 static void 00207 gaussian_filter_impl(double scale, double width, 00208 WorkType *&filter, double &filter_width) { 00209 double fscale; 00210 if (scale < 1.0) { 00211 // If we are compressing the image, we want to expand the range of 00212 // the filter function to prevent dropping below the Nyquist rate. 00213 // Hence, we divide by scale (to make fscale larger). 00214 fscale = 1.0 / scale; 00215 } else { 00216 00217 // If we are expanding the image, we want to increase the granularity 00218 // of the filter function since we will need to access fractional cel 00219 // values. Hence, we multiply by scale (to make fscale larger). 00220 fscale = scale; 00221 } 00222 double sigma = width/2; 00223 filter_width = 3.0 * sigma; 00224 int actual_width = (int)cceil((filter_width+1) * fscale); 00225 00226 // G(x, y) = (1/(2 pi sigma^2)) * exp( - (x^2 + y^2) / (2 sigma^2)) 00227 00228 // (We can throw away the initial factor, since these weights will all 00229 // be normalized; and we're only computing a 1-dimensional function, 00230 // so we can ignore the y^2.) 00231 00232 filter = (WorkType *)PANDA_MALLOC_ARRAY(actual_width * sizeof(WorkType)); 00233 double div = 2*sigma*sigma; 00234 00235 for (int i=0; i<actual_width; i++) { 00236 double x = i/fscale; 00237 filter[i] = (WorkType)(filter_max * exp(-x*x / div)); 00238 // The highest value of the exp function in this range is always 1.0, 00239 // at index value 0. Thus, we scale the whole range by filter_max, 00240 // to produce a filter in the range [0..filter_max]. 00241 } 00242 } 00243 00244 00245 // We have a function, defined in pnm-image-filter-core.cxx, that will scale 00246 // an image in both X and Y directions for a particular channel, by setting 00247 // up the temporary matrix appropriately and calling the above functions. 00248 00249 // What we really need is a series of such functions, one for each channel, 00250 // and also one to scale by X first, and one to scale by Y first. This sounds 00251 // a lot like a C++ template: we want to compile the same function several 00252 // times to work on slightly different sorts of things each time. However, 00253 // the things we want to vary are the particular member functions of PNMImage 00254 // that we call (e.g. Red(), Green(), etc.), and we can't declare a template 00255 // of member functions, only of types. 00256 00257 // It's doable using templates. It would involve the declaration of 00258 // lots of silly little functor objects. This is much more compact 00259 // and no more difficult to read. 00260 00261 // The function in pnm-image-filter-core.cxx uses macros to access the member 00262 // functions of PNMImage. Hence, we only need to redefine those macros 00263 // with each instance of the function to cause each instance to operate on 00264 // the correct member. 00265 00266 00267 // These instances scale by X first, then by Y. 00268 00269 #define FUNCTION_NAME filter_red_xy 00270 #define ASIZE get_x_size 00271 #define BSIZE get_y_size 00272 #define GETVAL(a, b) get_red(a, b) 00273 #define SETVAL(a, b, v) set_red(a, b, v) 00274 #include "pnm-image-filter-core.cxx" 00275 #undef SETVAL 00276 #undef GETVAL 00277 #undef BSIZE 00278 #undef ASIZE 00279 #undef FUNCTION_NAME 00280 00281 #define FUNCTION_NAME filter_green_xy 00282 #define ASIZE get_x_size 00283 #define BSIZE get_y_size 00284 #define GETVAL(a, b) get_green(a, b) 00285 #define SETVAL(a, b, v) set_green(a, b, v) 00286 #include "pnm-image-filter-core.cxx" 00287 #undef SETVAL 00288 #undef GETVAL 00289 #undef BSIZE 00290 #undef ASIZE 00291 #undef FUNCTION_NAME 00292 00293 #define FUNCTION_NAME filter_blue_xy 00294 #define ASIZE get_x_size 00295 #define BSIZE get_y_size 00296 #define GETVAL(a, b) get_blue(a, b) 00297 #define SETVAL(a, b, v) set_blue(a, b, v) 00298 #include "pnm-image-filter-core.cxx" 00299 #undef SETVAL 00300 #undef GETVAL 00301 #undef BSIZE 00302 #undef ASIZE 00303 #undef FUNCTION_NAME 00304 00305 #define FUNCTION_NAME filter_gray_xy 00306 #define ASIZE get_x_size 00307 #define BSIZE get_y_size 00308 #define GETVAL(a, b) get_bright(a, b) 00309 #define SETVAL(a, b, v) set_xel(a, b, v) 00310 #include "pnm-image-filter-core.cxx" 00311 #undef SETVAL 00312 #undef GETVAL 00313 #undef BSIZE 00314 #undef ASIZE 00315 #undef FUNCTION_NAME 00316 00317 #define FUNCTION_NAME filter_alpha_xy 00318 #define ASIZE get_x_size 00319 #define BSIZE get_y_size 00320 #define GETVAL(a, b) get_alpha(a, b) 00321 #define SETVAL(a, b, v) set_alpha(a, b, v) 00322 #include "pnm-image-filter-core.cxx" 00323 #undef SETVAL 00324 #undef GETVAL 00325 #undef BSIZE 00326 #undef ASIZE 00327 #undef FUNCTION_NAME 00328 00329 00330 // These instances scale by Y first, then by X. 00331 00332 #define FUNCTION_NAME filter_red_yx 00333 #define ASIZE get_y_size 00334 #define BSIZE get_x_size 00335 #define GETVAL(a, b) get_red(b, a) 00336 #define SETVAL(a, b, v) set_red(b, a, v) 00337 #include "pnm-image-filter-core.cxx" 00338 #undef SETVAL 00339 #undef GETVAL 00340 #undef BSIZE 00341 #undef ASIZE 00342 #undef FUNCTION_NAME 00343 00344 #define FUNCTION_NAME filter_green_yx 00345 #define ASIZE get_y_size 00346 #define BSIZE get_x_size 00347 #define GETVAL(a, b) get_green(b, a) 00348 #define SETVAL(a, b, v) set_green(b, a, v) 00349 #include "pnm-image-filter-core.cxx" 00350 #undef SETVAL 00351 #undef GETVAL 00352 #undef BSIZE 00353 #undef ASIZE 00354 #undef FUNCTION_NAME 00355 00356 #define FUNCTION_NAME filter_blue_yx 00357 #define ASIZE get_y_size 00358 #define BSIZE get_x_size 00359 #define GETVAL(a, b) get_blue(b, a) 00360 #define SETVAL(a, b, v) set_blue(b, a, v) 00361 #include "pnm-image-filter-core.cxx" 00362 #undef SETVAL 00363 #undef GETVAL 00364 #undef BSIZE 00365 #undef ASIZE 00366 #undef FUNCTION_NAME 00367 00368 #define FUNCTION_NAME filter_gray_yx 00369 #define ASIZE get_y_size 00370 #define BSIZE get_x_size 00371 #define GETVAL(a, b) get_bright(b, a) 00372 #define SETVAL(a, b, v) set_xel(b, a, v) 00373 #include "pnm-image-filter-core.cxx" 00374 #undef SETVAL 00375 #undef GETVAL 00376 #undef BSIZE 00377 #undef ASIZE 00378 #undef FUNCTION_NAME 00379 00380 #define FUNCTION_NAME filter_alpha_yx 00381 #define ASIZE get_y_size 00382 #define BSIZE get_x_size 00383 #define GETVAL(a, b) get_alpha(b, a) 00384 #define SETVAL(a, b, v) set_alpha(b, a, v) 00385 #include "pnm-image-filter-core.cxx" 00386 #undef SETVAL 00387 #undef GETVAL 00388 #undef BSIZE 00389 #undef ASIZE 00390 #undef FUNCTION_NAME 00391 00392 00393 // filter_image pulls everything together, and filters one image into 00394 // another. Both images can be the same with no ill effects. 00395 static void 00396 filter_image(PNMImage &dest, const PNMImage &source, 00397 double width, FilterFunction *make_filter) { 00398 00399 // We want to scale by the smallest destination axis first, for a 00400 // slight performance gain. 00401 00402 if (dest.get_x_size() <= dest.get_y_size()) { 00403 if (dest.is_grayscale() || source.is_grayscale()) { 00404 filter_gray_xy(dest, source, width, make_filter); 00405 } else { 00406 filter_red_xy(dest, source, width, make_filter); 00407 filter_green_xy(dest, source, width, make_filter); 00408 filter_blue_xy(dest, source, width, make_filter); 00409 } 00410 00411 if (dest.has_alpha() && source.has_alpha()) { 00412 filter_alpha_xy(dest, source, width, make_filter); 00413 } 00414 00415 } else { 00416 if (dest.is_grayscale() || source.is_grayscale()) { 00417 filter_gray_yx(dest, source, width, make_filter); 00418 } else { 00419 filter_red_yx(dest, source, width, make_filter); 00420 filter_green_yx(dest, source, width, make_filter); 00421 filter_blue_yx(dest, source, width, make_filter); 00422 } 00423 00424 if (dest.has_alpha() && source.has_alpha()) { 00425 filter_alpha_yx(dest, source, width, make_filter); 00426 } 00427 } 00428 } 00429 00430 00431 00432 //////////////////////////////////////////////////////////////////// 00433 // Function: PNMImage::box_filter_from 00434 // Access: Public 00435 // Description: Makes a resized copy of the indicated image into this 00436 // one using the indicated filter. The image to be 00437 // copied is squashed and stretched to match the 00438 // dimensions of the current image, applying the 00439 // appropriate filter to perform the stretching. 00440 //////////////////////////////////////////////////////////////////// 00441 void PNMImage:: 00442 box_filter_from(double width, const PNMImage ©) { 00443 filter_image(*this, copy, width, &box_filter_impl); 00444 } 00445 00446 //////////////////////////////////////////////////////////////////// 00447 // Function: PNMImage::gaussian_filter_from 00448 // Access: Public 00449 // Description: Makes a resized copy of the indicated image into this 00450 // one using the indicated filter. The image to be 00451 // copied is squashed and stretched to match the 00452 // dimensions of the current image, applying the 00453 // appropriate filter to perform the stretching. 00454 //////////////////////////////////////////////////////////////////// 00455 void PNMImage:: 00456 gaussian_filter_from(double width, const PNMImage ©) { 00457 filter_image(*this, copy, width, &gaussian_filter_impl); 00458 } 00459 00460 00461 // 00462 // The following functions are support for quick_box_filter(). 00463 // 00464 00465 INLINE void 00466 box_filter_xel(const PNMImage &image, 00467 int x, int y, double x_contrib, double y_contrib, 00468 double &red, double &grn, double &blu, double &alpha, 00469 double &pixel_count) { 00470 double contrib = x_contrib * y_contrib; 00471 red += image.get_red_val(x, y) * contrib; 00472 grn += image.get_green_val(x, y) * contrib; 00473 blu += image.get_blue_val(x, y) * contrib; 00474 if (image.has_alpha()) { 00475 alpha += image.get_alpha_val(x, y) * contrib; 00476 } 00477 00478 pixel_count += contrib; 00479 } 00480 00481 00482 INLINE void 00483 box_filter_line(const PNMImage &image, 00484 double x0, int y, double x1, double y_contrib, 00485 double &red, double &grn, double &blu, double &alpha, 00486 double &pixel_count) { 00487 int x = (int)x0; 00488 // Get the first (partial) xel 00489 box_filter_xel(image, x, y, (double)(x+1)-x0, y_contrib, 00490 red, grn, blu, alpha, pixel_count); 00491 00492 int x_last = (int)x1; 00493 if (x < x_last) { 00494 x++; 00495 while (x < x_last) { 00496 // Get each consecutive (complete) xel 00497 box_filter_xel(image, x, y, 1.0, y_contrib, 00498 red, grn, blu, alpha, pixel_count); 00499 x++; 00500 } 00501 00502 // Get the final (partial) xel 00503 double x_contrib = x1 - (double)x_last; 00504 if (x_contrib > 0.0001) { 00505 box_filter_xel(image, x, y, x_contrib, y_contrib, 00506 red, grn, blu, alpha, pixel_count); 00507 } 00508 } 00509 } 00510 00511 static void 00512 box_filter_region(const PNMImage &image, 00513 double x0, double y0, double x1, double y1, 00514 xel &result, xelval &alpha_result) { 00515 double red = 0.0, grn = 0.0, blu = 0.0, alpha = 0.0; 00516 double pixel_count = 0.0; 00517 00518 assert(y0 >=0 && y1 >=0); 00519 00520 int y = (int)y0; 00521 // Get the first (partial) row 00522 box_filter_line(image, x0, y, x1, (double)(y+1)-y0, 00523 red, grn, blu, alpha, pixel_count); 00524 00525 int y_last = (int)y1; 00526 if (y < y_last) { 00527 y++; 00528 while (y < y_last) { 00529 // Get each consecutive (complete) row 00530 box_filter_line(image, x0, y, x1, 1.0, 00531 red, grn, blu, alpha, pixel_count); 00532 y++; 00533 } 00534 00535 // Get the final (partial) row 00536 double y_contrib = y1 - (double)y_last; 00537 if (y_contrib > 0.0001) { 00538 box_filter_line(image, x0, y, x1, y_contrib, 00539 red, grn, blu, alpha, pixel_count); 00540 } 00541 } 00542 00543 PPM_ASSIGN(result, 00544 (xelval)(red / pixel_count + 0.5), 00545 (xelval)(grn / pixel_count + 0.5), 00546 (xelval)(blu / pixel_count + 0.5)); 00547 00548 alpha_result = (xelval)(alpha / pixel_count + 0.5); 00549 } 00550 00551 //////////////////////////////////////////////////////////////////// 00552 // Function: PNMImage::quick_filter_from 00553 // Access: Public 00554 // Description: Resizes from the given image, with a fixed radius of 00555 // 0.5. This is a very specialized and simple algorithm 00556 // that doesn't handle dropping below the Nyquist rate 00557 // very well, but is quite a bit faster than the more 00558 // general box_filter(), above. If borders are 00559 // specified, they will further restrict the size of the 00560 // resulting image. There's no point in using 00561 // quick_box_filter() on a single image. 00562 //////////////////////////////////////////////////////////////////// 00563 void PNMImage:: 00564 quick_filter_from(const PNMImage &from, int xborder, int yborder) { 00565 int from_xs = from.get_x_size(); 00566 int from_ys = from.get_y_size(); 00567 00568 int to_xs = get_x_size() - xborder; 00569 int to_ys = get_y_size() - yborder; 00570 00571 int to_xoff = xborder / 2; 00572 int to_yoff = yborder / 2; 00573 00574 double from_x0, from_x1, from_y0, from_y1; 00575 int to_x, to_y; 00576 00577 double x_scale = (double)from_xs / (double)to_xs; 00578 double y_scale = (double)from_ys / (double)to_ys; 00579 00580 from_y0 = max(0, -to_yoff) * y_scale; 00581 for (to_y = max(0, -to_yoff); 00582 to_y < min(to_ys, get_y_size()-to_yoff); 00583 to_y++) { 00584 from_y1 = (to_y+1) * y_scale; 00585 00586 from_x0 = max(0, -to_xoff) * x_scale; 00587 for (to_x = max(0, -to_xoff); 00588 to_x < min(to_xs, get_x_size()-to_xoff); 00589 to_x++) { 00590 from_x1 = (to_x+1) * x_scale; 00591 00592 // Now the box from (from_x0, from_y0) - (from_x1, from_y1) 00593 // but not including (from_x1, from_y1) maps to the pixel (to_x, to_y). 00594 xelval alpha_result; 00595 box_filter_region(from, 00596 from_x0, from_y0, from_x1, from_y1, 00597 (*this)[to_yoff + to_y][to_xoff + to_x], 00598 alpha_result); 00599 if (has_alpha()) { 00600 set_alpha_val(to_xoff+to_x, to_yoff+to_y, alpha_result); 00601 } 00602 00603 from_x0 = from_x1; 00604 } 00605 from_y0 = from_y1; 00606 Thread::consider_yield(); 00607 } 00608 }