diff --git a/src/libOpenImageIO/imagebufalgo_test.cpp b/src/libOpenImageIO/imagebufalgo_test.cpp index 42367ae3d2..1bcb3b93fb 100644 --- a/src/libOpenImageIO/imagebufalgo_test.cpp +++ b/src/libOpenImageIO/imagebufalgo_test.cpp @@ -738,7 +738,125 @@ test_zover() -// Test ImageBuf::resample +// Same pixels as `src`, but with a gap between scanlines and, optionally, +// a gap between pixels. `storage` owns the memory and must outlive the result. +// Row padding is legal for a source the tables handle -- contiguous_scanline() +// permits it -- so it has to be read through the real scanline stride. Pixel +// padding is not, and selects the per-pixel fallback instead. +static ImageBuf +padded_copy(const ImageBuf& src, bool pad_pixels, + std::unique_ptr& storage) +{ + const ImageSpec& spec(src.spec()); + stride_t chansize = stride_t(spec.format.size()); + stride_t xstride = chansize * (spec.nchannels + (pad_pixels ? 1 : 0)); + stride_t ystride = xstride * (spec.width + 1); + storage.reset(new char[size_t(ystride) * spec.height]); + memset(storage.get(), 0, size_t(ystride) * spec.height); + ImageBuf padded(spec, storage.get(), xstride, ystride); + ImageBufAlgo::paste(padded, spec.x, spec.y, spec.z, 0, src); + return padded; +} + + + +// The tables are a pure optimization: they must reproduce the per-pixel path +// exactly, so this compares the two rather than checking against a tolerance, +// which would hide the one-pixel shift a mistake in the mapping arithmetic +// actually causes. Pixel padding is what selects the fallback -- there is no +// runtime switch to flip, the dispatch is by capability. +static void +check_against_fallback(const ImageBuf& src, const ImageSpec& dstspec, + bool interpolate) +{ + std::unique_ptr fast_mem, slow_mem; + ImageBuf fast_src = padded_copy(src, false, fast_mem); + ImageBuf slow_src = padded_copy(src, true, slow_mem); + ImageBuf fast(dstspec), slow(dstspec); + OIIO_CHECK_ASSERT(ImageBufAlgo::resample(fast, fast_src, interpolate)); + OIIO_CHECK_ASSERT(ImageBufAlgo::resample(slow, slow_src, interpolate)); + if (memcmp(fast.localpixels(), slow.localpixels(), dstspec.image_bytes()) + == 0) + return; + auto cr = ImageBufAlgo::compare(fast, slow, 0.0f, 0.0f); + Strutil::print(" interp={} {}x{} nch={} type={}: maxerror={:g} " + "({} of {} values differ)\n", + int(interpolate), dstspec.width, dstspec.height, + dstspec.nchannels, dstspec.format, cr.maxerror, cr.nfail, + size_t(dstspec.width) * dstspec.height * dstspec.nchannels); + OIIO_CHECK_ASSERT(false && "resample paths disagree"); +} + + + +// Nearest resampling must pick the same source pixel as the mapping in the +// OIIO docs, and must be black wherever the destination maps outside the +// source data window. The crop case is what distinguishes a data window from +// a full window, and is not exercised by an ordinary resize. +void +test_resample_correctness() +{ + std::cout << "test resample correctness\n"; + + // Source pixel (x,y) holds the value x + 10*y, so a resampled pixel + // identifies exactly which source pixel it came from. + const int srcsize = 4; + ImageBuf src(ImageSpec(srcsize, srcsize, 1, TypeFloat)); + for (ImageBuf::Iterator it(src); !it.done(); ++it) + it[0] = float(it.x() + 10 * it.y()); + + // Both a magnification and a reduction, checked against the mapping + // evaluated independently here. + for (int dstsize : { 8, 2 }) { + ImageBuf dst(ImageSpec(dstsize, dstsize, 1, TypeFloat)); + OIIO_CHECK_ASSERT(ImageBufAlgo::resample(dst, src, false)); + for (ImageBuf::ConstIterator it(dst); !it.done(); ++it) { + int sx = int((it.x() + 0.5f) / dstsize * srcsize); + int sy = int((it.y() + 0.5f) / dstsize * srcsize); + OIIO_CHECK_EQUAL(it[0], float(sx + 10 * sy)); + } + } + + // A source whose data window is a 2x2 crop of a 4x4 full window. The + // destination covers the whole full window, so its outer ring must come + // out black rather than clamped or out of bounds. + ImageSpec cropspec(2, 2, 1, TypeFloat); + cropspec.x = cropspec.y = 1; + cropspec.full_x = cropspec.full_y = 0; + cropspec.full_width = cropspec.full_height = 4; + ImageBuf crop(cropspec); + for (ImageBuf::Iterator it(crop); !it.done(); ++it) + it[0] = float(it.x() + 10 * it.y()); + + ImageBuf dst(ImageSpec(4, 4, 1, TypeFloat)); + OIIO_CHECK_ASSERT(ImageBufAlgo::resample(dst, crop, false)); + for (ImageBuf::ConstIterator it(dst); !it.done(); ++it) { + bool inside = it.x() >= 1 && it.x() <= 2 && it.y() >= 1 && it.y() <= 2; + float expected = inside ? float(it.x() + 10 * it.y()) : 0.0f; + OIIO_CHECK_EQUAL(it[0], expected); + } + + // Against the per-pixel path, for both sampling modes: uint8 because its + // rescaling to [0,1] is the conversion most easily got wrong, float + // because it is the one where a mapping error cannot hide in rounding. + // One reduction and one magnification is enough -- the mapping is the + // same arithmetic in both directions. + for (auto type : { TypeUInt8, TypeFloat }) { + float top[4] = { 0.1f, 0.2f, 0.3f, 0.4f }; + float bottom[4] = { 0.9f, 0.8f, 0.7f, 0.6f }; + ImageBuf grad(ImageSpec(97, 61, 4, type)); + ImageBufAlgo::fill(grad, cspan(top), cspan(bottom)); + for (auto wh : + { std::pair(53, 79), std::pair(211, 43) }) + for (bool interp : { false, true }) + check_against_fallback(grad, + ImageSpec(wh.first, wh.second, 4, type), + interp); + } +} + + + void test_resample() { @@ -748,7 +866,7 @@ test_resample() Benchmarker bench; bench.trials(ntrials); bench.iterations(iterations); - bench.units(Benchmarker::Unit::ms); + bench.units(Benchmarker::Unit::us); ImageSpec spec_hd_rgba_f(1920, 1080, 4, TypeFloat); ImageSpec spec_hd_rgba_u8(1920, 1080, 4, TypeUInt8); @@ -771,6 +889,12 @@ test_resample() [&]() { ImageBufAlgo::resample(smallu8, buf_hd_rgba_f, false); }); bench(" IBA::resample HD->1024x512 rgba u8->u8 no interp ", [&]() { ImageBufAlgo::resample(smallu8, buf_hd_rgba_u8, false); }); + + // A nearest magnification: the reductions above cannot show the case + // where consecutive destination pixels read the same source pixel. + ImageBuf bigu8(ImageSpec(3840, 2160, 4, TypeUInt8)); + bench(" IBA::resample HD->4K rgba u8->u8 no interp ", + [&]() { ImageBufAlgo::resample(bigu8, buf_hd_rgba_u8, false); }); } @@ -1811,6 +1935,7 @@ main(int argc, char** argv) test_over(TypeFloat); test_over(TypeHalf); test_zover(); + test_resample_correctness(); test_resample(); test_compare(); test_isConstantColor(); diff --git a/src/libOpenImageIO/imagebufalgo_xform.cpp b/src/libOpenImageIO/imagebufalgo_xform.cpp index b469da6173..16ec4d8c67 100644 --- a/src/libOpenImageIO/imagebufalgo_xform.cpp +++ b/src/libOpenImageIO/imagebufalgo_xform.cpp @@ -8,6 +8,7 @@ #include #include +#include #include @@ -1266,6 +1267,233 @@ resample_deep(ImageBuf& dst, const ImageBuf& src, bool interpolate, ROI roi, +// A resampling filter's destination-to-source mapping is axis-separable and +// monotonically increasing, so it can be tabulated per axis instead of +// evaluated per pixel, and the destination range that lands inside the source +// data window is a single contiguous interval instead of a per-pixel range +// test. `pos` holds the source position of the filter for destination index +// `begin + i`; the entries are valid over the half-open range [begin, end), +// and outside it the destination is black. +// +// This carries no filter-specific state, so a wider filter can reuse it by +// tabulating its leftmost tap here and pairing it with a matching table of +// weights. +// For a two-tap filter `pos1` holds the second tap and `frac` its weight; both +// are empty for a single-tap filter. The two tables are only adjacent in the +// interior -- at the edges the clamp collapses them onto the same pixel -- so +// the second tap is stored rather than derived, which keeps the sampling loop +// branch-free right across the border. +struct AxisMap { + std::vector pos; + std::vector pos1; + std::vector frac; + int begin = 0; + int end = 0; +}; + + + +static AxisMap +build_axis_map(int dst_begin, int dst_end, float dst_full_origin, + float dst_full_size, float src_full_origin, float src_full_size, + int src_data_begin, int src_data_end, int stride) +{ + AxisMap map; + map.begin = dst_begin; + float dst_pixel_size = 1.0f / dst_full_size; + map.pos.reserve(size_t(dst_end - dst_begin)); + for (int d = dst_begin; d < dst_end; ++d) { + float t = (d - dst_full_origin + 0.5f) * dst_pixel_size; + int s = ifloor(src_full_origin + t * src_full_size); + if (s < src_data_begin) + continue; + if (s >= src_data_end) + break; + if (map.pos.empty()) + map.begin = d; + map.pos.push_back((s - src_data_begin) * stride); + } + map.end = map.begin + int(map.pos.size()); + return map; +} + + + +static AxisMap +build_axis_map_bilinear(int dst_begin, int dst_end, float dst_full_origin, + float dst_full_size, float src_full_origin, + float src_full_size, int src_data_begin, + int src_data_end, int stride) +{ + AxisMap map; + map.begin = dst_begin; + map.end = dst_end; + float dst_pixel_size = 1.0f / dst_full_size; + size_t n = size_t(dst_end - dst_begin); + map.pos.reserve(n); + map.pos1.reserve(n); + map.frac.reserve(n); + for (int d = dst_begin; d < dst_end; ++d) { + // Kept as separate statements mirroring resample_scalar() plus the + // half-pixel shift interppixel() applies internally: folding them into + // one expression lets the compiler contract the multiply-add, which + // shifts results by an ulp and lands on the wrong side of an integer + // rounding boundary often enough to matter. + float t = (d - dst_full_origin + 0.5f) * dst_pixel_size; + float src_xf = src_full_origin + t * src_full_size; + float pos = src_xf - 0.5f; + int s; + float f = floorfrac(pos, &s); + int s0 = clamp(s, src_data_begin, src_data_end - 1); + int s1 = clamp(s + 1, src_data_begin, src_data_end - 1); + map.pos.push_back((s0 - src_data_begin) * stride); + map.pos1.push_back((s1 - src_data_begin) * stride); + map.frac.push_back(f); + } + return map; +} + + + +static bool +data_window_is_full_window(const ImageSpec& spec) +{ + return spec.x == spec.full_x && spec.y == spec.full_y + && spec.width == spec.full_width && spec.height == spec.full_height; +} + + + +template +static bool +axis_map_usable(const ImageBuf& src, const ROI& roi) +{ + return src.localpixels() && src.contiguous_scanline() + && src.spec().format.basetype == BaseTypeFromC::value + && src.spec().depth == 1 && roi.chend <= src.spec().nchannels + && src.scanline_stride() % stride_t(sizeof(SRCTYPE)) == 0; +} + + + +template +static bool +resample_nearest(ImageBuf& dst, const ImageBuf& src, ROI roi, int nthreads) +{ + const ImageSpec& srcspec(src.spec()); + const ImageSpec& dstspec(dst.spec()); + int src_xstride = srcspec.nchannels; + int src_ystride = int(src.scanline_stride() / stride_t(sizeof(SRCTYPE))); + + // Built here rather than inside the parallel region: the maps depend only + // on the full windows, so every thread would otherwise rebuild the same + // tables, and allocating in a parallel region is worth avoiding anyway. + AxisMap xmap = build_axis_map(roi.xbegin, roi.xend, float(dstspec.full_x), + float(dstspec.full_width), + float(srcspec.full_x), + float(srcspec.full_width), srcspec.x, + srcspec.x + srcspec.width, src_xstride); + AxisMap ymap = build_axis_map(roi.ybegin, roi.yend, float(dstspec.full_y), + float(dstspec.full_height), + float(srcspec.full_y), + float(srcspec.full_height), srcspec.y, + srcspec.y + srcspec.height, src_ystride); + + const SRCTYPE* srcpixels = (const SRCTYPE*)src.localpixels(); + + ImageBufAlgo::parallel_image(roi, nthreads, [&](ROI roi) { + ImageBuf::Iterator out(dst, roi); + for (int y = roi.ybegin; y < roi.yend; ++y) { + bool row_inside = y >= ymap.begin && y < ymap.end; + const SRCTYPE* srcrow = row_inside + ? srcpixels + ymap.pos[y - ymap.begin] + : nullptr; + // When the row itself is outside the source, the whole scanline is + // black, which is the same as an empty inside-range. + int inside_begin = row_inside ? xmap.begin : roi.xend; + int inside_end = row_inside ? xmap.end : roi.xend; + int x = roi.xbegin; + for (; x < inside_begin; ++x, ++out) + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = 0.0f; + for (; x < inside_end; ++x, ++out) { + const SRCTYPE* p = srcrow + xmap.pos[x - xmap.begin]; + // convert_type is what the iterator's operator[] applies, and + // it rescales integer types to [0,1]; reading the raw value + // instead would silently change every integer image. + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = convert_type(p[c]); + } + for (; x < roi.xend; ++x, ++out) + for (int c = roi.chbegin; c < roi.chend; ++c) + out[c] = 0.0f; + } + }); + return true; +} + + + +template +static bool +resample_bilinear(ImageBuf& dst, const ImageBuf& src, ROI roi, int nthreads) +{ + const ImageSpec& srcspec(src.spec()); + const ImageSpec& dstspec(dst.spec()); + int src_xstride = srcspec.nchannels; + int src_ystride = int(src.scanline_stride() / stride_t(sizeof(SRCTYPE))); + + AxisMap xmap = build_axis_map_bilinear( + roi.xbegin, roi.xend, float(dstspec.full_x), float(dstspec.full_width), + float(srcspec.full_x), float(srcspec.full_width), srcspec.x, + srcspec.x + srcspec.width, src_xstride); + AxisMap ymap = build_axis_map_bilinear( + roi.ybegin, roi.yend, float(dstspec.full_y), float(dstspec.full_height), + float(srcspec.full_y), float(srcspec.full_height), srcspec.y, + srcspec.y + srcspec.height, src_ystride); + + const SRCTYPE* srcpixels = (const SRCTYPE*)src.localpixels(); + + ImageBufAlgo::parallel_image(roi, nthreads, [&](ROI roi) { + int nc = roi.chend - roi.chbegin; + // One scratch block per task, outside the pixel loop: the four taps + // have to be gathered into float before bilerp() sees them. + float* p0 = OIIO_ALLOCA(float, size_t(nc) * 5); + float* p1 = p0 + nc; + float* p2 = p0 + 2 * nc; + float* p3 = p0 + 3 * nc; + float* result = p0 + 4 * nc; + ImageBuf::Iterator out(dst, roi); + for (int y = roi.ybegin; y < roi.yend; ++y) { + size_t yi = size_t(y - ymap.begin); + const SRCTYPE* rowtop = srcpixels + ymap.pos[yi]; + const SRCTYPE* rowbot = srcpixels + ymap.pos1[yi]; + float wy = ymap.frac[yi]; + for (int x = roi.xbegin; x < roi.xend; ++x, ++out) { + size_t xi = size_t(x - xmap.begin); + int xl = xmap.pos[xi]; + int xr = xmap.pos1[xi]; + float wx = xmap.frac[xi]; + for (int c = 0; c < nc; ++c) { + int ch = roi.chbegin + c; + p0[c] = convert_type(rowtop[xl + ch]); + p1[c] = convert_type(rowtop[xr + ch]); + p2[c] = convert_type(rowbot[xl + ch]); + p3[c] = convert_type(rowbot[xr + ch]); + } + // The same bilerp() resample_scalar() uses, so the weights are + // associated identically and the two agree bit for bit. + bilerp(p0, p1, p2, p3, wx, wy, nc, result); + for (int c = 0; c < nc; ++c) + out[roi.chbegin + c] = result[c]; + } + } + }); + return true; +} + + + #if OIIO_USE_HWY template static bool @@ -1441,6 +1669,13 @@ resample_(ImageBuf& dst, const ImageBuf& src, bool interpolate, ROI roi, nthreads); #endif + if (axis_map_usable(src, roi) + && (!interpolate || data_window_is_full_window(src.spec()))) + return interpolate ? resample_bilinear(dst, src, roi, + nthreads) + : resample_nearest(dst, src, roi, + nthreads); + return resample_scalar(dst, src, interpolate, roi, nthreads); } diff --git a/testsuite/oiiotool-xform/ref/out.txt b/testsuite/oiiotool-xform/ref/out.txt index ec2b6f66b6..ebb6137c7c 100644 --- a/testsuite/oiiotool-xform/ref/out.txt +++ b/testsuite/oiiotool-xform/ref/out.txt @@ -1,5 +1,7 @@ Comparing "resample.tif" and "ref/resample.tif" PASS +Comparing "resample-nearest.tif" and "ref/resample-nearest.tif" +PASS Comparing "resize.tif" and "ref/resize.tif" PASS Comparing "resize2.tif" and "ref/resize2.tif" diff --git a/testsuite/oiiotool-xform/ref/resample-nearest.tif b/testsuite/oiiotool-xform/ref/resample-nearest.tif new file mode 100644 index 0000000000..e68dde266b Binary files /dev/null and b/testsuite/oiiotool-xform/ref/resample-nearest.tif differ diff --git a/testsuite/oiiotool-xform/run.py b/testsuite/oiiotool-xform/run.py index f662e15462..f7f4eb684f 100755 --- a/testsuite/oiiotool-xform/run.py +++ b/testsuite/oiiotool-xform/run.py @@ -40,6 +40,10 @@ def make_test_pattern1 (filename, xres=288, yres=216) : # test resample command += oiiotool ("../common/grid.tif --resample 128x128 -o resample.tif") +# test resample with interpolation off (nearest). Separate from the case above +# because the two use entirely different code paths. +command += oiiotool ("../common/grid.tif --resample:interp=0 128x128 -o resample-nearest.tif") + # test resize command += oiiotool ("../common/grid.tif --resize 256x256 -o resize.tif") command += oiiotool ("../common/grid.tif --resize 25% -o resize2.tif") @@ -141,7 +145,7 @@ def make_test_pattern1 (filename, xres=288, yres=216) : # Outputs to check against references outputs = [ - "resample.tif", "resize.tif", "resize2.tif", + "resample.tif", "resample-nearest.tif", "resize.tif", "resize2.tif", "resize64.tif", "resize512.tif", "resized-offset.exr", "resizefrom.tif", "resizefromto.tif", "resizefromtooffset.tif",