diff --git a/src/nyx/features/contour.cpp b/src/nyx/features/contour.cpp index df318595..6fdfc3a9 100644 --- a/src/nyx/features/contour.cpp +++ b/src/nyx/features/contour.cpp @@ -707,7 +707,9 @@ void ContourFeature::buildRegularContour_nontriv (LR& r, const Fsettings& s) //==== Prepare the contour image WriteImageMatrix_nontriv borderImage("borderImage", r.label); - paddedImage.allocate(width + 2, height + 2, 0); + // allocate() writes from the file's current position, so it is called on the fresh border image + // only; on paddedImage it would zero every slot after the last pixel written above + borderImage.allocate(width + 2, height + 2, 0); // Set entire image to WHITE for (int y = 0; y < (height + 2); y++) diff --git a/src/nyx/phase3.cpp b/src/nyx/phase3.cpp index cc56a63b..6e2fc5d5 100644 --- a/src/nyx/phase3.cpp +++ b/src/nyx/phase3.cpp @@ -144,13 +144,28 @@ namespace Nyxus r.aux_max = 0; // Iterate ROI's tiles and scan pixels - size_t nth = env.theImLoader.get_num_tiles_hor(), - ntv = env.theImLoader.get_num_tiles_vert(); - for (unsigned int row = 0; row < nth; row++) - for (unsigned int col = 0; col < ntv; col++) + size_t ntHor = env.theImLoader.get_num_tiles_hor(), // tiles across a row + ntVert = env.theImLoader.get_num_tiles_vert(), // tiles down a column + th = env.theImLoader.get_tile_height(), + tw = env.theImLoader.get_tile_width(); + for (size_t row = 0; row < ntVert; row++) + for (size_t col = 0; col < ntHor; col++) { - unsigned int tileIdx = row * ntv + col; - env.theImLoader.load_tile(tileIdx); + // row-major, as the loader splits an index into tile_idx / ntHor and tile_idx % ntHor + size_t tileIdx = row * ntHor + col; + // load_tile() refuses only an index outside the grid; a tile the reader fails on + // arrives zero-filled and is not reported here + if (! env.theImLoader.load_tile(tileIdx)) + { + std::string erm = "Error fetching tile row=" + std::to_string(row) + " col=" + std::to_string(col); + r.raw_pixels_NT.clear(); + env.theImLoader.close(); +#ifdef WITH_PYTHON_H + throw std::runtime_error (erm); +#endif + std::cerr << erm << "\n"; + return false; + } auto& dataI = env.theImLoader.get_int_tile_buffer(); auto& dataL = env.theImLoader.get_seg_tile_buffer(); for (unsigned long i = 0; i < env.theImLoader.get_tile_size(); i++) @@ -163,12 +178,8 @@ namespace Nyxus // Pixel intensity and global position auto intens = dataI[i]; - size_t row = tileIdx / env.theImLoader.get_num_tiles_hor(), - col = tileIdx / env.theImLoader.get_num_tiles_hor(), - th = env.theImLoader.get_tile_height(), - tw = env.theImLoader.get_tile_width(); - int y = row * th + i / tw, - x = col * tw + i % tw; + StatsInt y = (StatsInt) (row * th + i / tw), + x = (StatsInt) (col * tw + i % tw); // Track the streamed extrema (see the aux_min/aux_max reset above) r.aux_min = (std::min)(r.aux_min, (PixIntens)intens); diff --git a/tests/python/test_2d_ooc_invariant.py b/tests/python/test_2d_ooc_invariant.py index 75ed4520..56311af5 100644 --- a/tests/python/test_2d_ooc_invariant.py +++ b/tests/python/test_2d_ooc_invariant.py @@ -400,3 +400,76 @@ def test_2d_ooc_2d_two_rois_no_state_leak_invariant(tmp_path): if np.isnan(p) != np.isnan(q) or abs(p - q) > 1e-6 * max(abs(p), abs(q), 1.0) + 1e-9: bad.append((lab, c, p, q)) assert not bad, "Zernike/erosion/ellipse out-of-core values diverge from in-RAM across ROIs: %r" % (bad[:8],) + + +# A 160 x 96 image at 32 x 32 tiles: 5 tiles across, 3 down. The grid is not square, so a scan that +# swaps the horizontal and vertical tile counts, or derives a tile's column from the wrong one, puts +# pixels at different coordinates than the untiled copy does. +TILED_W, TILED_H, TILE = 160, 96, 32 + + +def _make_tiled_untiled_pairs(tmp_path): + """The same intensity and mask arrays written twice: as tiled TIFFs and as strip TIFFs. + + Label 1 is an ellipse inside the single tile at grid row 0, column 3, off the grid's diagonal. + Label 2 spans tiles: a column at x 36..50 reaching down to y 60, and an arm to its right at + y 40..50 that runs to x 120. Tile by tile, the arm's last pixel is the last one streamed, and the + column's rows below the arm come before it, which is the reverse of their raster order.""" + yy, xx = np.mgrid[0:TILED_H, 0:TILED_W] + inten = (1 + xx + 7 * yy).astype(np.uint32) + mask = np.zeros((TILED_H, TILED_W), np.uint32) + mask[((xx - 112) / 12.0) ** 2 + ((yy - 15) / 9.0) ** 2 <= 1.0] = 1 + mask[4:61, 36:51] = 2 + mask[40:51, 51:121] = 2 + assert set(np.unique(mask[0:TILE, 3 * TILE:4 * TILE])) == {0, 1} + assert not (mask == 1)[:, :3 * TILE].any() and not (mask == 1)[TILE:, :].any() + + dirs = {} + for kind, kw in (("tiled", {"tile": (TILE, TILE)}), ("untiled", {})): + intdir = tmp_path / kind / "int" + segdir = tmp_path / kind / "seg" + intdir.mkdir(parents=True) + segdir.mkdir(parents=True) + tifffile.imwrite(str(intdir / "img.tif"), inten, **kw) + tifffile.imwrite(str(segdir / "img.tif"), mask, **kw) + with tifffile.TiffFile(str(intdir / "img.tif")) as tf: + assert tf.pages[0].is_tiled == (kind == "tiled") + dirs[kind] = (str(intdir) + os.sep, str(segdir) + os.sep) + return dirs + + +def test_2d_ooc_2d_tiled_matches_untiled_invariant(tmp_path): + """Out-of-core featurization of a tiled TIFF must equal that of a strip TIFF holding the same + arrays, every column of both ROIs. + + Both sides run out-of-core, so the pixels and the code path are the same and only the order in + which the pixel cloud is filled differs: tile by tile against row by row. That isolates the + out-of-core tile scan from any out-of-core vs in-RAM difference. What it discriminates: a scan + that places a tile's pixels at the wrong column fails on label 1, and a contour that depends on + the order the cloud was filled in fails on label 2.""" + dirs = _make_tiled_untiled_pairs(tmp_path) + + frames = {} + for kind, (intdir, segdir) in dirs.items(): + nyx = nyxus.Nyxus(ALL_2D_FEATURE_GROUPS) + _set_ram_limit_mb(nyx, 0) # 0 -> every ROI oversized -> out-of-core + frames[kind] = nyx.featurize_directory(intdir, segdir).sort_values("ROI_label").reset_index(drop=True) + + df_t, df_u = frames["tiled"], frames["untiled"] + assert list(df_t["ROI_label"]) == [1, 2] and list(df_u["ROI_label"]) == [1, 2] + cols, _ = _feature_cols(df_u) + assert cols and _feature_cols(df_t)[0] == cols + + # Tolerance: same build, same pixels, same formulas; only the summation order differs. The + # absolute floor covers moments that are zero analytically and come out as float noise. + bad = [] + for row in range(df_u.shape[0]): + lab = int(df_u["ROI_label"].iloc[row]) + for c in cols: + p = float(df_t[c].iloc[row]) + q = float(df_u[c].iloc[row]) + if np.isnan(p) and np.isnan(q): + continue + if np.isnan(p) != np.isnan(q) or abs(p - q) > 1e-6 * max(abs(p), abs(q), 1.0) + 1e-9: + bad.append((lab, c, p, q)) + assert not bad, "out-of-core features of a tiled TIFF diverge from its untiled copy: %r" % (bad[:8],)