PfmFile::box_filter() etc., improved forward_distort().

panda/src/pnmimage/pfmFile.I

@@ -144,12 +144,16 @@ set_point2(int x, int y, const LVecBase2f &point) {
     _table[(y * _x_size + x)] = point[0];
     break;
 
+  case 2:
+    *(LPoint2f *)&_table[(y * _x_size + x) * _num_channels] = point;
+    break;
+
   case 3:
-    (*(LPoint3f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1], point[2]);
+    (*(LPoint3f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1], 0.0);
     break;
 
-  case 2:
-    *(LPoint2f *)&_table[(y * _x_size + x) * _num_channels] = point;
+  case 4:
+    (*(LPoint4f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1], 0.0, 0.0);
     break;
   } 
 }
@@ -214,10 +218,17 @@ set_point(int x, int y, const LVecBase3f &point) {
     _table[(y * _x_size + x)] = point[0];
     break;
 
+  case 2:
+    (*(LPoint2f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1]);
+    break;
+
   case 3:
-  case 4:
     *(LPoint3f *)&_table[(y * _x_size + x) * _num_channels] = point;
     break;
+
+  case 4:
+    (*(LPoint4f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1], 0.0f, 0.0f);
+    break;
   } 
 }
 
@@ -281,6 +292,10 @@ set_point4(int x, int y, const LVecBase4f &point) {
     _table[(y * _x_size + x)] = point[0];
     break;
 
+  case 2:
+    (*(LPoint2f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1]);
+    break;
+
   case 3:
     (*(LPoint3f *)&_table[(y * _x_size + x) * _num_channels]).set(point[0], point[1], point[2]);
     break;
@@ -320,6 +335,36 @@ modify_point4(int x, int y) {
   return *(LPoint4f *)&_table[(y * _x_size + x) * _num_channels];
 }
 
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::fill
+//       Access: Published
+//  Description: Fills the table with all of the same value.
+////////////////////////////////////////////////////////////////////
+INLINE void PfmFile::
+fill(PN_float32 value) {
+  fill(LPoint4f(value, 0.0f, 0.0f, 0.0f));
+}
+
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::fill
+//       Access: Published
+//  Description: Fills the table with all of the same value.
+////////////////////////////////////////////////////////////////////
+INLINE void PfmFile::
+fill(const LPoint2f &value) {
+  fill(LPoint4f(value[0], value[1], 0.0f, 0.0f));
+}
+
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::fill
+//       Access: Published
+//  Description: Fills the table with all of the same value.
+////////////////////////////////////////////////////////////////////
+INLINE void PfmFile::
+fill(const LPoint3f &value) {
+  fill(LPoint4f(value[0], value[1], value[2], 0.0f));
+}
+
 ////////////////////////////////////////////////////////////////////
 //     Function: PfmFile::calc_autocrop
 //       Access: Published

panda/src/pnmimage/pfmFile.cxx

@@ -470,6 +470,56 @@ store_mask(PNMImage &pnmimage) const {
   return true;
 }
 
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::fill
+//       Access: Published
+//  Description: Fills the table with all of the same value.
+////////////////////////////////////////////////////////////////////
+void PfmFile::
+fill(const LPoint4f &value) {
+  switch (_num_channels) {
+  case 1:
+    {
+      for (int yi = 0; yi < _y_size; ++yi) {
+        for (int xi = 0; xi < _x_size; ++xi) {
+          _table[(yi * _x_size + xi)] = value[0];
+        }
+      }
+    }
+    break;
+
+  case 2:
+    {
+      for (int yi = 0; yi < _y_size; ++yi) {
+        for (int xi = 0; xi < _x_size; ++xi) {
+          (*(LPoint2f *)&_table[(yi * _x_size + xi) * _num_channels]).set(value[0], value[1]);
+        }
+      }
+    }
+    break;
+
+  case 3:
+    {
+      for (int yi = 0; yi < _y_size; ++yi) {
+        for (int xi = 0; xi < _x_size; ++xi) {
+          (*(LPoint3f *)&_table[(yi * _x_size + xi) * _num_channels]).set(value[0], value[1], value[2]);
+        }
+      }
+    }
+    break;
+
+  case 4:
+    {
+      for (int yi = 0; yi < _y_size; ++yi) {
+        for (int xi = 0; xi < _x_size; ++xi) {
+          *(LPoint4f *)&_table[(yi * _x_size + xi) * _num_channels] = value;
+        }
+      }
+    }
+    break;
+  } 
+}
+
 ////////////////////////////////////////////////////////////////////
 //     Function: PfmFile::calc_average_point
 //       Access: Published
@@ -736,41 +786,73 @@ resize(int new_x_size, int new_y_size) {
     return;
   }
 
-  int new_size = new_x_size * new_y_size * _num_channels;
+  PfmFile result;
+  result.clear(new_x_size, new_y_size, _num_channels);
+
+  if (new_x_size < _x_size && new_y_size < _y_size) {
+    // If we're downscaling, we can use quick_filter, which is faster.
+    result.quick_filter_from(*this);
+
+  } else {
+    // Otherwise, we should use box_filter(), which is more general.
+    result.box_filter_from(0.5, *this);
+  }
+
+  _table.swap(result._table);
+  _x_size = new_x_size;
+  _y_size = new_y_size;
+}
+
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::quick_filter_from
+//       Access: Public
+//  Description: Resizes from the given image, with a fixed radius of
+//               0.5. This is a very specialized and simple algorithm
+//               that doesn't handle dropping below the Nyquist rate
+//               very well, but is quite a bit faster than the more
+//               general box_filter(), above.
+////////////////////////////////////////////////////////////////////
+void PfmFile::
+quick_filter_from(const PfmFile &from) {
+  if (_x_size == 0 || _y_size == 0) {
+    return;
+  }
 
-  // We allocate a little bit bigger to allow safe overflow.
   Table new_data;
-  new_data.reserve(new_size + 4);
+  new_data.reserve(_table.size());
 
   PN_float32 from_x0, from_x1, from_y0, from_y1;
 
+  int orig_x_size = from.get_x_size();
+  int orig_y_size = from.get_y_size();
+
   PN_float32 x_scale = 1.0;
   PN_float32 y_scale = 1.0;
 
-  if (new_x_size > 1) {
-    x_scale = (PN_float32)_x_size / (PN_float32)new_x_size;
+  if (_x_size > 1) {
+    x_scale = (PN_float32)orig_x_size / (PN_float32)_x_size;
   }
-  if (new_y_size > 1) {
-    y_scale = (PN_float32)_y_size / (PN_float32)new_y_size;
+  if (_y_size > 1) {
+    y_scale = (PN_float32)orig_y_size / (PN_float32)_y_size;
   }
 
   switch (_num_channels) {
   case 1:
     {
       from_y0 = 0.0;
-      for (int to_y = 0; to_y < new_y_size; ++to_y) {
+      for (int to_y = 0; to_y < _y_size; ++to_y) {
         from_y1 = (to_y + 1.0) * y_scale;
-        from_y1 = min(from_y1, (PN_float32) _y_size);
+        from_y1 = min(from_y1, (PN_float32)orig_y_size);
         
         from_x0 = 0.0;
-        for (int to_x = 0; to_x < new_x_size; ++to_x) {
+        for (int to_x = 0; to_x < _x_size; ++to_x) {
           from_x1 = (to_x + 1.0) * x_scale;
-          from_x1 = min(from_x1, (PN_float32) _x_size);
+          from_x1 = min(from_x1, (PN_float32)orig_x_size);
           
           // Now the box from (from_x0, from_y0) - (from_x1, from_y1)
           // but not including (from_x1, from_y1) maps to the pixel (to_x, to_y).
           PN_float32 result;
-          box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
+          from.box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
           new_data.push_back(result);
           
           from_x0 = from_x1;
@@ -783,19 +865,19 @@ resize(int new_x_size, int new_y_size) {
   case 3:
     {
       from_y0 = 0.0;
-      for (int to_y = 0; to_y < new_y_size; ++to_y) {
+      for (int to_y = 0; to_y < _y_size; ++to_y) {
         from_y1 = (to_y + 1.0) * y_scale;
-        from_y1 = min(from_y1, (PN_float32) _y_size);
+        from_y1 = min(from_y1, (PN_float32)orig_y_size);
         
         from_x0 = 0.0;
-        for (int to_x = 0; to_x < new_x_size; ++to_x) {
+        for (int to_x = 0; to_x < _x_size; ++to_x) {
           from_x1 = (to_x + 1.0) * x_scale;
-          from_x1 = min(from_x1, (PN_float32) _x_size);
+          from_x1 = min(from_x1, (PN_float32)orig_x_size);
           
           // Now the box from (from_x0, from_y0) - (from_x1, from_y1)
           // but not including (from_x1, from_y1) maps to the pixel (to_x, to_y).
           LPoint3f result;
-          box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
+          from.box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
           new_data.push_back(result[0]);
           new_data.push_back(result[1]);
           new_data.push_back(result[2]);
@@ -810,19 +892,19 @@ resize(int new_x_size, int new_y_size) {
   case 4:
     {
       from_y0 = 0.0;
-      for (int to_y = 0; to_y < new_y_size; ++to_y) {
+      for (int to_y = 0; to_y < _y_size; ++to_y) {
         from_y1 = (to_y + 1.0) * y_scale;
-        from_y1 = min(from_y1, (PN_float32) _y_size);
+        from_y1 = min(from_y1, (PN_float32)orig_y_size);
         
         from_x0 = 0.0;
-        for (int to_x = 0; to_x < new_x_size; ++to_x) {
+        for (int to_x = 0; to_x < _x_size; ++to_x) {
           from_x1 = (to_x + 1.0) * x_scale;
-          from_x1 = min(from_x1, (PN_float32) _x_size);
+          from_x1 = min(from_x1, (PN_float32)orig_x_size);
           
           // Now the box from (from_x0, from_y0) - (from_x1, from_y1)
           // but not including (from_x1, from_y1) maps to the pixel (to_x, to_y).
           LPoint4f result;
-          box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
+          from.box_filter_region(result, from_x0, from_y0, from_x1, from_y1);
           new_data.push_back(result[0]);
           new_data.push_back(result[1]);
           new_data.push_back(result[2]);
@@ -844,10 +926,8 @@ resize(int new_x_size, int new_y_size) {
   new_data.push_back(0.0);
   new_data.push_back(0.0);
 
-  nassertv(new_data.size() == new_size + 4);
+  nassertv(new_data.size() == _table.size());
   _table.swap(new_data);
-  _x_size = new_x_size;
-  _y_size = new_y_size;
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -963,53 +1043,72 @@ xform(const LMatrix4f &transform) {
 //               current file is replaced with the point from the same
 //               file at (x,y), where (x,y) is the point value from
 //               the dist map.
+//
+//               If scale_factor is not 1, it should be a value > 1,
+//               and it specifies the factor to upscale the working
+//               table while processing, to reduce artifacts from
+//               integer truncation.
 ////////////////////////////////////////////////////////////////////
 void PfmFile::
-forward_distort(const PfmFile &dist) {
-  PfmFile result;
+forward_distort(const PfmFile &dist, PN_float32 scale_factor) {
+  int working_x_size = (int)cceil(_x_size * scale_factor);
+  int working_y_size = (int)cceil(_y_size * scale_factor);
+
+  working_x_size = max(working_x_size, dist.get_x_size());
+  working_y_size = max(working_y_size, dist.get_y_size());
+
+  const PfmFile *source_p = this;
+  PfmFile scaled_source;
+  if ((*this).get_x_size() != working_x_size || (*this).get_y_size() != working_y_size) {
+    // Rescale the source file as needed.
+    scaled_source = (*this);
+    scaled_source.resize(working_x_size, working_y_size);
+    source_p = &scaled_source;
+  }
+
   const PfmFile *dist_p = &dist;
   PfmFile scaled_dist;
-  if (dist.get_x_size() != _x_size || dist.get_y_size() != _y_size) {
+  if (dist.get_x_size() != working_x_size || dist.get_y_size() != working_y_size) {
     // Rescale the dist file as needed.
     scaled_dist = dist;
-    scaled_dist.resize(_x_size, _y_size);
+    scaled_dist.resize(working_x_size, working_y_size);
     dist_p = &scaled_dist;
   }
 
+  PfmFile result;
+  result.clear(working_x_size, working_y_size, _num_channels);
+    
   if (_has_no_data_value) {
-    result = *this;
-    for (int yi = 0; yi < _y_size; ++yi) {
-      for (int xi = 0; xi < _x_size; ++xi) {
-        result.set_point4(xi, yi, _no_data_value);
-      }
-    }
-
-  } else {
-    result.clear(_x_size, _y_size, _num_channels);
+    result.set_no_data_value(_no_data_value);
+    result.fill(_no_data_value);
   }
 
-  for (int yi = 0; yi < _y_size; ++yi) {
-    for (int xi = 0; xi < _x_size; ++xi) {
+  for (int yi = 0; yi < working_y_size; ++yi) {
+    for (int xi = 0; xi < working_x_size; ++xi) {
       if (!dist_p->has_point(xi, yi)) {
         continue;
       }
       LPoint2f uv = dist_p->get_point2(xi, yi);
-      int dist_xi = (int)cfloor(uv[0] * (double)_x_size);
-      int dist_yi = (int)cfloor(uv[1] * (double)_y_size);
-      if (dist_xi < 0 || dist_xi >= _x_size || 
-          dist_yi < 0 || dist_yi >= _y_size) {
+      int dist_xi = (int)cfloor(uv[0] * (PN_float32)working_x_size);
+      int dist_yi = (int)cfloor(uv[1] * (PN_float32)working_y_size);
+      if (dist_xi < 0 || dist_xi >= working_x_size || 
+          dist_yi < 0 || dist_yi >= working_y_size) {
         continue;
       }
 
-      if (!has_point(dist_xi, dist_yi)) {
+      if (!source_p->has_point(dist_xi, dist_yi)) {
         continue;
       }
 
-      result.set_point(xi, yi, get_point(dist_xi, dist_yi));
+      result.set_point(xi, yi, source_p->get_point(dist_xi, dist_yi));
     }
   }
 
-  (*this) = result;
+  // Resize to the target size for completion.
+  result.resize(_x_size, _y_size);
+
+  nassertv(result._table.size() == _table.size());
+  _table.swap(result._table);
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -1022,51 +1121,75 @@ forward_distort(const PfmFile &dist) {
 //               current file is replaced with the point from the same
 //               file at (u,v), where (x,y) is the point value from
 //               the dist map.
+//
+//               If scale_factor is not 1, it should be a value > 1,
+//               and it specifies the factor to upscale the working
+//               table while processing, to reduce artifacts from
+//               integer truncation.
 ////////////////////////////////////////////////////////////////////
 void PfmFile::
-reverse_distort(const PfmFile &dist) {
-  PfmFile result;
+reverse_distort(const PfmFile &dist, PN_float32 scale_factor) {
+  int working_x_size = (int)cceil(_x_size * scale_factor);
+  int working_y_size = (int)cceil(_y_size * scale_factor);
+
+  working_x_size = max(working_x_size, dist.get_x_size());
+  working_y_size = max(working_y_size, dist.get_y_size());
+
+  const PfmFile *source_p = this;
+  PfmFile scaled_source;
+  if ((*this).get_x_size() != working_x_size || (*this).get_y_size() != working_y_size) {
+    // Rescale the source file as needed.
+    scaled_source = (*this);
+    scaled_source.resize(working_x_size, working_y_size);
+    source_p = &scaled_source;
+  }
+
   const PfmFile *dist_p = &dist;
   PfmFile scaled_dist;
-  if (dist.get_x_size() != _x_size || dist.get_y_size() != _y_size) {
+  if (dist.get_x_size() != working_x_size || dist.get_y_size() != working_y_size) {
     // Rescale the dist file as needed.
     scaled_dist = dist;
-    scaled_dist.resize(_x_size, _y_size);
+    scaled_dist.resize(working_x_size, working_y_size);
     dist_p = &scaled_dist;
   }
 
+  PfmFile result;
+  result.clear(working_x_size, working_y_size, _num_channels);
+    
   if (_has_no_data_value) {
-    result = *this;
-    for (int yi = 0; yi < _y_size; ++yi) {
-      for (int xi = 0; xi < _x_size; ++xi) {
-        result.set_point4(xi, yi, _no_data_value);
-      }
-    }
-  } else {
-    result.clear(_x_size, _y_size, _num_channels);
+    result.set_no_data_value(_no_data_value);
+    result.fill(_no_data_value);
   }
 
-  for (int yi = 0; yi < _y_size; ++yi) {
-    for (int xi = 0; xi < _x_size; ++xi) {
-      if (!has_point(xi, yi)) {
+  for (int yi = 0; yi < working_y_size; ++yi) {
+    for (int xi = 0; xi < working_x_size; ++xi) {
+      if (!source_p->has_point(xi, yi)) {
         continue;
       }
       if (!dist_p->has_point(xi, yi)) {
         continue;
       }
       LPoint2f uv = dist_p->get_point2(xi, yi);
-      int dist_xi = (int)cfloor(uv[0] * (double)_x_size);
-      int dist_yi = (int)cfloor(uv[1] * (double)_y_size);
-      if (dist_xi < 0 || dist_xi >= _x_size || 
-          dist_yi < 0 || dist_yi >= _y_size) {
+      int dist_xi = (int)cfloor(uv[0] * (PN_float32)working_x_size);
+      int dist_yi = (int)cfloor(uv[1] * (PN_float32)working_y_size);
+      if (dist_xi < 0 || dist_xi >= working_x_size || 
+          dist_yi < 0 || dist_yi >= working_y_size) {
+        continue;
+      }
+
+      if (!source_p->has_point(dist_xi, dist_yi)) {
         continue;
       }
 
-      result.set_point(dist_xi, dist_yi, get_point(xi, yi));
+      result.set_point(dist_xi, dist_yi, source_p->get_point(xi, yi));
     }
   }
 
-  (*this) = result;
+  // Resize to the target size for completion.
+  result.resize(_x_size, _y_size);
+
+  nassertv(result._table.size() == _table.size());
+  _table.swap(result._table);
 }
 
 ////////////////////////////////////////////////////////////////////
@@ -1089,7 +1212,7 @@ merge(const PfmFile &other) {
 
   for (int yi = 0; yi < _y_size; ++yi) {
     for (int xi = 0; xi < _x_size; ++xi) {
-      if (!has_point(xi, yi)) {
+      if (!has_point(xi, yi) && other.has_point(xi, yi)) {
         set_point(xi, yi, other.get_point(xi, yi));
       }
     }
@@ -1180,6 +1303,47 @@ clear_to_texcoords(int x_size, int y_size) {
   }
 }
 
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::calc_tight_bounds
+//       Access: Published
+//  Description: Calculates the minimum and maximum vertices of all
+//               points within the table.  Assumes the table contains
+//               3-D points.
+//
+//               The return value is true if any points in the table,
+//               or false if none are.
+////////////////////////////////////////////////////////////////////
+bool PfmFile::
+calc_tight_bounds(LPoint3f &min_point, LPoint3f &max_point) const {
+  min_point.set(0.0f, 0.0f, 0.0f);
+  max_point.set(0.0f, 0.0f, 0.0f);
+
+  bool found_any = false;
+  for (int yi = 0; yi < _y_size; ++yi) {
+    for (int xi = 0; xi < _x_size; ++xi) {
+      if (!has_point(xi, yi)) {
+        continue;
+      }
+   
+      const LPoint3f &point = get_point(xi, yi);
+      if (!found_any) {
+        min_point = point;
+        max_point = point;
+        found_any = true;
+      } else {
+        min_point.set(min(min_point[0], point[0]),
+                      min(min_point[0], point[0]),
+                      min(min_point[0], point[0]));
+        max_point.set(max(max_point[0], point[0]),
+                      max(max_point[0], point[0]),
+                      max(max_point[0], point[0]));
+      }
+    }
+  }
+
+  return found_any;
+}
+
 ////////////////////////////////////////////////////////////////////
 //     Function: PfmFile::compute_planar_bounds
 //       Access: Published

panda/src/pnmimage/pfmFile.h

@@ -74,6 +74,11 @@ PUBLISHED:
   INLINE void set_point4(int x, int y, const LVecBase4d &point);
   INLINE LPoint4f &modify_point4(int x, int y);
 
+  INLINE void fill(PN_float32 value);
+  INLINE void fill(const LPoint2f &value);
+  INLINE void fill(const LPoint3f &value);
+  void fill(const LPoint4f &value);
+
   BLOCKING bool calc_average_point(LPoint3f &result, PN_float32 x, PN_float32 y, PN_float32 radius) const;
   BLOCKING bool calc_min_max(LVecBase3f &min_points, LVecBase3f &max_points) const;
   BLOCKING bool calc_autocrop(int &x_begin, int &x_end, int &y_begin, int &y_end) const;
@@ -92,17 +97,22 @@ PUBLISHED:
   INLINE const LPoint4f &get_no_data_value() const;
 
   BLOCKING void resize(int new_x_size, int new_y_size);
+  BLOCKING void box_filter_from(double radius, const PfmFile &copy);
+  BLOCKING void gaussian_filter_from(double radius, const PfmFile &copy);
+  BLOCKING void quick_filter_from(const PfmFile &copy);
+
   BLOCKING void reverse_rows();
   BLOCKING void flip(bool flip_x, bool flip_y, bool transpose);
   BLOCKING void xform(const LMatrix4f &transform);
   INLINE BLOCKING void xform(const LMatrix4d &transform);
-  BLOCKING void forward_distort(const PfmFile &dist);
-  BLOCKING void reverse_distort(const PfmFile &dist);
+  BLOCKING void forward_distort(const PfmFile &dist, PN_float32 scale_factor = 1.0);
+  BLOCKING void reverse_distort(const PfmFile &dist, PN_float32 scale_factor = 1.0);
   BLOCKING void merge(const PfmFile &other);
   BLOCKING void copy_channel(int to_channel, const PfmFile &other, int from_channel);
   BLOCKING void apply_crop(int x_begin, int x_end, int y_begin, int y_end);
   BLOCKING void clear_to_texcoords(int x_size, int y_size);
 
+  bool calc_tight_bounds(LPoint3f &min_point, LPoint3f &max_point) const;
   BLOCKING PT(BoundingHexahedron) compute_planar_bounds(const LPoint2f &center, PN_float32 point_dist, PN_float32 sample_radius, bool points_only) const;
   INLINE BLOCKING PT(BoundingHexahedron) compute_planar_bounds(const LPoint2d &center, PN_float32 point_dist, PN_float32 sample_radius, bool points_only) const;
   void compute_sample_point(LPoint3f &result,

panda/src/pnmimage/pnm-image-filter-core.cxx

@@ -18,8 +18,8 @@
 
 
 static void
-FUNCTION_NAME(PNMImage &dest, const PNMImage &source,
-              double width, FilterFunction *make_filter) {
+FUNCTION_NAME(IMAGETYPE &dest, const IMAGETYPE &source,
+              double width, FilterFunction *make_filter, int channel) {
   if (!dest.is_valid() || !source.is_valid()) {
     return;
   }
@@ -37,20 +37,22 @@ FUNCTION_NAME(PNMImage &dest, const PNMImage &source,
     matrix[a] = (StoreType *)PANDA_MALLOC_ARRAY(source.BSIZE() * sizeof(StoreType));
   }
 
-  // Now, scale the image in the A direction.
-  double scale = (double)dest.ASIZE() / (double)source.ASIZE();
+  // First, scale the image in the A direction.
+  double scale;
+  StoreType *temp_source, *temp_dest;
 
-  StoreType *temp_source = (StoreType *)PANDA_MALLOC_ARRAY(source.ASIZE() * sizeof(StoreType));
-  StoreType *temp_dest = (StoreType *)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType));
+  scale = (double)dest.ASIZE() / (double)source.ASIZE();
+  temp_source = (StoreType *)PANDA_MALLOC_ARRAY(source.ASIZE() * sizeof(StoreType));
+  temp_dest = (StoreType *)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType));
 
   WorkType *filter;
   double filter_width;
 
   make_filter(scale, width, filter, filter_width);
 
-  for (b=0; b<source.BSIZE(); b++) {
-    for (a=0; a<source.ASIZE(); a++) {
-      temp_source[a] = (StoreType)(source_max * source.GETVAL(a, b));
+  for (b = 0; b < source.BSIZE(); b++) {
+    for (a = 0; a < source.ASIZE(); a++) {
+      temp_source[a] = (StoreType)(source_max * source.GETVAL(a, b, channel));
     }
 
     filter_row(temp_dest, dest.ASIZE(),
@@ -58,7 +60,7 @@ FUNCTION_NAME(PNMImage &dest, const PNMImage &source,
                scale,
                filter, filter_width);
 
-    for (a=0; a<dest.ASIZE(); a++) {
+    for (a = 0; a < dest.ASIZE(); a++) {
       matrix[a][b] = temp_dest[a];
     }
   }
@@ -69,20 +71,18 @@ FUNCTION_NAME(PNMImage &dest, const PNMImage &source,
 
   // Now, scale the image in the B direction.
   scale = (double)dest.BSIZE() / (double)source.BSIZE();
-
   temp_dest = (StoreType *)PANDA_MALLOC_ARRAY(dest.BSIZE() * sizeof(StoreType));
 
   make_filter(scale, width, filter, filter_width);
 
-  for (a=0; a<dest.ASIZE(); a++) {
-
+  for (a = 0; a < dest.ASIZE(); a++) {
     filter_row(temp_dest, dest.BSIZE(),
                matrix[a], source.BSIZE(),
                scale,
                filter, filter_width);
 
-    for (b=0; b<dest.BSIZE(); b++) {
-      dest.SETVAL(a, b, (double)temp_dest[b]/(double)source_max);
+    for (b = 0; b < dest.BSIZE(); b++) {
+      dest.SETVAL(a, b, channel, (double)temp_dest[b]/(double)source_max);
     }
   }
 
@@ -91,7 +91,7 @@ FUNCTION_NAME(PNMImage &dest, const PNMImage &source,
 
   // Now, clean up our temp matrix and go home!
 
-  for (a=0; a<dest.ASIZE(); a++) {
+  for (a = 0; a < dest.ASIZE(); a++) {
     PANDA_FREE_ARRAY(matrix[a]);
   }
   PANDA_FREE_ARRAY(matrix);

panda/src/pnmimage/pnm-image-filter-sparse-core.cxx

@@ -0,0 +1,114 @@
+// Filename: pnm-image-filter-sparse-core.cxx
+// Created by:  drose (25Jan13)
+//
+////////////////////////////////////////////////////////////////////
+//
+// PANDA 3D SOFTWARE
+// Copyright (c) Carnegie Mellon University.  All rights reserved.
+//
+// All use of this software is subject to the terms of the revised BSD
+// license.  You should have received a copy of this license along
+// with this source code in a file named "LICENSE."
+//
+////////////////////////////////////////////////////////////////////
+
+// We map X and Y to A and B, because we might change our minds about which
+// is dominant, and we map get/set functions for the channel in question to
+// GETVAL/SETVAL.
+
+
+static void
+FUNCTION_NAME(IMAGETYPE &dest, const IMAGETYPE &source,
+              double width, FilterFunction *make_filter, int channel) {
+  if (!dest.is_valid() || !source.is_valid()) {
+    return;
+  }
+
+  // First, set up a 2-d column-major matrix of StoreTypes, big enough to hold
+  // the image xelvals scaled in the A direction only.  This will hold the
+  // adjusted xel data from our first pass.
+
+  typedef StoreType *StoreTypeP;
+  StoreType **matrix = (StoreType **)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType *));
+  StoreType **matrix_weight = (StoreType **)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType *));
+
+  int a, b;
+
+  for (a=0; a<dest.ASIZE(); a++) {
+    matrix[a] = (StoreType *)PANDA_MALLOC_ARRAY(source.BSIZE() * sizeof(StoreType));
+    matrix_weight[a] = (StoreType *)PANDA_MALLOC_ARRAY(source.BSIZE() * sizeof(StoreType));
+  }
+
+  // First, scale the image in the A direction.
+  double scale;
+  StoreType *temp_source, *temp_source_weight, *temp_dest, *temp_dest_weight;
+
+  scale = (double)dest.ASIZE() / (double)source.ASIZE();
+  temp_source = (StoreType *)PANDA_MALLOC_ARRAY(source.ASIZE() * sizeof(StoreType));
+  temp_source_weight = (StoreType *)PANDA_MALLOC_ARRAY(source.ASIZE() * sizeof(StoreType));
+  temp_dest = (StoreType *)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType));
+  temp_dest_weight = (StoreType *)PANDA_MALLOC_ARRAY(dest.ASIZE() * sizeof(StoreType));
+
+  WorkType *filter;
+  double filter_width;
+
+  make_filter(scale, width, filter, filter_width);
+  memset(temp_source_weight, 0, source.ASIZE() * sizeof(StoreType));
+
+  for (b = 0; b < source.BSIZE(); b++) {
+    for (a = 0; a < source.ASIZE(); a++) {
+      if (source.HASVAL(a, b)) {
+        temp_source[a] = (StoreType)(source_max * source.GETVAL(a, b, channel));
+        temp_source_weight[a] = filter_max;
+      }
+    }
+
+    filter_sparse_row(temp_dest, temp_dest_weight, dest.ASIZE(),
+                      temp_source, temp_source_weight, source.ASIZE(),
+                      scale,
+                      filter, filter_width);
+    
+    for (a = 0; a < dest.ASIZE(); a++) {
+      matrix[a][b] = temp_dest[a];
+      matrix_weight[a][b] = temp_dest_weight[a];
+    }
+  }
+
+  PANDA_FREE_ARRAY(temp_source); 
+  PANDA_FREE_ARRAY(temp_source_weight);
+  PANDA_FREE_ARRAY(temp_dest);
+  PANDA_FREE_ARRAY(temp_dest_weight);
+  PANDA_FREE_ARRAY(filter);
+
+  // Now, scale the image in the B direction.
+  scale = (double)dest.BSIZE() / (double)source.BSIZE();
+  temp_dest = (StoreType *)PANDA_MALLOC_ARRAY(dest.BSIZE() * sizeof(StoreType));
+  temp_dest_weight = (StoreType *)PANDA_MALLOC_ARRAY(dest.BSIZE() * sizeof(StoreType));
+
+  make_filter(scale, width, filter, filter_width);
+
+  for (a = 0; a < dest.ASIZE(); a++) {
+    filter_sparse_row(temp_dest, temp_dest_weight, dest.BSIZE(),
+                      matrix[a], matrix_weight[a], source.BSIZE(),
+                      scale,
+                      filter, filter_width);
+
+    for (b = 0; b < dest.BSIZE(); b++) {
+      dest.SETVAL(a, b, channel, (double)temp_dest[b]/(double)source_max);
+    }
+  }
+
+  PANDA_FREE_ARRAY(temp_dest);
+  PANDA_FREE_ARRAY(temp_dest_weight);
+  PANDA_FREE_ARRAY(filter);
+
+  // Now, clean up our temp matrix and go home!
+
+  for (a = 0; a < dest.ASIZE(); a++) {
+    PANDA_FREE_ARRAY(matrix[a]);
+    PANDA_FREE_ARRAY(matrix_weight[a]);
+  }
+  PANDA_FREE_ARRAY(matrix);
+  PANDA_FREE_ARRAY(matrix_weight);
+}
+

panda/src/pnmimage/pnm-image-filter.cxx

@@ -37,6 +37,7 @@
 #include "cmath.h"
 
 #include "pnmImage.h"
+#include "pfmFile.h"
 
 // WorkType is an abstraction that allows the filtering process to be
 // recompiled to use either floating-point or integer arithmetic.  On SGI
@@ -121,8 +122,8 @@ filter_row(StoreType dest[], int dest_len,
   // have a non-zero offset.
   int offset = (int)cfloor(iscale*0.5);
 
-  for (int dest_x=0; dest_x<dest_len; dest_x++) {
-    double center = (dest_x-offset)/scale;
+  for (int dest_x = 0; dest_x < dest_len; dest_x++) {
+    double center = (dest_x - offset) / scale;
 
     // left and right are the starting and ending ranges of the radius of
     // interest of the filter function.  We need to apply the filter to each
@@ -142,19 +143,19 @@ filter_row(StoreType dest[], int dest_len,
     // This loop is broken into two pieces--the left of center and the right
     // of center--so we don't have to incur the overhead of calling fabs()
     // each time through the loop.
-    for (source_x=left; source_x<right_center; source_x++) {
-      index = (int)(iscale*(center-source_x));
+    for (source_x = left; source_x < right_center; source_x++) {
+      index = (int)(iscale * (center - source_x));
       net_value += filter[index] * source[source_x];
       net_weight += filter[index];
     }
 
-    for (; source_x<=right; source_x++) {
-      index = (int)(iscale*(source_x-center));
+    for (; source_x <= right; source_x++) {
+      index = (int)(iscale * (source_x - center));
       net_value += filter[index] * source[source_x];
       net_weight += filter[index];
     }
 
-    if (net_weight>0) {
+    if (net_weight > 0) {
       dest[dest_x] = (StoreType)(net_value / net_weight);
     } else {
       dest[dest_x] = 0;
@@ -163,6 +164,68 @@ filter_row(StoreType dest[], int dest_len,
   Thread::consider_yield();
 }
 
+// As above, but we also accept an array of weight values per
+// element, to support scaling a sparse array (as in a PfmFile).
+static void
+filter_sparse_row(StoreType dest[], StoreType dest_weight[], int dest_len,
+                  const StoreType source[], const StoreType source_weight[], int source_len,
+                  double scale,                    //  == dest_len / source_len
+                  const WorkType filter[],
+                  double filter_width) {
+  // If we are expanding the row (scale>1.0), we need to look at a fractional
+  // granularity.  Hence, we scale our filter index by scale.  If we are
+  // compressing (scale<1.0), we don't need to fiddle with the filter index, so
+  // we leave it at one.
+  double iscale = max(scale, 1.0);
+
+  // Similarly, if we are expanding the row, we want to start the new row at
+  // the far left edge of the original pixel, not in the center.  So we will
+  // have a non-zero offset.
+  int offset = (int)cfloor(iscale*0.5);
+
+  for (int dest_x = 0; dest_x < dest_len; dest_x++) {
+    double center = (dest_x - offset) / scale;
+
+    // left and right are the starting and ending ranges of the radius of
+    // interest of the filter function.  We need to apply the filter to each
+    // value in this range.
+    int left = max((int)cfloor(center - filter_width), 0);
+    int right = min((int)cceil(center + filter_width), source_len-1);
+
+    // right_center is the point just to the right of the center.  This
+    // allows us to flip the sign of the offset when we cross the center point.
+    int right_center = (int)cceil(center);
+
+    WorkType net_weight = 0;
+    WorkType net_value = 0;
+
+    int index, source_x;
+
+    // This loop is broken into two pieces--the left of center and the right
+    // of center--so we don't have to incur the overhead of calling fabs()
+    // each time through the loop.
+    for (source_x = left; source_x < right_center; source_x++) {
+      index = (int)(iscale * (center - source_x));
+      net_value += filter[index] * source[source_x] * source_weight[source_x];
+      net_weight += filter[index] * source_weight[source_x];
+    }
+
+    for (; source_x <= right; source_x++) {
+      index = (int)(iscale * (source_x - center));
+      net_value += filter[index] * source[source_x] * source_weight[source_x];
+      net_weight += filter[index] * source_weight[source_x];
+    }
+
+    if (net_weight > 0) {
+      dest[dest_x] = (StoreType)(net_value / net_weight);
+    } else {
+      dest[dest_x] = 0;
+    }
+    dest_weight[dest_x] = (StoreType)net_weight;
+  }
+  Thread::consider_yield();
+}
+
 
 // The various filter functions are called before each axis scaling to build
 // an kernel array suitable for the given scaling factor.  Given a scaling
@@ -186,20 +249,20 @@ box_filter_impl(double scale, double width,
     // the filter function to prevent dropping below the Nyquist rate.
     // Hence, we divide by scale.
     fscale = 1.0 / scale;
-  } else {
 
+  } else {
     // If we are expanding the image, we want to increase the granularity
     // of the filter function since we will need to access fractional cel
     // values.  Hence, we multiply by scale.
     fscale = scale;
   }
   filter_width = width;
-  int actual_width = (int)cceil((filter_width+1) * fscale);
+  int actual_width = (int)cceil((filter_width + 1) * fscale);
 
   filter = (WorkType *)PANDA_MALLOC_ARRAY(actual_width * sizeof(WorkType));
 
-  for (int i=0; i<actual_width; i++) {
-    filter[i] = (i<=filter_width*fscale) ? filter_max : 0;
+  for (int i = 0; i < actual_width; i++) {
+    filter[i] = (i <= filter_width * fscale) ? filter_max : 0;
   }
 }
 
@@ -219,9 +282,10 @@ gaussian_filter_impl(double scale, double width,
     // values.  Hence, we multiply by scale (to make fscale larger).
     fscale = scale;
   }
+
   double sigma = width/2;
   filter_width = 3.0 * sigma;
-  int actual_width = (int)cceil((filter_width+1) * fscale);
+  int actual_width = (int)cceil((filter_width + 1) * fscale);
 
   // G(x, y) = (1/(2 pi sigma^2)) * exp( - (x^2 + y^2) / (2 sigma^2))
 
@@ -230,10 +294,10 @@ gaussian_filter_impl(double scale, double width,
   // so we can ignore the y^2.)
 
   filter = (WorkType *)PANDA_MALLOC_ARRAY(actual_width * sizeof(WorkType));
-  double div = 2*sigma*sigma;
+  double div = 2 * sigma * sigma;
 
-  for (int i=0; i<actual_width; i++) {
-    double x = i/fscale;
+  for (int i = 0; i < actual_width; i++) {
+    double x = i / fscale;
     filter[i] = (WorkType)(filter_max * exp(-x*x / div));
     // The highest value of the exp function in this range is always 1.0,
     // at index value 0.  Thus, we scale the whole range by filter_max,
@@ -267,126 +331,146 @@ gaussian_filter_impl(double scale, double width,
 // These instances scale by X first, then by Y.
 
 #define FUNCTION_NAME filter_red_xy
+#define IMAGETYPE PNMImage
 #define ASIZE get_x_size
 #define BSIZE get_y_size
-#define GETVAL(a, b) get_red(a, b)
-#define SETVAL(a, b, v) set_red(a, b, v)
+#define GETVAL(a, b, channel) get_red(a, b)
+#define SETVAL(a, b, channel, v) set_red(a, b, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_green_xy
+#define IMAGETYPE PNMImage
 #define ASIZE get_x_size
 #define BSIZE get_y_size
-#define GETVAL(a, b) get_green(a, b)
-#define SETVAL(a, b, v) set_green(a, b, v)
+#define GETVAL(a, b, channel) get_green(a, b)
+#define SETVAL(a, b, channel, v) set_green(a, b, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_blue_xy
+#define IMAGETYPE PNMImage
 #define ASIZE get_x_size
 #define BSIZE get_y_size
-#define GETVAL(a, b) get_blue(a, b)
-#define SETVAL(a, b, v) set_blue(a, b, v)
+#define GETVAL(a, b, channel) get_blue(a, b)
+#define SETVAL(a, b, channel, v) set_blue(a, b, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_gray_xy
+#define IMAGETYPE PNMImage
 #define ASIZE get_x_size
 #define BSIZE get_y_size
-#define GETVAL(a, b) get_bright(a, b)
-#define SETVAL(a, b, v) set_xel(a, b, v)
+#define GETVAL(a, b, channel) get_bright(a, b)
+#define SETVAL(a, b, channel, v) set_xel(a, b, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_alpha_xy
+#define IMAGETYPE PNMImage
 #define ASIZE get_x_size
 #define BSIZE get_y_size
-#define GETVAL(a, b) get_alpha(a, b)
-#define SETVAL(a, b, v) set_alpha(a, b, v)
+#define GETVAL(a, b, channel) get_alpha(a, b)
+#define SETVAL(a, b, channel, v) set_alpha(a, b, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 
 // These instances scale by Y first, then by X.
 
 #define FUNCTION_NAME filter_red_yx
+#define IMAGETYPE PNMImage
 #define ASIZE get_y_size
 #define BSIZE get_x_size
-#define GETVAL(a, b) get_red(b, a)
-#define SETVAL(a, b, v) set_red(b, a, v)
+#define GETVAL(a, b, channel) get_red(b, a)
+#define SETVAL(a, b, channel, v) set_red(b, a, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_green_yx
+#define IMAGETYPE PNMImage
 #define ASIZE get_y_size
 #define BSIZE get_x_size
-#define GETVAL(a, b) get_green(b, a)
-#define SETVAL(a, b, v) set_green(b, a, v)
+#define GETVAL(a, b, channel) get_green(b, a)
+#define SETVAL(a, b, channel, v) set_green(b, a, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_blue_yx
+#define IMAGETYPE PNMImage
 #define ASIZE get_y_size
 #define BSIZE get_x_size
-#define GETVAL(a, b) get_blue(b, a)
-#define SETVAL(a, b, v) set_blue(b, a, v)
+#define GETVAL(a, b, channel) get_blue(b, a)
+#define SETVAL(a, b, channel, v) set_blue(b, a, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_gray_yx
+#define IMAGETYPE PNMImage
 #define ASIZE get_y_size
 #define BSIZE get_x_size
-#define GETVAL(a, b) get_bright(b, a)
-#define SETVAL(a, b, v) set_xel(b, a, v)
+#define GETVAL(a, b, channel) get_bright(b, a)
+#define SETVAL(a, b, channel, v) set_xel(b, a, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 #define FUNCTION_NAME filter_alpha_yx
+#define IMAGETYPE PNMImage
 #define ASIZE get_y_size
 #define BSIZE get_x_size
-#define GETVAL(a, b) get_alpha(b, a)
-#define SETVAL(a, b, v) set_alpha(b, a, v)
+#define GETVAL(a, b, channel) get_alpha(b, a)
+#define SETVAL(a, b, channel, v) set_alpha(b, a, v)
 #include "pnm-image-filter-core.cxx"
 #undef SETVAL
 #undef GETVAL
 #undef BSIZE
 #undef ASIZE
+#undef IMAGETYPE
 #undef FUNCTION_NAME
 
 
@@ -399,30 +483,35 @@ filter_image(PNMImage &dest, const PNMImage &source,
   // We want to scale by the smallest destination axis first, for a
   // slight performance gain.
 
+  // In the PNMImage case (unlike the PfmFile case), the channel
+  // parameter is not used.  We *could* use it to avoid the
+  // replication of quite so many functions, but we replicate them
+  // anyway, for another tiny performance gain.
+
   if (dest.get_x_size() <= dest.get_y_size()) {
     if (dest.is_grayscale() || source.is_grayscale()) {
-      filter_gray_xy(dest, source, width, make_filter);
+      filter_gray_xy(dest, source, width, make_filter, 0);
     } else {
-      filter_red_xy(dest, source, width, make_filter);
-      filter_green_xy(dest, source, width, make_filter);
-      filter_blue_xy(dest, source, width, make_filter);
+      filter_red_xy(dest, source, width, make_filter, 0);
+      filter_green_xy(dest, source, width, make_filter, 0);
+      filter_blue_xy(dest, source, width, make_filter, 0);
     }
 
     if (dest.has_alpha() && source.has_alpha()) {
-      filter_alpha_xy(dest, source, width, make_filter);
+      filter_alpha_xy(dest, source, width, make_filter, 0);
     }
 
   } else {
     if (dest.is_grayscale() || source.is_grayscale()) {
-      filter_gray_yx(dest, source, width, make_filter);
+      filter_gray_yx(dest, source, width, make_filter, 0);
     } else {
-      filter_red_yx(dest, source, width, make_filter);
-      filter_green_yx(dest, source, width, make_filter);
-      filter_blue_yx(dest, source, width, make_filter);
+      filter_red_yx(dest, source, width, make_filter, 0);
+      filter_green_yx(dest, source, width, make_filter, 0);
+      filter_blue_yx(dest, source, width, make_filter, 0);
     }
 
     if (dest.has_alpha() && source.has_alpha()) {
-      filter_alpha_yx(dest, source, width, make_filter);
+      filter_alpha_yx(dest, source, width, make_filter, 0);
     }
   }
 }
@@ -457,6 +546,137 @@ gaussian_filter_from(double width, const PNMImage &copy) {
   filter_image(*this, copy, width, &gaussian_filter_impl);
 }
 
+// Now we do it again, this time for PfmFile.  In this case we also
+// need to support the sparse variants, since PfmFiles can be
+// incomplete.  However, we don't need to have a different function
+// for each channel.
+
+#define FUNCTION_NAME filter_pfm_xy
+#define IMAGETYPE PfmFile
+#define ASIZE get_x_size
+#define BSIZE get_y_size
+#define GETVAL(a, b, channel) get_component(a, b, channel)
+#define SETVAL(a, b, channel, v) set_component(a, b, channel, v)
+#include "pnm-image-filter-core.cxx"
+#undef SETVAL
+#undef GETVAL
+#undef BSIZE
+#undef ASIZE
+#undef IMAGETYPE
+#undef FUNCTION_NAME
+
+#define FUNCTION_NAME filter_pfm_yx
+#define IMAGETYPE PfmFile
+#define ASIZE get_y_size
+#define BSIZE get_x_size
+#define GETVAL(a, b, channel) get_component(b, a, channel)
+#define SETVAL(a, b, channel, v) set_component(b, a, channel, v)
+#include "pnm-image-filter-core.cxx"
+#undef SETVAL
+#undef GETVAL
+#undef BSIZE
+#undef ASIZE
+#undef IMAGETYPE
+#undef FUNCTION_NAME
+
+
+#define FUNCTION_NAME filter_pfm_sparse_xy
+#define IMAGETYPE PfmFile
+#define ASIZE get_x_size
+#define BSIZE get_y_size
+#define HASVAL(a, b) has_point(a, b)
+#define GETVAL(a, b, channel) get_component(a, b, channel)
+#define SETVAL(a, b, channel, v) set_component(a, b, channel, v)
+#include "pnm-image-filter-sparse-core.cxx"
+#undef SETVAL
+#undef GETVAL
+#undef HASVAL
+#undef BSIZE
+#undef ASIZE
+#undef IMAGETYPE
+#undef FUNCTION_NAME
+
+#define FUNCTION_NAME filter_pfm_sparse_yx
+#define IMAGETYPE PfmFile
+#define ASIZE get_y_size
+#define BSIZE get_x_size
+#define HASVAL(a, b) has_point(b, a)
+#define GETVAL(a, b, channel) get_component(b, a, channel)
+#define SETVAL(a, b, channel, v) set_component(b, a, channel, v)
+#include "pnm-image-filter-sparse-core.cxx"
+#undef SETVAL
+#undef GETVAL
+#undef HASVAL
+#undef BSIZE
+#undef ASIZE
+#undef IMAGETYPE
+#undef FUNCTION_NAME
+
+
+// filter_image pulls everything together, and filters one image into
+// another.  Both images can be the same with no ill effects.
+static void
+filter_image(PfmFile &dest, const PfmFile &source,
+             double width, FilterFunction *make_filter) {
+  int num_channels = min(dest.get_num_channels(), source.get_num_channels());
+
+  if (source.has_no_data_value()) {
+    // We need to use the sparse variant.
+    if (dest.get_x_size() <= dest.get_y_size()) {
+      for (int ci = 0; ci < num_channels; ++ci) {
+        filter_pfm_sparse_xy(dest, source, width, make_filter, ci);
+      }
+      
+    } else {
+      for (int ci = 0; ci < num_channels; ++ci) {
+        filter_pfm_sparse_yx(dest, source, width, make_filter, ci);
+      }
+    }
+  } else {
+    // We can use the slightly faster fully-specified variant.
+    if (dest.get_x_size() <= dest.get_y_size()) {
+      for (int ci = 0; ci < num_channels; ++ci) {
+        filter_pfm_xy(dest, source, width, make_filter, ci);
+      }
+      
+    } else {
+      for (int ci = 0; ci < num_channels; ++ci) {
+        filter_pfm_yx(dest, source, width, make_filter, ci);
+      }
+    }
+  }
+}
+
+
+
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::box_filter_from
+//       Access: Public
+//  Description: Makes a resized copy of the indicated image into this
+//               one using the indicated filter.  The image to be
+//               copied is squashed and stretched to match the
+//               dimensions of the current image, applying the
+//               appropriate filter to perform the stretching.
+////////////////////////////////////////////////////////////////////
+void PfmFile::
+box_filter_from(double width, const PfmFile &copy) {
+  filter_image(*this, copy, width, &box_filter_impl);
+}
+
+////////////////////////////////////////////////////////////////////
+//     Function: PfmFile::gaussian_filter_from
+//       Access: Public
+//  Description: Makes a resized copy of the indicated image into this
+//               one using the indicated filter.  The image to be
+//               copied is squashed and stretched to match the
+//               dimensions of the current image, applying the
+//               appropriate filter to perform the stretching.
+////////////////////////////////////////////////////////////////////
+void PfmFile::
+gaussian_filter_from(double width, const PfmFile &copy) {
+  filter_image(*this, copy, width, &gaussian_filter_impl);
+}
+
 
 //
 // The following functions are support for quick_box_filter().