ImageBufAlgo: Initial implementation of st_warp transform function

src/include/OpenImageIO/imagebufalgo.h

@@ -790,6 +790,68 @@ bool OIIO_API warp (ImageBuf &dst, const ImageBuf &src, M33fParam M,
 /// @}
 
 
+/// @defgroup st_warp (st_warp: warp an image using per-pixel st coordinates)
+/// @{
+///
+/// Warp the `src` image using "st" coordinates from a secondary `stbuf` image.
+///
+/// Each pixel in the `stbuf` image is used as a normalized image-space
+/// coordinate in the `src` image, which is then sampled at that position using
+/// the given reconstruction filter to produce an output pixel.
+///
+/// The transform is only defined over the area of the `stbuf` image, and thus
+/// the given `roi` argument will be intersected with its geometry.
+///
+/// \b NOTE: The current behavior of this transform is modeled to match Nuke's
+/// STMap node.
+///
+/// @param dst
+///             The output ImageBuf. If an initialized buffer is provided, its
+///             full-size dimensions must match those of `stbuf`.
+/// @param src
+///             The source ImageBuf to warp.
+/// @param stbuf
+///             The ImageBuf holding the st coordinates. This must be holding
+///             a floating-point pixel data type.
+/// @param chan_s
+///             The index of the "s" channel in the `stbuf` image. This defaults
+///             to its first channel.
+/// @param chan_t
+///             The index of the "t" channel in the `stbuf` image. This defaults
+///             to its second channel.
+/// @param flip_s
+///             Whether to mirror the "s" coordinate along the horizontal axis
+///             when computing source pixel positions. This is useful if the
+///             coordinates are defined in terms of a different image origin 
+///             than OpenImageIO's.
+/// @param flip_t
+///             Whether to mirror the "t" coordinate along the vertical axis
+///             when computing source pixel positions. This is useful if the
+///             coordinates are defined in terms of a different image origin 
+///             than OpenImageIO's.
+
+ImageBuf OIIO_API st_warp (const ImageBuf &src, const ImageBuf& stbuf,
+                           string_view filtername=string_view(),
+                           float filterwidth=0.0f, int chan_s=0, int chan_t=1,
+                           bool flip_s=false, bool flip_t=false, ROI roi={},
+                           int nthreads=0);
+ImageBuf OIIO_API st_warp (const ImageBuf &src, const ImageBuf& stbuf,
+                           Filter2D *filter, int chan_s=0, int chan_t=1,
+                           bool flip_s=false, bool flip_t=false, ROI roi={},
+                           int nthreads=0);
+bool OIIO_API st_warp (ImageBuf &dst, const ImageBuf &src,
+                       const ImageBuf& stbuf,
+                       string_view filtername=string_view(),
+                       float filterwidth=0.0f, int chan_s=0, int chan_t=1,
+                       bool flip_s=false, bool flip_t=false, ROI roi={},
+                       int nthreads=0);
+bool OIIO_API st_warp (ImageBuf &dst, const ImageBuf &src,
+                       const ImageBuf& stbuf, Filter2D *filter, int chan_s=0,
+                       int chan_t=1, bool flip_s=false, bool flip_t=false,
+                       ROI roi={}, int nthreads=0);
+/// @}
+
+
 /// Compute per-pixel sum `A + B`, returning the result image.
 ///
 /// `A` and `B` may each either be an `ImageBuf&`, or a `cspan<float>`

src/libOpenImageIO/imagebufalgo_xform.cpp

@@ -205,6 +205,30 @@ filtered_sample(const ImageBuf& src, float s, float t, float dsdx, float dtdx,
 
 
 
+static std::shared_ptr<Filter2D>
+get_warp_filter(string_view filtername_, float filterwidth, ImageBuf& dst)
+{
+    // Set up a shared pointer with custom deleter to make sure any
+    // filter we allocate here is properly destroyed.
+    std::shared_ptr<Filter2D> filter((Filter2D*)nullptr, Filter2D::destroy);
+    std::string filtername = filtername_.size() ? filtername_ : "lanczos3";
+    for (int i = 0, e = Filter2D::num_filters(); i < e; ++i) {
+        FilterDesc fd;
+        Filter2D::get_filterdesc(i, &fd);
+        if (fd.name == filtername) {
+            float w = filterwidth > 0.0f ? filterwidth : fd.width;
+            filter.reset(Filter2D::create(filtername, w, w));
+            break;
+        }
+    }
+    if (!filter) {
+        dst.errorfmt("Filter \"{}\" not recognized", filtername);
+    }
+    return filter;
+}
+
+
+
 template<typename DSTTYPE, typename SRCTYPE>
 static bool
 warp_(ImageBuf& dst, const ImageBuf& src, const Imath::M33f& M,
@@ -287,29 +311,16 @@ ImageBufAlgo::warp(ImageBuf& dst, const ImageBuf& src, M33fParam M,
 
 bool
 ImageBufAlgo::warp(ImageBuf& dst, const ImageBuf& src, M33fParam M,
-                   string_view filtername_, float filterwidth,
+                   string_view filtername, float filterwidth,
                    bool recompute_roi, ImageBuf::WrapMode wrap, ROI roi,
                    int nthreads)
 {
     // Set up a shared pointer with custom deleter to make sure any
     // filter we allocate here is properly destroyed.
-    std::shared_ptr<Filter2D> filter((Filter2D*)NULL, Filter2D::destroy);
-    std::string filtername = filtername_.size() ? filtername_ : "lanczos3";
-    for (int i = 0, e = Filter2D::num_filters(); i < e; ++i) {
-        FilterDesc fd;
-        Filter2D::get_filterdesc(i, &fd);
-        if (fd.name == filtername) {
-            float w = filterwidth > 0.0f ? filterwidth : fd.width;
-            float h = filterwidth > 0.0f ? filterwidth : fd.width;
-            filter.reset(Filter2D::create(filtername, w, h));
-            break;
-        }
-    }
+    auto filter = get_warp_filter(filtername, filterwidth, dst);
     if (!filter) {
-        dst.errorfmt("Filter \"{}\" not recognized", filtername);
-        return false;
+        return false;  // error issued in get_warp_filter
     }
-
     return warp(dst, src, M, filter.get(), recompute_roi, wrap, roi, nthreads);
 }
 
@@ -1144,4 +1155,255 @@ ImageBufAlgo::rotate(const ImageBuf& src, float angle, string_view filtername,
 }
 
 
+
+template<typename DSTTYPE, typename SRCTYPE, typename STTYPE>
+static bool
+st_warp_(ImageBuf& dst, const ImageBuf& src, const ImageBuf& stbuf, int chan_s,
+         int chan_t, bool flip_s, bool flip_t, Filter2D* filter, ROI roi,
+         int nthreads)
+{
+    OIIO_DASSERT(filter);
+    OIIO_DASSERT(dst.spec().nchannels >= roi.chend);
+
+    ImageBufAlgo::parallel_image(roi, nthreads, [&](ROI roi) {
+        const ImageSpec& srcspec(src.spec());
+        const ImageSpec& dstspec(dst.spec());
+        const int src_width  = srcspec.full_width;
+        const int src_height = srcspec.full_height;
+
+        const float xscale = float(dstspec.full_width) / src_width;
+        const float yscale = float(dstspec.full_height) / src_height;
+
+        const int xbegin = src.xbegin();
+        const int xend   = src.xend();
+        const int ybegin = src.ybegin();
+        const int yend   = src.yend();
+
+        // The horizontal and vertical filter radii, in source pixels.
+        // We will sample and filter the source over
+        //   [x-filterrad_x, x+filterrad_x] X [y-filterrad_y,y+filterrad_y].
+        const int filterrad_x = (int)ceilf(filter->width() / 2.0f / xscale);
+        const int filterrad_y = (int)ceilf(filter->height() / 2.0f / yscale);
+
+        ImageBuf::ConstIterator<SRCTYPE> src_iter(src);
+        ImageBuf::ConstIterator<STTYPE> st_iter(stbuf, roi);
+        ImageBuf::Iterator<DSTTYPE> out_iter(dst, roi);
+
+        // Accumulation buffer for filter samples, typed to maintain the
+        // necessary precision.
+        typedef typename Accum_t<DSTTYPE>::type Acc_t;
+        const int nchannels = roi.chend - roi.chbegin;
+        Acc_t* sample_accum = OIIO_ALLOCA(Acc_t, nchannels);
+
+        // The ST buffer defines the output dimensions, and thus the bounds of
+        // the outer loop.
+        // XXX: Sampling of the source buffer can be entirely random, so there
+        // are probably some opportunities for optimization in here...
+        for (; !st_iter.done(); ++st_iter) {
+            // Look up source coordinates from ST channels.
+            // We don't care about faithfully maintaining `STTYPE`: Half will be
+            // promoted accurately, and double isn't really useful, since filter
+            // lookups don't support that (excessive) level of precision.
+            float src_s = st_iter[chan_s];
+            float src_t = st_iter[chan_t];
+
+            if (flip_s) {
+                src_s = 1.0f - src_s;
+            }
+            if (flip_t) {
+                src_t = 1.0f - src_t;
+            }
+
+            const float src_x = src_s * src_width;
+            const float src_y = src_t * src_height;
+
+            // Set up source iterator range
+            const int x_min = clamp((int)floorf(src_x - filterrad_x), xbegin,
+                                    xend);
+            const int x_max = clamp((int)ceilf(src_x + filterrad_x), xbegin,
+                                    xend);
+            const int y_min = clamp((int)floorf(src_y - filterrad_y), ybegin,
+                                    yend);
+            const int y_max = clamp((int)ceilf(src_y + filterrad_y), ybegin,
+                                    yend);
+
+            src_iter.rerange(x_min, x_max + 1, y_min, y_max + 1, 0, 1);
+
+            memset(sample_accum, 0, nchannels * sizeof(Acc_t));
+            float total_weight = 0.0f;
+            for (; !src_iter.done(); ++src_iter) {
+                const float weight = (*filter)(src_iter.x() - src_x + 0.5f,
+                                               src_iter.y() - src_y + 0.5f);
+                total_weight += weight;
+                for (int idx = 0, chan = roi.chbegin; chan < roi.chend;
+                     ++chan, ++idx) {
+                    sample_accum[idx] += src_iter[chan] * weight;
+                }
+            }
+
+            if (total_weight > 0.0f) {
+                for (int idx = 0, chan = roi.chbegin; chan < roi.chend;
+                     ++chan, ++idx) {
+                    out_iter[chan] = sample_accum[idx] / total_weight;
+                }
+            } else {
+                for (int chan = roi.chbegin; chan < roi.chend; ++chan) {
+                    out_iter[chan] = 0;
+                }
+            }
+            ++out_iter;
+        }
+    });  // end of parallel_image
+    return true;
+}
+
+
+
+static bool
+check_st_warp_args(ImageBuf& dst, const ImageBuf& src, const ImageBuf& stbuf,
+                   int chan_s, int chan_t, ROI& roi)
+{
+    // Validate ST buffer
+    if (!stbuf.initialized()) {
+        dst.error("ImageBufAlgo::st_warp : Uninitialized ST buffer");
+        return false;
+    }
+
+    const ImageSpec& stSpec(stbuf.spec());
+    // XXX: Wanted to use `uint32_t` for channel indices, but I don't want to
+    // break from the rest of the API and introduce a bunch of compile warnings.
+    if (chan_s >= stSpec.nchannels) {
+        dst.errorfmt("ImageBufAlgo::st_warp : Out-of-range S channel index: {}",
+                     chan_s);
+        return false;
+    }
+    if (chan_t >= stSpec.nchannels) {
+        dst.errorfmt("ImageBufAlgo::st_warp : Out-of-range T channel index: {}",
+                     chan_t);
+        return false;
+    }
+    // N.B. We currently require floating-point ST channels
+    if (!stSpec.format.is_floating_point()) {
+        dst.error("ImageBufAlgo::st_warp : ST buffer must be holding floating-"
+                  "point data");
+        return false;
+    }
+
+    if (dst.initialized()) {
+        // XXX: Currently requiring the pixel scales of ST and a preallocated
+        // output buffer to match.
+        if (dst.spec().full_width != stSpec.full_width
+            || dst.spec().full_height != stSpec.full_height) {
+            dst.error("ImageBufAlgo::st_warp : Output and ST buffers must have "
+                      "the same full width and height");
+            return false;
+        }
+    }
+
+    // Prep the dest spec using the channels from `src`, and the intersection of
+    // `roi` and the ROI from the ST buffer (since the ST warp is only defined
+    // for pixels in the ST buffer). We grab a copy of the input ROI before
+    // `IBAprep`, since we want to intersect it ourselves.
+    ROI inputROI(roi);
+    bool res
+        = ImageBufAlgo::IBAprep(roi, &dst, &src,
+                                ImageBufAlgo::IBAprep_NO_SUPPORT_VOLUME
+                                    | ImageBufAlgo::IBAprep_NO_COPY_ROI_FULL);
+    if (res) {
+        roi         = roi_intersection(inputROI, stSpec.roi());
+        roi.chbegin = std::max(roi.chbegin, 0);
+        roi.chend   = std::min(roi.chend, src.spec().nchannels);
+
+        ImageSpec destSpec(dst.spec());
+        destSpec.set_roi(roi);
+        destSpec.set_roi_full(stSpec.roi_full());
+        // XXX: It would be nice to be able to tell `IBAprep` to skip the buffer
+        // initialization. Worth adding a new prep flag?
+        dst.reset(destSpec);
+    }
+    return res;
+}
+
+
+
+bool
+ImageBufAlgo::st_warp(ImageBuf& dst, const ImageBuf& src, const ImageBuf& stbuf,
+                      Filter2D* filter, int chan_s, int chan_t, bool flip_s,
+                      bool flip_t, ROI roi, int nthreads)
+{
+    pvt::LoggedTimer logtime("IBA::st_warp");
+
+    if (!check_st_warp_args(dst, src, stbuf, chan_s, chan_t, roi)) {
+        return false;
+    }
+
+    // Set up a shared pointer with custom deleter to make sure any
+    // filter we allocate here is properly destroyed.
+    std::shared_ptr<Filter2D> filterptr((Filter2D*)nullptr, Filter2D::destroy);
+    if (!filter) {
+        // If a null filter was provided, fall back to a reasonable default.
+        filterptr.reset(Filter2D::create("lanczos3", 6.0f, 6.0f));
+        filter = filterptr.get();
+    }
+
+    bool ok;
+    OIIO_DISPATCH_COMMON_TYPES3(ok, "st_warp", st_warp_, dst.spec().format,
+                                src.spec().format, stbuf.spec().format, dst,
+                                src, stbuf, chan_s, chan_t, flip_s, flip_t,
+                                filter, roi, nthreads);
+    return ok;
+}
+
+
+
+bool
+ImageBufAlgo::st_warp(ImageBuf& dst, const ImageBuf& src, const ImageBuf& stbuf,
+                      string_view filtername, float filterwidth, int chan_s,
+                      int chan_t, bool flip_s, bool flip_t, ROI roi,
+                      int nthreads)
+{
+    // Set up a shared pointer with custom deleter to make sure any
+    // filter we allocate here is properly destroyed.
+    auto filter = get_warp_filter(filtername, filterwidth, dst);
+    if (!filter) {
+        return false;  // Error issued in `get_warp_filter`.
+    }
+    return st_warp(dst, src, stbuf, filter.get(), chan_s, chan_t, flip_s,
+                   flip_t, roi, nthreads);
+}
+
+
+
+ImageBuf
+ImageBufAlgo::st_warp(const ImageBuf& src, const ImageBuf& stbuf,
+                      Filter2D* filter, int chan_s, int chan_t, bool flip_s,
+                      bool flip_t, ROI roi, int nthreads)
+{
+    ImageBuf result;
+    bool ok = st_warp(result, src, stbuf, filter, chan_s, chan_t, flip_s,
+                      flip_t, roi, nthreads);
+    if (!ok && !result.has_error()) {
+        result.error("ImageBufAlgo::st_warp : Unknown error");
+    }
+    return result;
+}
+
+
+
+ImageBuf
+ImageBufAlgo::st_warp(const ImageBuf& src, const ImageBuf& stbuf,
+                      string_view filtername, float filterwidth, int chan_s,
+                      int chan_t, bool flip_s, bool flip_t, ROI roi,
+                      int nthreads)
+{
+    ImageBuf result;
+    bool ok = st_warp(result, src, stbuf, filtername, filterwidth, chan_s,
+                      chan_t, flip_s, flip_t, roi, nthreads);
+    if (!ok && !result.has_error()) {
+        result.error("ImageBufAlgo::st_warp : Unknown error");
+    }
+    return result;
+}
+
+
 OIIO_NAMESPACE_END