Panda3D

pnm-image-filter.cxx

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 &copy) {
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 &copy) {
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 }
 All Classes Functions Variables Enumerations