Dev #19

Merged
zhangyijian merged 53 commits from dev into main 2026-09-02 21:39:04 +08:00
64 changed files with 527 additions and 7591 deletions
Showing only changes of commit 28e6c7c8a9 - Show all commits
+5 -1
View File
@@ -20,7 +20,7 @@ from paths import NUITKA_COMPILE_DIR, NUITKA_WHEEL_DIR, REPO_ROOT, repo_path
from pyproject_utils import bundle_items, extras_require, get_version, install_requires
from verify import verify_built_wheel
from wheel_build import publish_wheel_to_sdk, run_bdist_wheel, write_wheel_setup
from wheel_layout import copy_bundle_items, package_data_from_bundle
from wheel_layout import copy_bundle_items, copy_orisys_blend, package_data_from_bundle
def _reset_dir(path) -> None:
@@ -51,6 +51,10 @@ def build() -> None:
ext_suffix = copy_compiled_orisys(pkg_dir, ext_file)
wheel_root = NUITKA_WHEEL_DIR / "wheel_root"
blend_count = copy_orisys_blend(REPO_ROOT, wheel_root)
print(f"Staged orisys_blend Python modules: {blend_count}")
# Bundle rules live in pyproject.toml → [[tool.orisys.bundle.items]]
bundle = bundle_items(REPO_ROOT)
copied_counts = copy_bundle_items(REPO_ROOT, pkg_dir, bundle)
+1
View File
@@ -19,6 +19,7 @@ NUITKA_INCLUDE_PACKAGES = (
)
NUITKA_NOFOLLOW_IMPORTS = (
"orisys_blend",
"cupy",
"cupyx",
"pyopencl",
+1 -1
View File
@@ -31,7 +31,7 @@ setup(
version="{version}",
description="Orisys - Optical Tactile Sensor SDK (Nuitka Compiled)",
author="Orisys",
packages=['orisys', 'orisys.examples'],
packages=['orisys', 'orisys.examples', 'orisys_blend'],
package_dir={{'': 'wheel_root'}},
package_data={package_data},
include_package_data=True,
+13
View File
@@ -83,9 +83,22 @@ def package_data_from_bundle(ext_suffix: str, items: list[dict]) -> dict[str, li
package_data: dict[str, list[str]] = {
"orisys": [f"*{ext_suffix}"],
"orisys.examples": ["*.py", "README.md", "viz/**/*.py", "pipeline/**/*.py", "gui/**/*.py"],
"orisys_blend": ["*.py"],
}
for item in items:
if item["name"] == "examples":
continue
package_data["orisys"].extend(package_data_paths_for_item(item))
return package_data
def copy_orisys_blend(repo_root: Path, wheel_root: Path) -> int:
"""Copy pure-Python orisys_blend (Numba kernels) next to the Nuitka orisys .pyd."""
src = repo_root / "src" / "orisys_blend"
dst = wheel_root / "orisys_blend"
if dst.exists():
shutil.rmtree(dst)
shutil.copytree(src, dst)
count = sum(1 for p in dst.rglob("*.py") if p.is_file())
print(f"Copied orisys_blend ({count} .py files) to {dst}")
return count
+3 -1
View File
@@ -18,7 +18,9 @@
"nominal_ov_pos": 100,
"n": [3,3],
"row_tol": 0.6,
"box_size": [250, 250]
"box_size": [250, 250],
"blend":"numba",
"blend_list": ["numpy", "numba","none"]
},
"tracker": {
"adaptive": true,
+3 -1
View File
@@ -18,7 +18,9 @@
"nominal_ov_pos": 44,
"n": [3,3],
"row_tol": 0.6,
"box_size": [98, 98]
"box_size": [98, 98],
"blend":"numpy",
"blend_list": ["numpy", "numba","none"]
},
"tracker": {
"adaptive": true,
+3 -1
View File
@@ -18,7 +18,9 @@
"nominal_ov_pos": 67,
"n": [3,3],
"row_tol": 0.6,
"box_size": [180, 180]
"box_size": [180, 180],
"blend":"numpy",
"blend_list": ["numpy", "numba","none"]
},
"tracker": {
"adaptive": true,
+3 -1
View File
@@ -18,7 +18,9 @@
"nominal_ov_pos": 67,
"n": [3,3],
"row_tol": 0.6,
"box_size": [180, 180]
"box_size": [180, 180],
"blend":"numpy",
"blend_list": ["numpy", "numba","none"]
},
"tracker": {
"adaptive": true,
-44
View File
@@ -1,44 +0,0 @@
{
"name": "prototype",
"version": "1.0",
"camera": {
"resolution_raw": [1920, 1080],
"fps": 30,
"units": 15,
"units_used": 9,
"resolution_stitch": [365, 365],
"stitch_res_list": [365,263,123,63]
},
"stitcher": {
"min_area_frac": 0.01,
"max_area_frac": 0.04,
"max_aspect_ratio": 100,
"index": [1,2,3,6,7,8,11,12,13],
"margin": 10,
"nominal_ov_pos": 100,
"n": [3,3],
"row_tol": 0.6,
"box_size": [220, 220]
},
"tracker": {
"adaptive": true,
"lower_threshold": 12.0,
"upper_threshold": 18.0,
"loop_closure_interval": 1,
"algorithm": "dis",
"dis_preset": "medium",
"dis_mode": 3,
"algorithm_list": ["raft", "dis"],
"hhd_downsample": 3
},
"contact": {
"min_contour_area": 3000,
"flow_threshold":1,
"div_threshold":0.02,
"expansion_threshold": 0.01,
"kernel_size": [5,5],
"denoise_sigma": 15
}
}
-44
View File
@@ -1,44 +0,0 @@
{
"name": "prototype",
"version": "1.0",
"camera": {
"resolution_raw": [854, 480],
"fps": 30,
"units": 15,
"units_used": 9,
"resolution_stitch": [123, 123],
"stitch_res_list": [365,263,123,63]
},
"stitcher": {
"min_area_frac": 0.01,
"max_area_frac": 0.06,
"max_aspect_ratio": 100,
"index": [1,2,3,6,7,8,11,12,13],
"margin": 4,
"nominal_ov_pos": 44,
"n": [3,3],
"row_tol": 0.6,
"box_size": [98, 98]
},
"tracker": {
"adaptive": true,
"lower_threshold": 12.0,
"upper_threshold": 18.0,
"loop_closure_interval": 1,
"algorithm": "dis",
"dis_preset": "medium",
"dis_mode": 2,
"algorithm_list": ["raft", "dis"],
"hhd_downsample": 3
},
"contact": {
"min_contour_area": 3000,
"flow_threshold":1,
"div_threshold":0.02,
"expansion_threshold": 0.01,
"kernel_size": [5,5],
"denoise_sigma": 6
}
}
-44
View File
@@ -1,44 +0,0 @@
{
"name": "prototype",
"version": "1.0",
"camera": {
"resolution_raw": [1280, 720],
"fps": 30,
"units": 15,
"units_used": 9,
"resolution_stitch": [263, 263],
"stitch_res_list": [365,263,123,63]
},
"stitcher": {
"min_area_frac": 0.01,
"max_area_frac": 0.04,
"max_aspect_ratio": 25,
"index": [1,2,3,6,7,8,11,12,13],
"margin": 10,
"nominal_ov_pos": 67,
"n": [3,3],
"row_tol": 0.6,
"box_size": [180, 180]
},
"tracker": {
"adaptive": true,
"lower_threshold": 12.0,
"upper_threshold": 18.0,
"loop_closure_interval": 1,
"algorithm": "dis",
"dis_preset": "medium",
"dis_mode": 2,
"algorithm_list": ["raft", "dis"],
"hhd_downsample": 3
},
"contact": {
"min_contour_area": 3000,
"flow_threshold":1,
"div_threshold":0.02,
"expansion_threshold": 0.01,
"kernel_size": [5,5],
"denoise_sigma": 10
}
}
-43
View File
@@ -1,43 +0,0 @@
{
"name": "prototype",
"version": "1.0",
"camera": {
"resolution_raw": [1920, 1080],
"fps": 30,
"units": 15,
"units_used": 9,
"resolution_stitch": [384, 384]
},
"stitcher": {
"min_area_frac": 0.01,
"max_area_frac": 0.04,
"max_aspect_ratio": 100,
"index": [1,2,3,6,7,8,11,12,13],
"margin": 10,
"nominal_ov_pos": 100,
"n": [3,3],
"row_tol": 0.6,
"box_size": [220, 220]
},
"tracker": {
"adaptive": true,
"lower_threshold": 12.0,
"upper_threshold": 18.0,
"loop_closure_interval": 1,
"algorithm": "raft",
"dis_preset": "medium",
"dis_mode": 3,
"algorithm_list": ["raft", "dis"],
"hhd_downsample": 3
},
"contact": {
"min_contour_area": 3000,
"flow_threshold":1.5,
"div_threshold":0.08,
"expansion_threshold": 0.01,
"kernel_size": [5,5],
"denoise_sigma": 15
}
}
+1
View File
@@ -22,6 +22,7 @@ cuda13 = ["cupy-cuda13x>=13.0.0"]
# Backward compatibility: 'cuda' defaults to CUDA 12.x
cuda = ["cupy-cuda12x>=12.3.0"]
gui = ["pyqt5>=5.15", "pyqtgraph>=0.13"]
numba = ["numba>=0.59.0"]
[tool.setuptools]
package-dir = { "" = "src" }
+432 -118
View File
@@ -5,45 +5,92 @@ from scipy import ndimage
import cv2
try:
from numba import njit
from orisys_blend.kernel import _NUMBA_AVAILABLE, blend_kernel as _blend_kernel_numba
except Exception:
_NUMBA_AVAILABLE = False
_blend_kernel_numba = None
_NUMBA_AVAILABLE = True
except ImportError:
# Re-export availability from the pure-Python kernel package when import succeeded
# but numba itself was missing (blend_kernel is None).
if _blend_kernel_numba is None:
_NUMBA_AVAILABLE = False
def njit(*args, **kwargs):
def decorator(fn):
return fn
return decorator
def _blend_kernel_numpy(blend_patch, wsum_patch, cov_patch, src_patch, w_patch):
"""Incremental weighted-average composite (pure NumPy, no JIT).
Uses contiguous arithmetic instead of boolean fancy indexing. Callers must
zero ``blend`` / ``wsum`` on uncovered pixels (as the stitcher reset does);
then the covered lerp also paints first-touch pixels correctly when ``w > 0``.
"""
w = w_patch if w_patch.dtype == np.float32 else np.asarray(w_patch, dtype=np.float32)
src = src_patch if src_patch.dtype == np.float32 else np.asarray(src_patch, dtype=np.float32)
uncovered = cov_patch == 0
if uncovered.all():
# First paint over a cleared patch: direct copy.
np.copyto(blend_patch, src)
np.copyto(wsum_patch, w)
cov_patch.fill(1)
return
# Contiguous lerp. With zeroed uncovered pixels this also equals a first paint
# when w > 0 (a ≈ 1 → blend becomes src).
ws = wsum_patch + w
a = w / np.maximum(ws, np.float32(1e-6))
# Channel loop keeps temps HxW instead of one HxWx3 gather/scatter.
for c in range(src.shape[2]):
bc = blend_patch[..., c]
bc += (src[..., c] - bc) * a
np.copyto(wsum_patch, ws)
# Match numba on rare first-touch pixels with ~0 weight (lerp would leave 0).
if uncovered.any():
need_fix = uncovered & (w < np.float32(1e-6))
if need_fix.any():
blend_patch[need_fix] = src[need_fix]
wsum_patch[need_fix] = w[need_fix]
cov_patch.fill(1)
@njit(cache=False)
def _blend_kernel(blend_patch, wsum_patch, cov_patch, src_patch, w_patch):
"""Fused incremental weighted-average composite (numba-accelerated)."""
sh = src_patch.shape[0]
sw = src_patch.shape[1]
for y in range(sh):
for x in range(sw):
w = w_patch[y, x]
s0 = np.float32(src_patch[y, x, 0])
s1 = np.float32(src_patch[y, x, 1])
s2 = np.float32(src_patch[y, x, 2])
if cov_patch[y, x] == 0:
blend_patch[y, x, 0] = s0
blend_patch[y, x, 1] = s1
blend_patch[y, x, 2] = s2
wsum_patch[y, x] = w
else:
ws = wsum_patch[y, x] + w
if ws < 1e-6:
ws = 1e-6
a = w / ws
blend_patch[y, x, 0] += (s0 - blend_patch[y, x, 0]) * a
blend_patch[y, x, 1] += (s1 - blend_patch[y, x, 1]) * a
blend_patch[y, x, 2] += (s2 - blend_patch[y, x, 2]) * a
wsum_patch[y, x] = ws
cov_patch[y, x] = 1
# Backward-compatible alias used by older call sites / tests.
_blend_kernel = _blend_kernel_numpy
# JSON ``stitcher.blend`` / ``blend_list``: none = fast paste; numpy/numba = feathered.
_BLEND_MODES = frozenset({"none", "numpy", "numba"})
_DEFAULT_BLEND_MODE = "none"
def _resolve_blend_mode(config_stitcher=None):
"""Resolve stitcher blend mode from env, then JSON ``blend``.
Returns ``(mode, fell_back_from_numba)`` where mode is one of
``none`` / ``numpy`` / ``numba``.
Priority: ``ORISYS_STITCH_BLEND``, ``ORISYS_STITCH_BLEND_BACKEND``,
then ``config_stitcher["blend"]`` / legacy ``blend_backend``, then ``none``.
"""
requested = None
for key in ("ORISYS_STITCH_BLEND", "ORISYS_STITCH_BLEND_BACKEND"):
env = os.environ.get(key, "").strip().lower()
if env in _BLEND_MODES:
requested = env
break
if requested is None and isinstance(config_stitcher, dict):
for key in ("blend", "blend_backend"):
cfg = str(config_stitcher.get(key, "")).strip().lower()
if cfg in _BLEND_MODES:
requested = cfg
break
if requested is None:
requested = _DEFAULT_BLEND_MODE
if requested == "numba" and not _NUMBA_AVAILABLE:
return "numpy", True
return requested, False
class Stitcher:
@@ -85,25 +132,100 @@ class Stitcher:
self._feature_stitch_print_every = 30
self._feature_stitch_pad = 5
self._feature_stitch_match_brightness = False
self._blend_mode, fell_back = _resolve_blend_mode(config_stitcher)
# Feather path uses numpy/numba; ``none`` keeps a numpy fallback if forced on.
self._blend_backend = "numpy" if self._blend_mode == "none" else self._blend_mode
self._numba_failed = False
if self.isstitch and not _NUMBA_AVAILABLE:
if fell_back and self.isstitch:
print(
"Warning: numba is not installed; feature stitching will run a slow "
"Python blend path (~1 FPS). Install numba for real-time performance."
"Warning: stitcher blend=numba requested but numba is "
"unavailable; using numpy blend path."
)
elif self.isstitch and self.verbose:
print(f"Stitch blend: {self._blend_mode}")
if self.verbose:
print(f"Stitcher initialized with config: {self.config_stitcher}")
def configure_feature_stitch(self, *, match_brightness=None, pad=None, print_every=None):
"""Update runtime options for feature-based stitching."""
@property
def blend_mode(self):
"""Configured blend mode: ``none``, ``numpy``, or ``numba``."""
return self._blend_mode
def configure_feature_stitch(
self,
*,
match_brightness=None,
pad=None,
print_every=None,
blend=None,
blend_backend=None,
):
"""Update runtime options for feature-based stitching.
``blend`` accepts ``none`` / ``numpy`` / ``numba`` (same as JSON ``stitcher.blend``).
``blend_backend`` is a legacy alias for the same values.
"""
if match_brightness is not None:
self._feature_stitch_match_brightness = bool(match_brightness)
if pad is not None:
self._feature_stitch_pad = max(0, int(pad))
if print_every is not None:
self._feature_stitch_print_every = max(1, int(print_every))
mode = blend if blend is not None else blend_backend
if mode is not None:
mode = str(mode).strip().lower()
if mode not in _BLEND_MODES:
raise ValueError(
f"Unknown blend={mode!r}; expected one of {sorted(_BLEND_MODES)}"
)
if mode == "numba" and (not _NUMBA_AVAILABLE or self._numba_failed):
print("Warning: numba unavailable; using numpy blend.")
mode = "numpy"
self._blend_mode = mode
self._blend_backend = "numpy" if mode == "none" else mode
if mode == "numba":
self._numba_failed = False
def _resolve_call_blend(self, blend):
"""Map a ``stitch_rt_optimized(blend=...)`` arg to ``(use_feather, backend)``."""
if blend is None:
mode = self._blend_mode
elif isinstance(blend, bool):
if blend:
mode = self._blend_mode if self._blend_mode != "none" else "numpy"
else:
return False, self._blend_backend
else:
mode = str(blend).strip().lower()
if mode not in _BLEND_MODES:
raise ValueError(
f"Unknown blend={mode!r}; expected one of {sorted(_BLEND_MODES)}"
)
if mode == "numba" and not _NUMBA_AVAILABLE:
mode = "numpy"
if mode == "numba" and self._numba_failed:
mode = "numpy"
if mode == "none":
return False, self._blend_backend
return True, mode
def _mark_numba_failed(self, exc):
"""Stick to numpy after a Numba JIT/runtime failure (log once)."""
if self._numba_failed:
return
self._numba_failed = True
self._blend_backend = "numpy"
if self._blend_mode == "numba":
self._blend_mode = "numpy"
print(
f"Warning: numba blend failed ({type(exc).__name__}: {exc}); "
f"falling back to numpy."
)
def _get_blend_coef(self, ov_pos, dim, ndim):
key = (ov_pos, dim, ndim)
coef = self._blend_coef_cache.get(key)
@@ -117,13 +239,30 @@ class Stitcher:
self._blend_coef_cache[key] = blend_coef
return blend_coef
def stitch_rt_optimized(self, change_sv=False, s_mean=28., v_mean=180., blend=True):
feature_stitched = self.stitch_feature_optimized(
blend=blend,
match_brightness=self._feature_stitch_match_brightness,
pad=self._feature_stitch_pad,
print_every=self._feature_stitch_print_every,
)
def stitch_rt_optimized(self, change_sv=False, s_mean=28., v_mean=180., blend=None):
"""Stitch the current tile list.
``blend`` overrides JSON/env for this call:
- ``None`` (default): use configured ``stitcher.blend``
- ``True`` / ``False``: force feathered blend on/off
- ``\"none\"`` / ``\"numpy\"`` / ``\"numba\"``: force that mode for this call
"""
use_feather, blend_backend = self._resolve_call_blend(blend)
prev_backend = self._blend_backend
self._blend_backend = blend_backend
try:
feature_stitched = self.stitch_feature_optimized(
blend=use_feather,
match_brightness=self._feature_stitch_match_brightness,
pad=self._feature_stitch_pad,
print_every=self._feature_stitch_print_every,
)
finally:
# Keep sticky numpy fallback after a Numba failure; do not restore numba.
if self._numba_failed:
self._blend_backend = "numpy"
else:
self._blend_backend = prev_backend
if feature_stitched is not None:
return feature_stitched
@@ -455,7 +594,7 @@ class Stitcher:
return wy[:, None] * wx[None, :]
def _blend_tile_inplace(self, blend, wsum, cov, src, weight, y0, x0, scratch=None):
"""Composite a tile onto the canvas using numba-fused weighted average."""
"""Composite a tile onto the canvas using weighted average blending."""
h, w = src.shape[:2]
H, W = blend.shape[:2]
y0d, x0d = max(0, y0), max(0, x0)
@@ -467,13 +606,20 @@ class Stitcher:
sh, sw = y1d - y0d, x1d - x0d
src_f = src[sy0:sy0+sh, sx0:sx0+sw].astype(np.float32, copy=False)
_blend_kernel(
args = (
blend[y0d:y1d, x0d:x1d],
wsum[y0d:y1d, x0d:x1d],
cov[y0d:y1d, x0d:x1d],
src_f,
weight[sy0:sy0+sh, sx0:sx0+sw],
)
if self._blend_backend == "numba" and _blend_kernel_numba is not None:
try:
_blend_kernel_numba(*args)
return
except Exception as exc:
self._mark_numba_failed(exc)
_blend_kernel_numpy(*args)
def _mean_gray_bgr(self, bgr):
ch = cv2.mean(bgr)
@@ -535,6 +681,214 @@ class Stitcher:
shapes = tuple(tb.shape[:2] for tb in tiles_bgr_c)
return ref_h, ref_w, shapes, pad
def _ensure_feature_geom_cache(
self,
canvas_cache,
t_ref_bgr,
tiles_bgr_c,
H_use_list,
pad,
*,
allocate_blend_bufs=True,
):
"""Populate/reuse fixed-window geometry in canvas_cache. Returns geom_cache_hit."""
h0c, w0c = t_ref_bgr.shape[:2]
geom_key = self._geom_cache_key(t_ref_bgr, tiles_bgr_c, pad)
geom_cache_hit = canvas_cache.get("geom_key") == geom_key
if geom_cache_hit:
if allocate_blend_bufs and (
"blend_buf" not in canvas_cache
or "wsum_buf" not in canvas_cache
or "cov_buf" not in canvas_cache
):
out_h = canvas_cache["out_h"]
out_w = canvas_cache["out_w"]
canvas_cache["blend_buf"] = np.zeros((out_h, out_w, 3), dtype=np.float32)
canvas_cache["wsum_buf"] = np.zeros((out_h, out_w), dtype=np.float32)
canvas_cache["cov_buf"] = np.zeros((out_h, out_w), dtype=np.uint8)
return True
ref_corners = np.array([[0, 0], [w0c, 0], [w0c, h0c], [0, h0c]], dtype=np.float64)
pts = [ref_corners]
for tile_bgr, H_use in zip(tiles_bgr_c, H_use_list):
if H_use is None:
continue
hi, wi = tile_bgr.shape[:2]
tx, ty = self._translation_from_H_to_ref(H_use)
pts.append(
np.array(
[[tx, ty], [tx + wi, ty], [tx + wi, ty + hi], [tx, ty + hi]],
dtype=np.float64,
)
)
all_xy = np.concatenate(pts, axis=0)
x_min_f = float(all_xy[:, 0].min()) - pad
y_min_f = float(all_xy[:, 1].min()) - pad
x_max_f = float(all_xy[:, 0].max()) + pad
y_max_f = float(all_xy[:, 1].max()) + pad
out_w = max(1, int(np.ceil(x_max_f - x_min_f)))
out_h = max(1, int(np.ceil(y_max_f - y_min_f)))
ref_y0 = int(np.floor(-y_min_f))
ref_x0 = int(np.floor(-x_min_f))
pastes = []
for H_use in H_use_list:
if H_use is None:
pastes.append(None)
continue
tx, ty = self._translation_from_H_to_ref(H_use)
pastes.append((int(np.floor(ty - y_min_f)), int(np.floor(tx - x_min_f))))
_dr_y0 = max(0, ref_y0)
_dr_x0 = max(0, ref_x0)
_dr_y1 = min(out_h, ref_y0 + h0c)
_dr_x1 = min(out_w, ref_x0 + w0c)
for _tb, _pa in zip(tiles_bgr_c, pastes):
if _pa is None:
continue
_py, _px = _pa
_th, _tw = _tb.shape[:2]
_dr_y0 = min(_dr_y0, max(0, _py))
_dr_x0 = min(_dr_x0, max(0, _px))
_dr_y1 = max(_dr_y1, min(out_h, _py + _th))
_dr_x1 = max(_dr_x1, min(out_w, _px + _tw))
canvas_cache["geom_key"] = geom_key
canvas_cache["out_w"] = out_w
canvas_cache["out_h"] = out_h
canvas_cache["ref_y0"] = ref_y0
canvas_cache["ref_x0"] = ref_x0
canvas_cache["pastes"] = pastes
canvas_cache["x_min_f"] = x_min_f
canvas_cache["y_min_f"] = y_min_f
canvas_cache["crop_rect"] = None
canvas_cache["dirty_rect"] = (_dr_y0, _dr_y1, _dr_x0, _dr_x1)
canvas_cache["weight_maps"] = {}
canvas_cache["uint8_buf"] = np.zeros((out_h, out_w, 3), dtype=np.uint8)
if allocate_blend_bufs:
canvas_cache["blend_buf"] = np.zeros((out_h, out_w, 3), dtype=np.float32)
canvas_cache["wsum_buf"] = np.zeros((out_h, out_w), dtype=np.float32)
canvas_cache["cov_buf"] = np.zeros((out_h, out_w), dtype=np.uint8)
else:
# Coverage buffer only needed once to compute crop_rect for paste path.
canvas_cache["cov_buf"] = np.zeros((out_h, out_w), dtype=np.uint8)
canvas_cache.pop("blend_buf", None)
canvas_cache.pop("wsum_buf", None)
return False
@staticmethod
def _paste_tile_clipped(canvas, tile, y0, x0, cov=None):
"""Paste tile into canvas with clipping; optionally mark coverage."""
th, tw = tile.shape[:2]
H, W = canvas.shape[:2]
y0d, x0d = max(0, y0), max(0, x0)
y1d, x1d = min(H, y0 + th), min(W, x0 + tw)
if y0d >= y1d or x0d >= x1d:
return
sy0, sx0 = y0d - y0, x0d - x0
sh, sw = y1d - y0d, x1d - x0d
np.copyto(canvas[y0d:y1d, x0d:x1d], tile[sy0:sy0 + sh, sx0:sx0 + sw])
if cov is not None:
cov[y0d:y1d, x0d:x1d] = 1
def _paste_feature_tiles_fast(
self,
t_ref_bgr,
tiles_bgr_c,
H_use_list,
canvas_cache,
*,
geom_cache_hit,
frame_idx=0,
print_every=30,
verbose=False,
t_profile0=0.0,
t_after_brightness=0.0,
t_after_geom=0.0,
):
"""Fixed-geometry uint8 paste path used when blend=False."""
do_profile = verbose and frame_idx % print_every == 0
out_w = canvas_cache["out_w"]
out_h = canvas_cache["out_h"]
ref_y0 = canvas_cache["ref_y0"]
ref_x0 = canvas_cache["ref_x0"]
pastes = canvas_cache["pastes"]
canvas = canvas_cache.get("uint8_buf")
if canvas is None or canvas.shape[0] != out_h or canvas.shape[1] != out_w:
canvas = np.zeros((out_h, out_w, 3), dtype=np.uint8)
canvas_cache["uint8_buf"] = canvas
dirty = canvas_cache.get("dirty_rect")
if geom_cache_hit and dirty is not None:
dy0, dy1, dx0, dx1 = dirty
canvas[dy0:dy1, dx0:dx1] = 0
elif not geom_cache_hit:
canvas.fill(0)
else:
canvas.fill(0)
t_after_reset = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
need_crop = canvas_cache.get("crop_rect") is None
cov = canvas_cache.get("cov_buf") if need_crop else None
if need_crop:
if cov is None or cov.shape != (out_h, out_w):
cov = np.zeros((out_h, out_w), dtype=np.uint8)
canvas_cache["cov_buf"] = cov
else:
if dirty is not None:
dy0, dy1, dx0, dx1 = dirty
cov[dy0:dy1, dx0:dx1] = 0
else:
cov.fill(0)
pasted_tiles = 0
for tile_bgr, H_use, paste in zip(tiles_bgr_c, H_use_list, pastes):
if H_use is None or paste is None:
continue
y0, x0 = paste
self._paste_tile_clipped(canvas, tile_bgr, y0, x0, cov=cov)
pasted_tiles += 1
self._paste_tile_clipped(canvas, t_ref_bgr, ref_y0, ref_x0, cov=cov)
t_after_paste = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
crop_rect = canvas_cache.get("crop_rect")
if crop_rect is None and cov is not None:
crop_rect = self._compute_final_crop_rect(cov)
canvas_cache["crop_rect"] = crop_rect
if crop_rect is not None:
y0, y1, x0, x1 = crop_rect
roi = canvas[y0:y1, x0:x1]
else:
roi = canvas
rh, rw = roi.shape[:2]
stitch_buf = canvas_cache.get("stitch_buf")
if stitch_buf is None or stitch_buf.shape[0] != rh or stitch_buf.shape[1] != rw:
stitch_buf = np.empty((rh, rw, 3), dtype=np.uint8)
canvas_cache["stitch_buf"] = stitch_buf
np.copyto(stitch_buf, roi)
t_after_finalize = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
if do_profile:
xmf = canvas_cache.get("x_min_f", 0.0)
ymf = canvas_cache.get("y_min_f", 0.0)
print(
f"stitch_feature_optimized: canvas out_w={out_w}, out_h={out_h}, "
f"x_min={xmf:.2f}, y_min={ymf:.2f} (fast_paste)"
)
print(
"stitch_feature_optimized steps: "
f"brightness={t_after_brightness - t_profile0:.6f}s, "
f"geom={t_after_geom - t_after_brightness:.6f}s ({'hit' if geom_cache_hit else 'miss'}), "
f"reset={t_after_reset - t_after_geom:.6f}s, "
f"paste={t_after_paste - t_after_reset:.6f}s ({pasted_tiles + 1} tiles incl ref), "
f"finalize={t_after_finalize - t_after_paste:.6f}s, "
f"total={t_after_finalize - t_profile0:.6f}s"
)
return stitch_buf
def _render_feature_tiles(
self,
t_ref_bgr,
@@ -579,77 +933,37 @@ class Stitcher:
tiles_bgr_c = tiles_bgr
t_after_brightness = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
geom_key = self._geom_cache_key(t_ref_bgr, tiles_bgr_c, pad)
geom_cache_hit = canvas_cache.get("geom_key") == geom_key
if not geom_cache_hit:
ref_corners = np.array([[0, 0], [w0c, 0], [w0c, h0c], [0, h0c]], dtype=np.float64)
pts = [ref_corners]
for tile_bgr, H_use in zip(tiles_bgr_c, H_use_list):
if H_use is None:
continue
hi, wi = tile_bgr.shape[:2]
tx, ty = self._translation_from_H_to_ref(H_use)
pts.append(
np.array(
[[tx, ty], [tx + wi, ty], [tx + wi, ty + hi], [tx, ty + hi]],
dtype=np.float64,
)
)
all_xy = np.concatenate(pts, axis=0)
x_min_f = float(all_xy[:, 0].min()) - pad
y_min_f = float(all_xy[:, 1].min()) - pad
x_max_f = float(all_xy[:, 0].max()) + pad
y_max_f = float(all_xy[:, 1].max()) + pad
out_w = max(1, int(np.ceil(x_max_f - x_min_f)))
out_h = max(1, int(np.ceil(y_max_f - y_min_f)))
ref_y0 = int(np.floor(-y_min_f))
ref_x0 = int(np.floor(-x_min_f))
pastes = []
for H_use in H_use_list:
if H_use is None:
pastes.append(None)
continue
tx, ty = self._translation_from_H_to_ref(H_use)
pastes.append((int(np.floor(ty - y_min_f)), int(np.floor(tx - x_min_f))))
_dr_y0 = max(0, ref_y0)
_dr_x0 = max(0, ref_x0)
_dr_y1 = min(out_h, ref_y0 + h0c)
_dr_x1 = min(out_w, ref_x0 + w0c)
for _tb, _pa in zip(tiles_bgr_c, pastes):
if _pa is None:
continue
_py, _px = _pa
_th, _tw = _tb.shape[:2]
_dr_y0 = min(_dr_y0, max(0, _py))
_dr_x0 = min(_dr_x0, max(0, _px))
_dr_y1 = max(_dr_y1, min(out_h, _py + _th))
_dr_x1 = max(_dr_x1, min(out_w, _px + _tw))
canvas_cache["geom_key"] = geom_key
canvas_cache["out_w"] = out_w
canvas_cache["out_h"] = out_h
canvas_cache["ref_y0"] = ref_y0
canvas_cache["ref_x0"] = ref_x0
canvas_cache["pastes"] = pastes
canvas_cache["x_min_f"] = x_min_f
canvas_cache["y_min_f"] = y_min_f
canvas_cache["crop_rect"] = None
canvas_cache["dirty_rect"] = (_dr_y0, _dr_y1, _dr_x0, _dr_x1)
canvas_cache["blend_buf"] = np.zeros((out_h, out_w, 3), dtype=np.float32)
canvas_cache["wsum_buf"] = np.zeros((out_h, out_w), dtype=np.float32)
canvas_cache["cov_buf"] = np.zeros((out_h, out_w), dtype=np.uint8)
canvas_cache["weight_maps"] = {}
else:
out_w = canvas_cache["out_w"]
out_h = canvas_cache["out_h"]
ref_y0 = canvas_cache["ref_y0"]
ref_x0 = canvas_cache["ref_x0"]
pastes = canvas_cache["pastes"]
geom_cache_hit = self._ensure_feature_geom_cache(
canvas_cache,
t_ref_bgr,
tiles_bgr_c,
H_use_list,
pad,
allocate_blend_bufs=bool(blend),
)
t_after_geom = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
if not blend:
return self._paste_feature_tiles_fast(
t_ref_bgr,
tiles_bgr_c,
H_use_list,
canvas_cache,
geom_cache_hit=geom_cache_hit,
frame_idx=frame_idx,
print_every=print_every,
verbose=verbose,
t_profile0=t_profile0,
t_after_brightness=t_after_brightness,
t_after_geom=t_after_geom,
)
out_w = canvas_cache["out_w"]
out_h = canvas_cache["out_h"]
ref_y0 = canvas_cache["ref_y0"]
ref_x0 = canvas_cache["ref_x0"]
pastes = canvas_cache["pastes"]
blend_buf = canvas_cache["blend_buf"]
wsum = canvas_cache["wsum_buf"]
cov = canvas_cache["cov_buf"]
@@ -667,7 +981,7 @@ class Stitcher:
cov.fill(0)
t_after_reset = cv2.getTickCount() / cv2.getTickFrequency() if do_profile else 0.0
blend_margin = 16 if blend else 0
blend_margin = 16
blended_tiles = 0
for tile_bgr, H_use, paste in zip(tiles_bgr_c, H_use_list, pastes):
+1 -1
View File
@@ -16,7 +16,7 @@ class CompoundEyeFrameProcessor:
def process(self, frame):
self.stitcher.get_stitch_unit(frame)
self.stitcher.extract_tiles()
return self.stitcher.stitch_rt_optimized(blend=False)
return self.stitcher.stitch_rt_optimized()
class SingleCameraFrameProcessor:
+10
View File
@@ -0,0 +1,10 @@
"""Pure-Python Numba stitch kernels kept outside the Nuitka-compiled orisys .pyd.
Nuitka cannot reliably JIT ``@njit`` functions baked into a compiled extension.
This package is shipped as source next to the Nuitka wheel so Numba can see
real CPython bytecode at runtime.
"""
from .kernel import blend_kernel
__all__ = ["blend_kernel"]
+51
View File
@@ -0,0 +1,51 @@
"""Fused stitch blend kernel (Numba JIT when available)."""
from __future__ import annotations
import numpy as np
try:
from numba import njit
_NUMBA_AVAILABLE = True
except ImportError:
_NUMBA_AVAILABLE = False
def njit(*args, **kwargs):
def decorator(fn):
return fn
return decorator
if _NUMBA_AVAILABLE:
@njit(cache=True)
def blend_kernel(blend_patch, wsum_patch, cov_patch, src_patch, w_patch):
"""Fused incremental weighted-average composite (numba-accelerated)."""
sh = src_patch.shape[0]
sw = src_patch.shape[1]
for y in range(sh):
for x in range(sw):
w = w_patch[y, x]
s0 = np.float32(src_patch[y, x, 0])
s1 = np.float32(src_patch[y, x, 1])
s2 = np.float32(src_patch[y, x, 2])
if cov_patch[y, x] == 0:
blend_patch[y, x, 0] = s0
blend_patch[y, x, 1] = s1
blend_patch[y, x, 2] = s2
wsum_patch[y, x] = w
else:
ws = wsum_patch[y, x] + w
if ws < 1e-6:
ws = 1e-6
a = w / ws
blend_patch[y, x, 0] += (s0 - blend_patch[y, x, 0]) * a
blend_patch[y, x, 1] += (s1 - blend_patch[y, x, 1]) * a
blend_patch[y, x, 2] += (s2 - blend_patch[y, x, 2]) * a
wsum_patch[y, x] = ws
cov_patch[y, x] = 1
else:
blend_kernel = None
-121
View File
@@ -1,121 +0,0 @@
# 3D Patch Generation Blender 工具
本文件夹包含 Blender 场景文件以及辅助脚本,用于导出网格体/相机数据,并将一张采集到的 2D 图像重新映射回 3D 表面。
## 主要文件
- `dd02c_v2.blend`
DD02C 曲面与相机配置所使用的主要 Blender 工程文件。
- `blender_export.py`
在 Blender 中运行,用于导出以下数据:
- 网格几何信息
- 与 UV 相关的 loop 数据
- 相机内参
- 相机位姿
- 预计算的像素到 3D 查找表
- `mapto3d.py`
读取导出的数据和一张采集图像,并生成带颜色的 PLY 格式 3D 点云。
## 导出 / 生成的数据
- `exported_data/`
由 `blender_export.py` 生成。通常应包含以下文件:
- `mesh.npz`
- `camera_intrinsics.npz`
- `world_to_cam.npy`
- `pixel_to_xyz.npz`
- `captured_image.png`
要映射回导出 3D 表面的输入图像。当前该图像由 Blender 渲染得到。
- `mapped_points.ply`
由 `mapto3d.py` 生成的输出点云文件。
- `texture/`
保存实验过程中生成的纹理图片。
## 工作流程
### 1. 准备 Blender 场景
在 Blender 中打开 `dd02c_v2.blend`。
请确认:
- 已选中目标网格
- 相机存在且位置正确
- 渲染分辨率与采集图像分辨率一致
### 2. 导出几何和相机数据
在 Blender 的里面上面转到“scripts” 模式,然后点击运行 `blender_export.py`。
该脚本会将所需数据导出到:
```text
exported_data/
```
脚本会写出以下内容:
- 网格顶点和面信息
- UV loop 数据
- OpenCV 风格的相机内参
- world-to-camera 变换矩阵
- 可选的逐像素 3D 查找表
### 3. 准备输入图像
将需要投影到 3D 表面的图像放到当前文件夹中,可以在Blender 渲染生成,例如:
```text
captured_image.png
```
`mapto3d.py` 目前在其 `__main__` 代码块中直接指定输入图像路径。
### 4. 将图像像素映射到 3D
运行:
```bash
python mapto3d.py
```
该脚本会:
- 加载 Blender 导出的数据
- 加载输入图像
- 在需要时对图像做上下翻转,以匹配导出的像素映射顺序
- 优先使用预计算的 `pixel_to_xyz.npz`
- 将结果保存为 `mapped_points.ply`
如果缺少预计算查找表,脚本会退回到逐像素射线求交,运行速度会慢很多。
### 5. 查看输出结果
生成的点云可以使用支持 PLY 的工具打开,例如:
- Open3D
- MeshLab
- CloudCompare
`mapto3d.py` 中也包含可选的 Open3D 可视化步骤。
## 说明
- `mapto3d.py` 目前假设导出的相机模型与输入图像分辨率一致。
- 脚本现在默认会对加载的图像做上下翻转,以保证映射颜色与 Blender 导出的像素顺序一致。
- 为了获得更好的性能,建议在 `exported_data/` 中保留 `pixel_to_xyz.npz`。
## 依赖
### Blender 侧
- Blender Python API:`bpy`
- `numpy`
### Python 侧
- `numpy`
- `opencv-python`
- `open3d`(可选,用于可视化)
@@ -1,147 +0,0 @@
import bpy
import numpy as np
import json
import os
from mathutils import Matrix, Vector
def export_camera_intrinsics(camera_obj, scene):
"""Return camera intrinsic matrix K (3x3) and distortion (none)."""
cam = camera_obj.data
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
# Sensor size and focal length
if cam.sensor_fit == 'HORIZONTAL':
sensor_width = cam.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam.lens # in mm
# Pixel sizes in mm
pixel_width = sensor_width / res_x
pixel_height = sensor_height / res_y
# Intrinsic matrix (OpenCV style: fx, fy in pixels, cx, cy principal point)
fx = focal / pixel_width
fy = focal / pixel_height
cx = res_x / 2.0
cy = res_y / 2.0
K = np.array([[fx, 0, cx],
[0, fy, cy],
[0, 0, 1]], dtype=np.float64)
return K, res_x, res_y
def export_mesh_data(obj, filepath):
"""Export mesh vertices, faces, and UVs to a compressed .npz file."""
me = obj.data
# Transform vertices to world space
vertices = np.array([obj.matrix_world @ v.co for v in me.vertices], dtype=np.float32)
faces = np.array([poly.vertices[:] for poly in me.polygons], dtype=np.int32)
# UVs (per loop, need mapping to vertices for simplicity we'll store per-face)
uv_layer = me.uv_layers.active
if uv_layer:
uv_data = np.array([uv_layer.data[i].uv for i in range(len(uv_layer.data))], dtype=np.float32)
loop_to_vertex = np.array([me.loops[i].vertex_index for i in range(len(me.loops))], dtype=np.int32)
else:
uv_data = np.zeros((0,2), dtype=np.float32)
loop_to_vertex = np.zeros(0, dtype=np.int32)
np.savez_compressed(filepath,
vertices=vertices,
faces=faces,
uv_data=uv_data,
loop_to_vertex=loop_to_vertex)
print(f"Mesh exported to {filepath}")
def export_camera_pose(camera_obj, filepath):
"""Export camera extrinsics (world-to-camera transform matrix)."""
# Camera world matrix: transforms from camera local to world
cam_to_world = np.array(camera_obj.matrix_world, dtype=np.float32)
# We need world-to-camera transform (inverse)
world_to_cam = np.linalg.inv(cam_to_world)
np.save(filepath, world_to_cam)
print(f"Camera pose exported to {filepath}")
def export_pixel_to_3d_lookup(obj, camera_obj, scene, filepath):
"""
Precompute for every pixel the 3D hit point (world coordinates)
and save as a 2D array (height, width, 3). This avoids raycasting later.
"""
depsgraph = bpy.context.evaluated_depsgraph_get()
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
# Get camera data
cam = camera_obj.data
if cam.sensor_fit == 'HORIZONTAL':
sensor_width = cam.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
# Preallocate array (height, width, 3) for XYZ world coordinates
xyz_map = np.full((res_y, res_x, 3), np.nan, dtype=np.float32)
origin = camera_obj.matrix_world.translation
cam_rot = camera_obj.matrix_world.to_quaternion()
for y in range(res_y):
# Normalized sensor coordinate v (0=bottom, 1=top)
v_ndc = (y + 0.5) / res_y
sy = (v_ndc - 0.5) * frame_height
for x in range(res_x):
u_ndc = (x + 0.5) / res_x
sx = (u_ndc - 0.5) * frame_width
direction_cam = Vector((sx, sy, -1.0))
direction_world = cam_rot @ direction_cam
direction_world.normalize()
hit, location, normal, face_index, hit_obj, matrix = scene.ray_cast(
depsgraph, origin, direction_world)
if hit and hit_obj == obj:
xyz_map[y, x, :] = location
np.savez_compressed(filepath, xyz_map=xyz_map)
print(f"Pixel-to-3D lookup saved to {filepath} (shape {res_y}x{res_x})")
def export_all():
# Configuration
output_dir = bpy.path.abspath("//exported_data/")
os.makedirs(output_dir, exist_ok=True)
obj = bpy.context.active_object # or specify name
cam = bpy.data.objects.get("Camera")
if not cam:
cam = bpy.context.scene.camera
if not obj or not cam:
raise RuntimeError("Select mesh and ensure camera exists.")
# 1. Mesh data
export_mesh_data(obj, os.path.join(output_dir, "mesh.npz"))
# 2. Camera intrinsics (K matrix) and resolution
K, w, h = export_camera_intrinsics(cam, bpy.context.scene)
np.savez(os.path.join(output_dir, "camera_intrinsics.npz"), K=K, width=w, height=h)
# 3. Camera pose (world-to-camera)
export_camera_pose(cam, os.path.join(output_dir, "world_to_cam.npy"))
# 4. Optional: precomputed pixel-to-3D map (fastest)
export_pixel_to_3d_lookup(obj, cam, bpy.context.scene,
os.path.join(output_dir, "pixel_to_xyz.npz"))
print("Export complete.")
if __name__ == "__main__":
export_all()
Binary file not shown.

Before

Width:  |  Height:  |  Size: 1.1 MiB

@@ -1,585 +0,0 @@
"""
Texture Generator for Arbitrary 3D Surfaces (Blender Script)
- Mode: "patches" (black/white checker patches) or "dots" (circular markers)
- Mesh UVs are filled white so that non‑visible parts become white.
- Texture background (no‑mesh areas) stays black.
"""
import bpy
import numpy as np
import os
import time
from mathutils import Vector
# ----------------------------------------------------------------------
# Common Utilities
# ----------------------------------------------------------------------
def draw_cross(pixels, width, height, cx_uv, cy_uv, length_px, thickness_px, color):
"""
Draw a cross (plus sign) in pixel space.
cx_uv, cy_uv: centre of cross in UV coordinates (0‑1), origin at bottom‑left.
length_px: total length of each arm (pixels).
thickness_px: line thickness (odd number gives symmetric result).
"""
# Convert UV centre to pixel coordinates (y from bottom)
cx = int(round(cx_uv * (width - 1)))
cy = int(round(cy_uv * (height - 1)))
half_len = length_px // 2
half_thick = thickness_px // 2
# Horizontal arm
x0 = max(0, cx - half_len)
x1 = min(width - 1, cx + half_len)
y0 = max(0, cy - half_thick)
y1 = min(height - 1, cy + half_thick)
for y in range(y0, y1 + 1):
for x in range(x0, x1 + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
# Vertical arm
x0 = max(0, cx - half_thick)
x1 = min(width - 1, cx + half_thick)
y0 = max(0, cy - half_len)
y1 = min(height - 1, cy + half_len)
for y in range(y0, y1 + 1):
for x in range(x0, x1 + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def fill_background(pixels, width, height, color):
"""Fill the entire image with a solid RGBA color (0-1 floats)."""
r, g, b, a = color
for y in range(height):
for x in range(width):
idx = (y * width + x) * 4
pixels[idx] = r
pixels[idx+1] = g
pixels[idx+2] = b
pixels[idx+3] = a
def fill_polygon(image, pixels, width, height, polygon, color):
"""Fill a polygon in the pixel array (flattened, RGBA)."""
poly = [(int(round(x)), int(round(y))) for x, y in polygon]
if len(poly) < 3:
return
min_x = max(0, min(p[0] for p in poly))
max_x = min(width-1, max(p[0] for p in poly))
min_y = max(0, min(p[1] for p in poly))
max_y = min(height-1, max(p[1] for p in poly))
for y in range(min_y, max_y + 1):
intersections = []
for i in range(len(poly)):
x1, y1 = poly[i]
x2, y2 = poly[(i+1) % len(poly)]
if y1 == y2:
continue
if (y1 <= y < y2) or (y2 <= y < y1):
t = (y - y1) / (y2 - y1)
x = x1 + t * (x2 - x1)
intersections.append(x)
intersections.sort()
for k in range(0, len(intersections), 2):
if k+1 >= len(intersections):
break
x_start = int(round(intersections[k]))
x_end = int(round(intersections[k+1]))
x_start = max(0, min(width-1, x_start))
x_end = max(0, min(width-1, x_end))
for x in range(x_start, x_end + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def fill_polygon_value(grid, width, height, polygon, value):
"""Fill a polygon into a 2D scalar grid."""
poly = [(int(round(x)), int(round(y))) for x, y in polygon]
if len(poly) < 3:
return
min_x = max(0, min(p[0] for p in poly))
max_x = min(width - 1, max(p[0] for p in poly))
min_y = max(0, min(p[1] for p in poly))
max_y = min(height - 1, max(p[1] for p in poly))
for y in range(min_y, max_y + 1):
intersections = []
for i in range(len(poly)):
x1, y1 = poly[i]
x2, y2 = poly[(i + 1) % len(poly)]
if y1 == y2:
continue
if (y1 <= y < y2) or (y2 <= y < y1):
t = (y - y1) / (y2 - y1)
x = x1 + t * (x2 - x1)
intersections.append(x)
intersections.sort()
for k in range(0, len(intersections), 2):
if k + 1 >= len(intersections):
break
x_start = int(round(intersections[k]))
x_end = int(round(intersections[k + 1]))
x_start = max(0, min(width - 1, x_start))
x_end = max(0, min(width - 1, x_end))
grid[y, x_start:x_end + 1] = value
def draw_ellipse(pixels, width, height, cx, cy, a, b, angle, color):
"""Filled ellipse with semi‑axes a,b (pixels), rotation angle (rad)."""
cos_a = np.cos(angle)
sin_a = np.sin(angle)
rx = abs(a * cos_a) + abs(b * sin_a)
ry = abs(a * sin_a) + abs(b * cos_a)
x_min = int(max(0, cx - rx))
x_max = int(min(width-1, cx + rx))
y_min = int(max(0, cy - ry))
y_max = int(min(height-1, cy + ry))
for y in range(y_min, y_max + 1):
for x in range(x_min, x_max + 1):
dx = x - cx
dy = y - cy
u = dx * cos_a + dy * sin_a
v = -dx * sin_a + dy * cos_a
if (u/a)**2 + (v/b)**2 <= 1.0:
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def get_uv_from_face(obj, face_index, location):
"""World‑space location → UV coordinate on the mesh face."""
me = obj.data
local_loc = obj.matrix_world.inverted() @ location
poly = me.polygons[face_index]
vert_indices = poly.vertices
loop_indices = poly.loop_indices
if len(vert_indices) < 3:
return None
v0 = me.vertices[vert_indices[0]].co
v1 = me.vertices[vert_indices[1]].co
v2 = me.vertices[vert_indices[2]].co
uv_layer = me.uv_layers.active.data
uv0 = uv_layer[loop_indices[0]].uv
uv1 = uv_layer[loop_indices[1]].uv
uv2 = uv_layer[loop_indices[2]].uv
def tri_area(a, b, c):
return 0.5 * ((b - a).cross(c - a)).length
total_area = tri_area(v0, v1, v2)
if total_area == 0:
return None
w0 = tri_area(local_loc, v1, v2) / total_area
w1 = tri_area(local_loc, v2, v0) / total_area
w2 = 1.0 - w0 - w1
uv = w0 * uv0 + w1 * uv1 + w2 * uv2
return Vector((uv.x, uv.y))
def raycast_from_camera(scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height):
"""Camera raycast through normalised sensor coordinates (0‑1)."""
sx = (u_ndc - 0.5) * frame_width
sy = (v_ndc - 0.5) * frame_height
direction_cam = Vector((sx, sy, -1.0))
direction_world = cam.matrix_world.to_quaternion() @ direction_cam
direction_world.normalize()
origin = cam.matrix_world.translation
hit, location, normal, face_index, hit_obj, matrix = scene.ray_cast(
depsgraph, origin, direction_world)
if hit and hit_obj == obj:
return location, face_index
return None, None
def rasterize_mesh_uvs(obj, pixels, tex_width, tex_height, color=(1,1,1,1)):
"""
Fill all triangles of the active UV map with 'color' (default white).
This marks every part of the mesh, visible or not, as white.
"""
me = obj.data
uv_layer = me.uv_layers.active
if not uv_layer:
print("Warning: No active UV layer – mesh rasterization skipped.")
return
uv_data = uv_layer.data
for poly in me.polygons:
# Get loop indices for this polygon
loops = poly.loop_indices
# Build UV polygon
uv_poly = []
for loop_idx in loops:
uv = uv_data[loop_idx].uv
uv_poly.append((uv.x, uv.y))
# Convert UV (0‑1) to pixel coordinates
pixel_poly = [(u * (tex_width - 1), v * (tex_height - 1)) for u, v in uv_poly]
fill_polygon(None, pixels, tex_width, tex_height, pixel_poly, color)
def rasterize_mesh_mask(obj, tex_width, tex_height):
"""Rasterize the active UV map to a boolean mesh mask."""
mask = np.zeros((tex_height, tex_width), dtype=bool)
me = obj.data
uv_layer = me.uv_layers.active
if not uv_layer:
return mask
uv_data = uv_layer.data
for poly in me.polygons:
uv_poly = []
for loop_idx in poly.loop_indices:
uv = uv_data[loop_idx].uv
uv_poly.append((uv.x, uv.y))
pixel_poly = [(u * (tex_width - 1), v * (tex_height - 1)) for u, v in uv_poly]
fill_polygon_value(mask, tex_width, tex_height, pixel_poly, True)
return mask
def paint_mask_color(pixels, width, height, mask, color):
"""Paint all True pixels of a boolean mask with an RGBA color."""
r, g, b, a = color
ys, xs = np.nonzero(mask)
for y, x in zip(ys, xs):
idx = (y * width + x) * 4
pixels[idx:idx+4] = (r, g, b, a)
def complete_checker_gaps(pixels, width, height, mesh_mask, patch_state, max_gap_px=6):
"""
Fill small white slivers inside the mesh by expanding nearby checker colors.
Only grows locally (max_gap_px) so non-visible mesh areas remain white.
patch_state: -1 = unassigned, 0 = black patch, 1 = white patch
"""
for _ in range(max_gap_px):
candidate = mesh_mask & (patch_state < 0)
if not np.any(candidate):
break
black_counts = np.zeros_like(patch_state, dtype=np.uint8)
white_counts = np.zeros_like(patch_state, dtype=np.uint8)
for dy in (-1, 0, 1):
for dx in (-1, 0, 1):
if dx == 0 and dy == 0:
continue
src_y0 = max(0, -dy)
src_y1 = height - max(0, dy)
src_x0 = max(0, -dx)
src_x1 = width - max(0, dx)
dst_y0 = max(0, dy)
dst_y1 = height - max(0, -dy)
dst_x0 = max(0, dx)
dst_x1 = width - max(0, -dx)
neighbor_vals = patch_state[src_y0:src_y1, src_x0:src_x1]
black_counts[dst_y0:dst_y1, dst_x0:dst_x1] += (neighbor_vals == 0)
white_counts[dst_y0:dst_y1, dst_x0:dst_x1] += (neighbor_vals == 1)
fillable = candidate & ((black_counts + white_counts) > 0)
if not np.any(fillable):
break
patch_state[fillable] = np.where(
black_counts[fillable] > white_counts[fillable], 0, 1
).astype(np.int8)
paint_mask_color(pixels, width, height, patch_state == 0, (0.0, 0.0, 0.0, 1.0))
paint_mask_color(pixels, width, height, patch_state == 1, (1.0, 1.0, 1.0, 1.0))
def checker_patch_color(i, j):
"""Return a black/white checker color for patch index (i, j)."""
if (i + j) % 2 == 0:
return np.array([1.0, 1.0, 1.0], dtype=np.float32)
return np.array([0.0, 0.0, 0.0], dtype=np.float32)
# ----------------------------------------------------------------------
# Mode 1: Colored Patches
# ----------------------------------------------------------------------
def generate_patch_texture(
obj, cam, scene, depsgraph,
nx=40, ny=30, r=0.05,
texture_width=2048, texture_height=2048,
output_path="//patch_texture.png",uv_centre=None, cross_length=40, cross_thickness=3, cross_color=(1,0,0,1),
gap_fill_px=6
):
cam_data = cam.data
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
if cam_data.sensor_fit == 'HORIZONTAL':
sensor_width = cam_data.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam_data.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam_data.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
patch_uv_polygons = []
patch_colors = np.zeros((ny, nx, 3), dtype=np.float32)
print(f"Raycasting {nx*ny} patches...")
for j in range(ny):
v_start = j / ny
v_end = (j+1) / ny
for i in range(nx):
u_start = i / nx
u_end = (i+1) / nx
corners_ndc = [(u_start, v_start), (u_end, v_start),
(u_end, v_end), (u_start, v_end)]
uv_poly = []
for (u_ndc, v_ndc) in corners_ndc:
loc, face_idx = raycast_from_camera(scene, depsgraph, cam, obj,
u_ndc, v_ndc, frame_width, frame_height)
if loc is not None:
uv = get_uv_from_face(obj, face_idx, loc)
if uv is not None:
uv_poly.append(uv)
if len(uv_poly) < 3:
patch_uv_polygons.append(None)
continue
patch_uv_polygons.append(uv_poly)
patch_colors[j, i] = checker_patch_color(i, j)
# Create image – start with black background
img_name = "PatchTexture"
if img_name in bpy.data.images:
bpy.data.images.remove(bpy.data.images[img_name])
image = bpy.data.images.new(img_name, width=texture_width, height=texture_height, alpha=True)
pixels = list(image.pixels)
mesh_mask = rasterize_mesh_mask(obj, texture_width, texture_height)
patch_state = np.full((texture_height, texture_width), -1, dtype=np.int8)
# Fill entire image with black (no mesh areas)
fill_background(pixels, texture_width, texture_height, (0,0,0,1))
# Rasterize all mesh UVs in white (visible + non‑visible parts)
paint_mask_color(pixels, texture_width, texture_height, mesh_mask, (1,1,1,1))
# Now draw the visible patches on top
idx = 0
for j in range(ny):
for i in range(nx):
uv_poly = patch_uv_polygons[idx]; idx += 1
if uv_poly is None: continue
pixel_poly = [(u*(texture_width-1), v*(texture_height-1)) for u,v in uv_poly]
color = patch_colors[j,i]
rgba = (color[0], color[1], color[2], 1.0)
fill_polygon(image, pixels, texture_width, texture_height, pixel_poly, rgba)
fill_polygon_value(
patch_state,
texture_width,
texture_height,
pixel_poly,
1 if np.mean(color) >= 0.5 else 0,
)
# Fill small steep-edge gaps from neighboring checker cells.
complete_checker_gaps(
pixels, texture_width, texture_height, mesh_mask, patch_state, max_gap_px=gap_fill_px
)
# (after drawing all patches)
if uv_centre is not None:
draw_cross(pixels, texture_width, texture_height,
uv_centre[0], uv_centre[1],
cross_length, cross_thickness, cross_color)
image.pixels[:] = pixels
image.update()
abs_path = bpy.path.abspath(output_path)
image.filepath_raw = abs_path
image.file_format = 'PNG'
image.save()
print(f"Patch texture saved to {abs_path}")
return image
# ----------------------------------------------------------------------
# Mode 2: Circular Dot Markers
# ----------------------------------------------------------------------
def generate_dot_texture(
obj, cam, scene, depsgraph,
dot_diameter_pixels=12,
spacing_pixels=40,
dot_color=(0,0,0,1), # RGBA 0-1
texture_width=2048, texture_height=2048,
output_path="//dot_texture.png",uv_centre=None, cross_length=40, cross_thickness=3, cross_color=(1,0,0,1)
):
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
cam_data = cam.data
if cam_data.sensor_fit == 'HORIZONTAL':
sensor_width = cam_data.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam_data.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam_data.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
margin = spacing_pixels // 2
x_centers = np.arange(margin, res_x, spacing_pixels)
y_centers = np.arange(margin, res_y, spacing_pixels)
offset_pix = max(2, dot_diameter_pixels // 4)
du_ndc = offset_pix / res_x
dv_ndc = offset_pix / res_y
dots = [] # (cx_uv, cy_uv, a_uv, b_uv, angle_uv)
print(f"Raycasting {len(x_centers)*len(y_centers)} dot centers...")
for y_pix in y_centers:
v_ndc = y_pix / res_y
for x_pix in x_centers:
u_ndc = x_pix / res_x
center_loc, face_idx = raycast_from_camera(scene, depsgraph, cam, obj,
u_ndc, v_ndc, frame_width, frame_height)
if center_loc is None: continue
center_uv = get_uv_from_face(obj, face_idx, center_loc)
if center_uv is None: continue
offsets = [
( u_ndc + du_ndc, v_ndc),
( u_ndc - du_ndc, v_ndc),
( u_ndc, v_ndc + dv_ndc),
( u_ndc, v_ndc - dv_ndc)
]
uv_offsets = []
valid = True
for (u_off, v_off) in offsets:
loc_off, f_idx_off = raycast_from_camera(scene, depsgraph, cam, obj,
u_off, v_off, frame_width, frame_height)
if loc_off is None:
valid = False; break
uv_off = get_uv_from_face(obj, f_idx_off, loc_off)
if uv_off is None:
valid = False; break
uv_offsets.append(uv_off)
if not valid: continue
scale_u = 1.0 / res_x
scale_v = 1.0 / res_y
du_dx = (uv_offsets[0] - uv_offsets[1]) / (2 * offset_pix * scale_u)
dv_dx = (uv_offsets[2] - uv_offsets[3]) / (2 * offset_pix * scale_v)
r_pix = dot_diameter_pixels / 2.0
a_uv = r_pix * du_dx.length * scale_u
b_uv = r_pix * dv_dx.length * scale_v
angle = np.arctan2(du_dx.y, du_dx.x) if du_dx.length > 0 else 0.0
dots.append((center_uv.x, center_uv.y, a_uv, b_uv, angle))
print(f"Found {len(dots)} valid dots.")
if uv_centre is not None:
draw_cross(pixels, texture_width, texture_height,
uv_centre[0], uv_centre[1],
cross_length, cross_thickness, cross_color)
# Create image
img_name = "DotTexture"
if img_name in bpy.data.images:
bpy.data.images.remove(bpy.data.images[img_name])
image = bpy.data.images.new(img_name, width=texture_width, height=texture_height, alpha=True)
pixels = list(image.pixels)
# Black background (no mesh)
fill_background(pixels, texture_width, texture_height, (0,0,0,1))
# White for all mesh areas
rasterize_mesh_uvs(obj, pixels, texture_width, texture_height, color=(1,1,1,1))
# Draw dots
for (cx_uv, cy_uv, a_uv, b_uv, angle_uv) in dots:
cx_px = cx_uv * (texture_width - 1)
cy_px = cy_uv * (texture_height - 1)
a_px = a_uv * texture_width
b_px = b_uv * texture_height
draw_ellipse(pixels, texture_width, texture_height,
cx_px, cy_px, a_px, b_px, angle_uv, dot_color)
image.pixels[:] = pixels
image.update()
abs_path = bpy.path.abspath(output_path)
image.filepath_raw = abs_path
image.file_format = 'PNG'
image.save()
print(f"Dot texture saved to {abs_path}")
return image
# ----------------------------------------------------------------------
# Main
# ----------------------------------------------------------------------
def main():
# ==================== CONFIGURATION ====================
MODE = "patches"
OBJ_NAME = "FingerMesh"
CAM_NAME = "Camera"
TEXTURE_WIDTH = 4096
TEXTURE_HEIGHT = 4096
OUTPUT_PATH = "//" + ("patch_texture.png" if MODE=="patches" else "dot_texture.png")
# Patch mode
PATCH_NX = 40; PATCH_NY = 40; COLOR_R = 0.05
# Dot mode
DOT_DIAMETER_PIX = 12; DOT_SPACING_PIX = 40; DOT_COLOR = (0,0,0,1)
# Cross
CROSS_LENGTH_PX = 40
CROSS_THICKNESS_PX = 3
CROSS_COLOR = (1, 0, 0, 1) # red cross
# ======================================================
start_time = time.time()
scene = bpy.context.scene
depsgraph = bpy.context.evaluated_depsgraph_get()
if OBJ_NAME:
obj = bpy.data.objects.get(OBJ_NAME)
else:
obj = bpy.context.active_object
if obj is None or obj.type != 'MESH':
raise ValueError("No valid mesh object found.")
cam = bpy.data.objects.get(CAM_NAME) if CAM_NAME else scene.camera
if cam is None or cam.type != 'CAMERA':
raise ValueError("No camera found.")
if not obj.data.uv_layers.active:
obj.data.uv_layers.new(name="UVMap")
# Compute UV centre
uv_layer = obj.data.uv_layers.active
uv_coords = np.array([uv_layer.data[i].uv.to_tuple() for i in range(len(uv_layer.data))])
uv_min = uv_coords.min(axis=0)
uv_max = uv_coords.max(axis=0)
uv_centre = (uv_min + uv_max) / 2.0
print(f"UV centre: ({uv_centre[0]:.4f}, {uv_centre[1]:.4f})")
# Compute UV bounding box centre
me = obj.data
uv_layer = me.uv_layers.active
uv_coords = np.array([uv_layer.data[i].uv.to_tuple() for i in range(len(uv_layer.data))])
uv_min = uv_coords.min(axis=0)
uv_max = uv_coords.max(axis=0)
uv_centre = (uv_min + uv_max) / 2.0
print(f"UV centre: {uv_centre[0]:.4f}, {uv_centre[1]:.4f}")
if MODE == "patches":
generate_patch_texture(
obj, cam, scene, depsgraph,
nx=PATCH_NX, ny=PATCH_NY, r=COLOR_R,
texture_width=TEXTURE_WIDTH, texture_height=TEXTURE_HEIGHT,
output_path=OUTPUT_PATH,
uv_centre=uv_centre,
cross_length=CROSS_LENGTH_PX,
cross_thickness=CROSS_THICKNESS_PX,
cross_color=CROSS_COLOR
)
elif MODE == "dots":
generate_dot_texture(
obj, cam, scene, depsgraph,
dot_diameter_pixels=DOT_DIAMETER_PIX,
spacing_pixels=DOT_SPACING_PIX,
dot_color=DOT_COLOR,
texture_width=TEXTURE_WIDTH, texture_height=TEXTURE_HEIGHT,
output_path=OUTPUT_PATH,
uv_centre=uv_centre,
cross_length=CROSS_LENGTH_PX,
cross_thickness=CROSS_THICKNESS_PX,
cross_color=CROSS_COLOR
)
print(f"Completed in {time.time()-start_time:.2f} seconds.")
if __name__ == "__main__":
main()
Binary file not shown.
Binary file not shown.
@@ -1,411 +0,0 @@
import argparse
import os
import time
import cv2
import numpy as np
def load_exported_data(data_dir):
"""Load Blender-exported mesh/camera files."""
mesh = np.load(os.path.join(data_dir, "mesh.npz"))
cam_data = np.load(os.path.join(data_dir, "camera_intrinsics.npz"))
xyz_map_path = os.path.join(data_dir, "pixel_to_xyz.npz")
if not os.path.exists(xyz_map_path):
raise FileNotFoundError(
f"Real-time video viewing requires precomputed pixel_to_xyz.npz: {xyz_map_path}"
)
xyz_map = np.load(xyz_map_path)["xyz_map"]
width = int(cam_data["width"])
height = int(cam_data["height"])
if xyz_map.shape[:2] != (height, width):
raise ValueError(
f"pixel_to_xyz shape {xyz_map.shape[:2]} does not match camera size {(height, width)}"
)
return {
"vertices": mesh["vertices"],
"faces": mesh["faces"],
"uv_data": mesh["uv_data"],
"loop_to_vertex": mesh["loop_to_vertex"],
"K": cam_data["K"],
"width": width,
"height": height,
"world_to_cam": np.load(os.path.join(data_dir, "world_to_cam.npy")),
"xyz_map": xyz_map,
}
def prepare_static_points(data, display_stride=4):
"""
Prepare fixed 3D points from the precomputed pixel-to-3D map.
The geometry does not change between video frames, so real-time viewing only
needs to update point colors.
"""
xyz_map = data["xyz_map"]
valid_mask = ~np.isnan(xyz_map[:, :, 0])
valid_points = xyz_map[valid_mask].reshape(-1, 3)
if len(valid_points) == 0:
raise ValueError("pixel_to_xyz.npz contains no valid mesh hit points.")
display_stride = max(1, int(display_stride))
display_indices = np.arange(len(valid_points))[::display_stride]
return valid_mask, display_indices, valid_points[display_indices]
def frame_to_rgb_colors(frame_bgr, data, valid_mask, display_indices, flip_vertical=True):
"""Convert a video frame to RGB point colors aligned with the exported map."""
target_w = int(data["width"])
target_h = int(data["height"])
h, w = frame_bgr.shape[:2]
if (w, h) != (target_w, target_h):
frame_bgr = cv2.resize(frame_bgr, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
frame_rgb = cv2.cvtColor(frame_bgr, cv2.COLOR_BGR2RGB)
if flip_vertical:
frame_rgb = np.flipud(frame_rgb).copy()
colors = frame_rgb.astype(np.float64) / 255.0
valid_colors = colors[valid_mask].reshape(-1, 3)
return valid_colors[display_indices]
def init_open3d_viewer(points, colors):
import open3d as o3d
vis = o3d.visualization.Visualizer()
vis.create_window(window_name="Real-time 3D Point Cloud", width=1280, height=720)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points.astype(np.float64))
pcd.colors = o3d.utility.Vector3dVector(colors.astype(np.float64))
vis.add_geometry(pcd)
render_option = vis.get_render_option()
render_option.point_size = 2.0
render_option.background_color = np.array([0.0, 0.0, 0.0])
vis.poll_events()
vis.update_renderer()
return vis, pcd, o3d
def save_colored_pointcloud(points, colors, output_path):
"""Save displayed points as a PLY file."""
valid_colors = (colors * 255).clip(0, 255).astype(np.uint8)
with open(output_path, "w", encoding="utf-8") as f:
f.write("ply\n")
f.write("format ascii 1.0\n")
f.write(f"element vertex {len(points)}\n")
f.write("property float x\n")
f.write("property float y\n")
f.write("property float z\n")
f.write("property uchar red\n")
f.write("property uchar green\n")
f.write("property uchar blue\n")
f.write("end_header\n")
for pt, col in zip(points, valid_colors):
f.write(f"{pt[0]} {pt[1]} {pt[2]} {col[0]} {col[1]} {col[2]}\n")
print(f"Saved point cloud to {output_path}")
def view_video_pointcloud(
video_path,
data,
display_stride=4,
frame_step=1,
flip_vertical=True,
loop=False,
save_snapshot_path=None,
):
"""Open a video and show it as a real-time colored 3D point cloud."""
cap = cv2.VideoCapture(video_path)
if not cap.isOpened():
raise FileNotFoundError(f"Could not open video: {video_path}")
valid_mask, display_indices, display_points = prepare_static_points(
data,
display_stride=display_stride,
)
ok, frame = cap.read()
if not ok:
cap.release()
raise RuntimeError(f"Could not read first frame from video: {video_path}")
display_colors = frame_to_rgb_colors(
frame,
data,
valid_mask,
display_indices,
flip_vertical=flip_vertical,
)
vis, pcd, o3d = init_open3d_viewer(display_points, display_colors)
print("Real-time point cloud viewer started.")
print("Close the Open3D window or press Ctrl+C in the terminal to stop.")
print(f"Displayed points: {len(display_points)} (stride={display_stride})")
frame_idx = 0
shown_frames = 0
t0 = time.perf_counter()
try:
while True:
alive = vis.poll_events()
if not alive:
break
ok, frame = cap.read()
if not ok:
if loop:
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)
continue
break
frame_idx += 1
if frame_idx % max(1, frame_step) != 0:
continue
display_colors = frame_to_rgb_colors(
frame,
data,
valid_mask,
display_indices,
flip_vertical=flip_vertical,
)
pcd.colors = o3d.utility.Vector3dVector(display_colors)
vis.update_geometry(pcd)
vis.update_renderer()
shown_frames += 1
if shown_frames % 30 == 0:
elapsed = max(time.perf_counter() - t0, 1e-6)
print(f"Displayed {shown_frames} frames, viewer FPS={shown_frames / elapsed:.1f}")
finally:
if save_snapshot_path:
save_colored_pointcloud(display_points, display_colors, save_snapshot_path)
cap.release()
vis.destroy_window()
def parse_args():
parser = argparse.ArgumentParser(description="View a video as a real-time 3D point cloud.")
parser.add_argument("--data-dir", default="./exported_data", help="Blender exported_data folder")
parser.add_argument("--video", default="20260602-094916.mp4", help="Input video path")
parser.add_argument("--display-stride", type=int, default=4,
help="Show every Nth valid point. Increase for faster viewing.")
parser.add_argument("--frame-step", type=int, default=1,
help="Process every Nth video frame.")
parser.add_argument("--no-flip-vertical", action="store_true",
help="Disable vertical flip before color mapping.")
parser.add_argument("--loop", action="store_true", help="Loop video playback.")
parser.add_argument("--save-snapshot", default="",
help="Optional PLY path for saving the last displayed frame.")
return parser.parse_args()
if __name__ == "__main__":
args = parse_args()
exported = load_exported_data(args.data_dir)
view_video_pointcloud(
args.video,
exported,
display_stride=args.display_stride,
frame_step=args.frame_step,
flip_vertical=not args.no_flip_vertical,
loop=args.loop,
save_snapshot_path=args.save_snapshot or None,
)
import numpy as np
import cv2
import os
def load_exported_data(data_dir):
"""Load all exported files."""
# Mesh
mesh = np.load(os.path.join(data_dir, "mesh.npz"))
vertices = mesh['vertices']
faces = mesh['faces']
uv_data = mesh['uv_data']
loop_to_vertex = mesh['loop_to_vertex']
# Camera intrinsics
cam_data = np.load(os.path.join(data_dir, "camera_intrinsics.npz"))
K = cam_data['K']
width = int(cam_data['width'])
height = int(cam_data['height'])
# Camera pose (world-to-camera)
world_to_cam = np.load(os.path.join(data_dir, "world_to_cam.npy"))
# Precomputed pixel-to-3D map (optional)
xyz_map_path = os.path.join(data_dir, "pixel_to_xyz.npz")
if os.path.exists(xyz_map_path):
xyz_map = np.load(xyz_map_path)['xyz_map']
else:
xyz_map = None
return {
'vertices': vertices,
'faces': faces,
'uv_data': uv_data,
'loop_to_vertex': loop_to_vertex,
'K': K,
'width': width,
'height': height,
'world_to_cam': world_to_cam,
'xyz_map': xyz_map
}
def ray_intersect_triangle(ray_origin, ray_dir, v0, v1, v2):
"""Möller–Trumbore intersection algorithm. Returns t (distance) or None."""
epsilon = 1e-8
edge1 = v1 - v0
edge2 = v2 - v0
h = np.cross(ray_dir, edge2)
a = np.dot(edge1, h)
if abs(a) < epsilon:
return None
f = 1.0 / a
s = ray_origin - v0
u = f * np.dot(s, h)
if u < 0.0 or u > 1.0:
return None
q = np.cross(s, edge1)
v = f * np.dot(ray_dir, q)
if v < 0.0 or u + v > 1.0:
return None
t = f * np.dot(edge2, q)
if t > epsilon:
return t
return None
def pixel_to_3d_ray(u, v, K, world_to_cam):
"""
Given pixel coordinates (u,v) and camera intrinsics K,
return a ray in world coordinates: origin and direction.
"""
# Invert K to get normalized image coordinates
fx = K[0,0]; fy = K[1,1]; cx = K[0,2]; cy = K[1,2]
x_norm = (u - cx) / fx
y_norm = (v - cy) / fy
# Ray in camera space: (x_norm, y_norm, 1.0)
ray_cam = np.array([x_norm, y_norm, 1.0], dtype=np.float64)
# Transform to world
cam_to_world = np.linalg.inv(world_to_cam)
R = cam_to_world[:3, :3]
t = cam_to_world[:3, 3]
ray_world_dir = R @ ray_cam
ray_world_dir /= np.linalg.norm(ray_world_dir)
ray_world_origin = t
return ray_world_origin, ray_world_dir
def map_image_to_3d(image_path, data, use_precomputed=True, flip_vertical=True):
"""
Map each pixel of the image to a 3D world point.
If use_precomputed and xyz_map exists, return that directly.
If `flip_vertical` is True, flip the loaded image top-to-bottom so its
row order matches the exported Blender camera/pixel mapping.
Otherwise, perform raycasting on the mesh.
Returns:
points: (H, W, 3) array of 3D points (NaN for misses).
colors: (H, W, 3) array of RGB colors from image.
"""
img = cv2.imread(image_path)
if img is None:
raise FileNotFoundError(f"Could not read image: {image_path}")
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # to RGB
if flip_vertical:
img = np.flipud(img).copy()
h, w = img.shape[:2]
if use_precomputed and data['xyz_map'] is not None:
# Use precomputed map (must match resolution)
if h != data['height'] or w != data['width']:
raise ValueError("Image resolution does not match precomputed map.")
points = data['xyz_map'].copy()
colors = img.astype(np.float32) / 255.0
return points, colors
# Otherwise, perform raycasting on the mesh
vertices = data['vertices']
faces = data['faces']
K = data['K']
world_to_cam = data['world_to_cam']
points = np.full((h, w, 3), np.nan, dtype=np.float32)
colors = img.astype(np.float32) / 255.0
# For speed, we can use a bounding volume hierarchy (BVH).
# Here we implement a simple brute-force over all faces (slow for large meshes).
# In practice, precomputing the xyz_map is strongly recommended.
for y in range(h):
for x in range(w):
origin, direction = pixel_to_3d_ray(x, y, K, world_to_cam)
best_t = np.inf
best_hit = None
for face in faces:
v0 = vertices[face[0]]
v1 = vertices[face[1]]
v2 = vertices[face[2]]
t = ray_intersect_triangle(origin, direction, v0, v1, v2)
if t is not None and t < best_t:
best_t = t
best_hit = origin + direction * t
if best_hit is not None:
points[y, x] = best_hit
return points, colors
def save_colored_pointcloud(points, colors, output_path):
"""Save as PLY file (colored point cloud)."""
# Flatten valid points
mask = ~np.isnan(points[:,:,0])
valid_points = points[mask].reshape(-1, 3)
valid_colors = (colors[mask].reshape(-1, 3) * 255).astype(np.uint8)
with open(output_path, 'w') as f:
f.write("ply\n")
f.write("format ascii 1.0\n")
f.write(f"element vertex {len(valid_points)}\n")
f.write("property float x\n")
f.write("property float y\n")
f.write("property float z\n")
f.write("property uchar red\n")
f.write("property uchar green\n")
f.write("property uchar blue\n")
f.write("end_header\n")
for pt, col in zip(valid_points, valid_colors):
f.write(f"{pt[0]} {pt[1]} {pt[2]} {col[0]} {col[1]} {col[2]}\n")
print(f"Saved point cloud to {output_path}")
# if __name__ == "__main__":
# data_dir = "./exported_data" # Path to exported files
# image_path = "captured_image_1.jpg" # Input camera image
# output_ply = "mapped_points.ply"
# data = load_exported_data(data_dir)
# points, colors = map_image_to_3d(
# image_path,
# data,
# use_precomputed=True,
# flip_vertical=True,
# )
# save_colored_pointcloud(points, colors, output_ply)
# # Optional: Visualize using Open3D
# import open3d as o3d
# pcd = o3d.geometry.PointCloud()
# mask = ~np.isnan(points[:,:,0])
# pcd.points = o3d.utility.Vector3dVector(points[mask].reshape(-1,3))
# pcd.colors = o3d.utility.Vector3dVector(colors[mask].reshape(-1,3))
# o3d.visualization.draw_geometries([pcd])
Binary file not shown.

Before

Width:  |  Height:  |  Size: 272 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 662 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 635 KiB

@@ -1,917 +0,0 @@
"""
Texture Generator for Arbitrary 3D Surfaces (Blender Script)
- Mode: "patches" (colored patches), "checker" (black/white patches),
or "dots" (circular markers)
- Mesh UVs are filled white so that non‑visible parts become white.
- Texture background (no‑mesh areas) stays black.
"""
import bpy
import numpy as np
import os
import time
from mathutils import Vector
# ----------------------------------------------------------------------
# Common Utilities
# ----------------------------------------------------------------------
def draw_cross(pixels, width, height, cx_uv, cy_uv, length_px, thickness_px, color):
"""
Draw a cross (plus sign) in pixel space.
cx_uv, cy_uv: centre of cross in UV coordinates (0‑1), origin at bottom‑left.
length_px: total length of each arm (pixels).
thickness_px: line thickness (odd number gives symmetric result).
"""
# Convert UV centre to pixel coordinates (y from bottom)
cx = int(round(cx_uv * (width - 1)))
cy = int(round(cy_uv * (height - 1)))
half_len = length_px // 2
half_thick = thickness_px // 2
# Horizontal arm
x0 = max(0, cx - half_len)
x1 = min(width - 1, cx + half_len)
y0 = max(0, cy - half_thick)
y1 = min(height - 1, cy + half_thick)
for y in range(y0, y1 + 1):
for x in range(x0, x1 + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
# Vertical arm
x0 = max(0, cx - half_thick)
x1 = min(width - 1, cx + half_thick)
y0 = max(0, cy - half_len)
y1 = min(height - 1, cy + half_len)
for y in range(y0, y1 + 1):
for x in range(x0, x1 + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def fill_background(pixels, width, height, color):
"""Fill the entire image with a solid RGBA color (0-1 floats)."""
r, g, b, a = color
for y in range(height):
for x in range(width):
idx = (y * width + x) * 4
pixels[idx] = r
pixels[idx+1] = g
pixels[idx+2] = b
pixels[idx+3] = a
def fill_polygon(image, pixels, width, height, polygon, color):
"""Fill a polygon in the pixel array (flattened, RGBA)."""
poly = [(int(round(x)), int(round(y))) for x, y in polygon]
if len(poly) < 3:
return
min_x = max(0, min(p[0] for p in poly))
max_x = min(width-1, max(p[0] for p in poly))
min_y = max(0, min(p[1] for p in poly))
max_y = min(height-1, max(p[1] for p in poly))
for y in range(min_y, max_y + 1):
intersections = []
for i in range(len(poly)):
x1, y1 = poly[i]
x2, y2 = poly[(i+1) % len(poly)]
if y1 == y2:
continue
if (y1 <= y < y2) or (y2 <= y < y1):
t = (y - y1) / (y2 - y1)
x = x1 + t * (x2 - x1)
intersections.append(x)
intersections.sort()
for k in range(0, len(intersections), 2):
if k+1 >= len(intersections):
break
x_start = int(round(intersections[k]))
x_end = int(round(intersections[k+1]))
x_start = max(0, min(width-1, x_start))
x_end = max(0, min(width-1, x_end))
for x in range(x_start, x_end + 1):
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def fill_polygon_value(grid, width, height, polygon, value):
"""Fill a polygon into a 2D scalar grid."""
poly = [(int(round(x)), int(round(y))) for x, y in polygon]
if len(poly) < 3:
return
min_x = max(0, min(p[0] for p in poly))
max_x = min(width - 1, max(p[0] for p in poly))
min_y = max(0, min(p[1] for p in poly))
max_y = min(height - 1, max(p[1] for p in poly))
for y in range(min_y, max_y + 1):
intersections = []
for i in range(len(poly)):
x1, y1 = poly[i]
x2, y2 = poly[(i + 1) % len(poly)]
if y1 == y2:
continue
if (y1 <= y < y2) or (y2 <= y < y1):
t = (y - y1) / (y2 - y1)
x = x1 + t * (x2 - x1)
intersections.append(x)
intersections.sort()
for k in range(0, len(intersections), 2):
if k + 1 >= len(intersections):
break
x_start = int(round(intersections[k]))
x_end = int(round(intersections[k + 1]))
x_start = max(0, min(width - 1, x_start))
x_end = max(0, min(width - 1, x_end))
grid[y, x_start:x_end + 1] = value
def unwrap_periodic_values(values):
"""
Unwrap periodic UV coordinates so points straddling the 0/1 seam stay local.
"""
arr = np.asarray(values, dtype=np.float64).copy()
if arr.size == 0:
return arr
if arr.max() - arr.min() > 0.5:
arr[arr < 0.5] += 1.0
return arr
def wrap_uv_polygon_for_raster(uv_poly):
"""
Return a seam-aware version of a UV polygon.
UV islands may cross the 0/1 seam. If we rasterize those raw coordinates as a
single polygon, Blender's repeated UVs turn into a long diagonal triangle
across the texture. We unwrap seam-crossing polygons first, then draw wrapped
copies in the texture tile(s) they overlap.
"""
if len(uv_poly) < 3:
return []
uv_arr = np.array([(float(u), float(v)) for u, v in uv_poly], dtype=np.float64)
uv_arr[:, 0] = unwrap_periodic_values(uv_arr[:, 0])
uv_arr[:, 1] = unwrap_periodic_values(uv_arr[:, 1])
return [tuple(pt) for pt in uv_arr]
def uv_polygon_to_pixel_polygons(uv_poly, tex_width, tex_height):
"""
Convert one UV polygon into one or more pixel polygons, duplicating seam-
crossing polygons into adjacent tiles so rasterization stays local.
"""
unwrapped = wrap_uv_polygon_for_raster(uv_poly)
if len(unwrapped) < 3:
return []
u_vals = [u for u, _ in unwrapped]
v_vals = [v for _, v in unwrapped]
u_shift_start = int(np.floor(min(u_vals)))
u_shift_end = int(np.floor(max(u_vals)))
v_shift_start = int(np.floor(min(v_vals)))
v_shift_end = int(np.floor(max(v_vals)))
pixel_polys = []
for u_shift in range(u_shift_start, u_shift_end + 1):
for v_shift in range(v_shift_start, v_shift_end + 1):
shifted = [(u - u_shift, v - v_shift) for u, v in unwrapped]
shifted_u = [u for u, _ in shifted]
shifted_v = [v for _, v in shifted]
if max(shifted_u) < 0.0 or min(shifted_u) > 1.0:
continue
if max(shifted_v) < 0.0 or min(shifted_v) > 1.0:
continue
pixel_polys.append([
(u * (tex_width - 1), v * (tex_height - 1))
for u, v in shifted
])
return pixel_polys
def draw_ellipse(pixels, width, height, cx, cy, a, b, angle, color):
"""Filled ellipse with semi‑axes a,b (pixels), rotation angle (rad)."""
cos_a = np.cos(angle)
sin_a = np.sin(angle)
rx = abs(a * cos_a) + abs(b * sin_a)
ry = abs(a * sin_a) + abs(b * cos_a)
x_min = int(max(0, cx - rx))
x_max = int(min(width-1, cx + rx))
y_min = int(max(0, cy - ry))
y_max = int(min(height-1, cy + ry))
for y in range(y_min, y_max + 1):
for x in range(x_min, x_max + 1):
dx = x - cx
dy = y - cy
u = dx * cos_a + dy * sin_a
v = -dx * sin_a + dy * cos_a
if (u/a)**2 + (v/b)**2 <= 1.0:
idx = (y * width + x) * 4
pixels[idx:idx+4] = color
def get_uv_from_face(obj, face_index, location):
"""World‑space location → UV coordinate on the mesh face."""
me = obj.data
local_loc = obj.matrix_world.inverted() @ location
poly = me.polygons[face_index]
vert_indices = poly.vertices
loop_indices = poly.loop_indices
if len(vert_indices) < 3:
return None
v0 = me.vertices[vert_indices[0]].co
v1 = me.vertices[vert_indices[1]].co
v2 = me.vertices[vert_indices[2]].co
uv_layer = me.uv_layers.active.data
uv0 = uv_layer[loop_indices[0]].uv
uv1 = uv_layer[loop_indices[1]].uv
uv2 = uv_layer[loop_indices[2]].uv
def tri_area(a, b, c):
return 0.5 * ((b - a).cross(c - a)).length
total_area = tri_area(v0, v1, v2)
if total_area == 0:
return None
w0 = tri_area(local_loc, v1, v2) / total_area
w1 = tri_area(local_loc, v2, v0) / total_area
w2 = 1.0 - w0 - w1
uv = w0 * uv0 + w1 * uv1 + w2 * uv2
return Vector((uv.x, uv.y))
def raycast_from_camera(scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height):
"""Camera raycast through normalised sensor coordinates (0‑1)."""
sx = (u_ndc - 0.5) * frame_width
sy = (v_ndc - 0.5) * frame_height
direction_cam = Vector((sx, sy, -1.0))
direction_world = cam.matrix_world.to_quaternion() @ direction_cam
direction_world.normalize()
origin = cam.matrix_world.translation
hit, location, normal, face_index, hit_obj, matrix = scene.ray_cast(
depsgraph, origin, direction_world)
if hit and hit_obj == obj:
return location, face_index
return None, None
def raycast_uv_from_camera(scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height):
"""Camera raycast through normalized sensor coordinates and return UV."""
loc, face_idx = raycast_from_camera(
scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height
)
if loc is None:
return None
return get_uv_from_face(obj, face_idx, loc)
def convex_hull_uv(points):
"""Return UV points ordered around their 2D convex hull."""
unwrapped = wrap_uv_polygon_for_raster([(p.x, p.y) for p in points])
unique = sorted(set(unwrapped))
if len(unique) < 3:
return None
def cross(o, a, b):
return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
lower = []
for p in unique:
while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
lower.pop()
lower.append(p)
upper = []
for p in reversed(unique):
while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
upper.pop()
upper.append(p)
hull = lower[:-1] + upper[:-1]
if len(hull) < 3:
return None
return [Vector((u, v)) for u, v in hull]
def sample_patch_uv_polygon(
scene, depsgraph, cam, obj,
u_start, u_end, v_start, v_end,
frame_width, frame_height,
edge_sample_grid=5
):
"""
Build a UV polygon for one camera-space patch.
Full interior cells use the four corners. Edge cells that only partially hit
the mesh are sampled more densely and converted to a convex hull so they are
not dropped just because some corners miss the silhouette.
"""
corners_ndc = [
(u_start, v_start),
(u_end, v_start),
(u_end, v_end),
(u_start, v_end),
]
corner_uvs = []
for u_ndc, v_ndc in corners_ndc:
uv = raycast_uv_from_camera(
scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height
)
if uv is not None:
corner_uvs.append(uv)
if len(corner_uvs) >= 3:
return corner_uvs
sampled_uvs = list(corner_uvs)
for sy in range(edge_sample_grid):
v_ndc = v_start + (v_end - v_start) * sy / (edge_sample_grid - 1)
for sx in range(edge_sample_grid):
u_ndc = u_start + (u_end - u_start) * sx / (edge_sample_grid - 1)
uv = raycast_uv_from_camera(
scene, depsgraph, cam, obj, u_ndc, v_ndc, frame_width, frame_height
)
if uv is not None:
sampled_uvs.append(uv)
return convex_hull_uv(sampled_uvs)
def rasterize_mesh_uvs(obj, pixels, tex_width, tex_height, color=(1,1,1,1)):
"""
Fill all triangles of the active UV map with 'color' (default white).
This marks every part of the mesh, visible or not, as white.
"""
me = obj.data
uv_layer = me.uv_layers.active
if not uv_layer:
print("Warning: No active UV layer – mesh rasterization skipped.")
return
uv_data = uv_layer.data
for poly in me.polygons:
# Get loop indices for this polygon
loops = poly.loop_indices
# Build UV polygon
uv_poly = []
for loop_idx in loops:
uv = uv_data[loop_idx].uv
uv_poly.append((uv.x, uv.y))
# Convert UV (0‑1) to pixel coordinates
for pixel_poly in uv_polygon_to_pixel_polygons(uv_poly, tex_width, tex_height):
fill_polygon(None, pixels, tex_width, tex_height, pixel_poly, color)
def rasterize_mesh_mask(obj, tex_width, tex_height):
"""Rasterize the active UV map to a boolean mesh mask."""
mask = np.zeros((tex_height, tex_width), dtype=bool)
me = obj.data
uv_layer = me.uv_layers.active
if not uv_layer:
return mask
uv_data = uv_layer.data
for poly in me.polygons:
uv_poly = []
for loop_idx in poly.loop_indices:
uv = uv_data[loop_idx].uv
uv_poly.append((uv.x, uv.y))
for pixel_poly in uv_polygon_to_pixel_polygons(uv_poly, tex_width, tex_height):
fill_polygon_value(mask, tex_width, tex_height, pixel_poly, True)
return mask
def paint_mask_color(pixels, width, height, mask, color):
"""Paint all True pixels of a boolean mask with an RGBA color."""
r, g, b, a = color
ys, xs = np.nonzero(mask)
for y, x in zip(ys, xs):
idx = (y * width + x) * 4
pixels[idx:idx+4] = (r, g, b, a)
def complete_checker_gaps(pixels, width, height, mesh_mask, patch_state, max_gap_px=6):
"""
Fill small white slivers inside the mesh by expanding nearby checker colors.
`patch_state`: -1 = unassigned, 0 = black patch, 1 = white patch.
"""
for _ in range(max_gap_px):
candidate = mesh_mask & (patch_state < 0)
if not np.any(candidate):
break
black_counts = np.zeros_like(patch_state, dtype=np.uint8)
white_counts = np.zeros_like(patch_state, dtype=np.uint8)
for dy in (-1, 0, 1):
for dx in (-1, 0, 1):
if dx == 0 and dy == 0:
continue
src_y0 = max(0, -dy)
src_y1 = height - max(0, dy)
src_x0 = max(0, -dx)
src_x1 = width - max(0, dx)
dst_y0 = max(0, dy)
dst_y1 = height - max(0, -dy)
dst_x0 = max(0, dx)
dst_x1 = width - max(0, -dx)
neighbor_vals = patch_state[src_y0:src_y1, src_x0:src_x1]
black_counts[dst_y0:dst_y1, dst_x0:dst_x1] += (neighbor_vals == 0)
white_counts[dst_y0:dst_y1, dst_x0:dst_x1] += (neighbor_vals == 1)
fillable = candidate & ((black_counts + white_counts) > 0)
if not np.any(fillable):
break
patch_state[fillable] = np.where(
black_counts[fillable] > white_counts[fillable], 0, 1
).astype(np.int8)
paint_mask_color(pixels, width, height, patch_state == 0, (0.0, 0.0, 0.0, 1.0))
paint_mask_color(pixels, width, height, patch_state == 1, (1.0, 1.0, 1.0, 1.0))
# ----------------------------------------------------------------------
# Color generation (for patch mode)
# ----------------------------------------------------------------------
def merge_intervals(intervals):
if not intervals: return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
merged = merge_intervals(forbidden)
allowed = []
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0: return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed: return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def checker_patch_color(i, j):
"""Return a black/white checker color for patch index (i, j)."""
if (i + j) % 2 == 0:
return np.array([1.0, 1.0, 1.0], dtype=np.float32)
return np.array([0.0, 0.0, 0.0], dtype=np.float32)
# ----------------------------------------------------------------------
# Mode 1: Colored Patches
# ----------------------------------------------------------------------
def generate_patch_texture(
obj, cam, scene, depsgraph,
nx=40, ny=30, r=0.05,
texture_width=2048, texture_height=2048,
output_path="//patch_texture.png", uv_centre=None, cross_length=40,
cross_thickness=3, cross_color=(1,0,0,1), edge_sample_grid=5
):
cam_data = cam.data
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
if cam_data.sensor_fit == 'HORIZONTAL':
sensor_width = cam_data.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam_data.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam_data.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
edge_sample_grid = max(2, int(edge_sample_grid))
patch_uv_polygons = []
patch_colors = np.zeros((ny, nx, 3), dtype=np.float32)
skipped_patches = 0
print(f"Raycasting {nx*ny} patches...")
for j in range(ny):
v_start = j / ny
v_end = (j+1) / ny
for i in range(nx):
u_start = i / nx
u_end = (i+1) / nx
uv_poly = sample_patch_uv_polygon(
scene, depsgraph, cam, obj,
u_start, u_end, v_start, v_end,
frame_width, frame_height,
edge_sample_grid=edge_sample_grid
)
if uv_poly is None:
patch_uv_polygons.append(None)
skipped_patches += 1
continue
patch_uv_polygons.append(uv_poly)
neighbors = []
if j > 0: neighbors.append(patch_colors[j-1, i])
if i > 0: neighbors.append(patch_colors[j, i-1])
color = np.random.rand(3) if len(neighbors)==0 else generate_valid_color(neighbors, r)
patch_colors[j, i] = color
print(f"Skipped {skipped_patches} patches with fewer than 3 mesh samples.")
# Create image – start with black background
img_name = "PatchTexture"
if img_name in bpy.data.images:
bpy.data.images.remove(bpy.data.images[img_name])
image = bpy.data.images.new(img_name, width=texture_width, height=texture_height, alpha=True)
pixels = list(image.pixels)
# Fill entire image with black (no mesh areas)
fill_background(pixels, texture_width, texture_height, (0,0,0,1))
# Rasterize all mesh UVs in white (visible + non‑visible parts)
rasterize_mesh_uvs(obj, pixels, texture_width, texture_height, color=(1,1,1,1))
# Now draw the visible patches on top
idx = 0
for j in range(ny):
for i in range(nx):
uv_poly = patch_uv_polygons[idx]; idx += 1
if uv_poly is None: continue
color = patch_colors[j,i]
rgba = (color[0], color[1], color[2], 1.0)
for pixel_poly in uv_polygon_to_pixel_polygons(uv_poly, texture_width, texture_height):
fill_polygon(image, pixels, texture_width, texture_height, pixel_poly, rgba)
# (after drawing all patches)
if uv_centre is not None:
draw_cross(pixels, texture_width, texture_height,
uv_centre[0], uv_centre[1],
cross_length, cross_thickness, cross_color)
image.pixels[:] = pixels
image.update()
abs_path = bpy.path.abspath(output_path)
image.filepath_raw = abs_path
image.file_format = 'PNG'
image.save()
print(f"Patch texture saved to {abs_path}")
return image
# ----------------------------------------------------------------------
# Mode 2: Black/White Checker Patches
# ----------------------------------------------------------------------
def generate_checker_texture(
obj, cam, scene, depsgraph,
nx=40, ny=40,
texture_width=2048, texture_height=2048,
output_path="//checker_texture.png", uv_centre=None, cross_length=40,
cross_thickness=3, cross_color=(1,0,0,1), edge_sample_grid=5,
gap_fill_px=6
):
cam_data = cam.data
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
if cam_data.sensor_fit == 'HORIZONTAL':
sensor_width = cam_data.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam_data.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam_data.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
edge_sample_grid = max(2, int(edge_sample_grid))
patch_uv_polygons = []
patch_colors = np.zeros((ny, nx, 3), dtype=np.float32)
skipped_patches = 0
print(f"Raycasting {nx*ny} checker patches...")
for j in range(ny):
v_start = j / ny
v_end = (j+1) / ny
for i in range(nx):
u_start = i / nx
u_end = (i+1) / nx
uv_poly = sample_patch_uv_polygon(
scene, depsgraph, cam, obj,
u_start, u_end, v_start, v_end,
frame_width, frame_height,
edge_sample_grid=edge_sample_grid
)
if uv_poly is None:
patch_uv_polygons.append(None)
skipped_patches += 1
continue
patch_uv_polygons.append(uv_poly)
patch_colors[j, i] = checker_patch_color(i, j)
print(f"Skipped {skipped_patches} checker patches with fewer than 3 mesh samples.")
img_name = "CheckerTexture"
if img_name in bpy.data.images:
bpy.data.images.remove(bpy.data.images[img_name])
image = bpy.data.images.new(img_name, width=texture_width, height=texture_height, alpha=True)
pixels = list(image.pixels)
mesh_mask = rasterize_mesh_mask(obj, texture_width, texture_height)
patch_state = np.full((texture_height, texture_width), -1, dtype=np.int8)
fill_background(pixels, texture_width, texture_height, (0,0,0,1))
paint_mask_color(pixels, texture_width, texture_height, mesh_mask, (1,1,1,1))
idx = 0
for j in range(ny):
for i in range(nx):
uv_poly = patch_uv_polygons[idx]; idx += 1
if uv_poly is None:
continue
color = patch_colors[j,i]
rgba = (color[0], color[1], color[2], 1.0)
patch_value = 1 if np.mean(color) >= 0.5 else 0
for pixel_poly in uv_polygon_to_pixel_polygons(uv_poly, texture_width, texture_height):
fill_polygon(image, pixels, texture_width, texture_height, pixel_poly, rgba)
fill_polygon_value(
patch_state,
texture_width,
texture_height,
pixel_poly,
patch_value,
)
complete_checker_gaps(
pixels, texture_width, texture_height, mesh_mask, patch_state, max_gap_px=gap_fill_px
)
if uv_centre is not None:
draw_cross(pixels, texture_width, texture_height,
uv_centre[0], uv_centre[1],
cross_length, cross_thickness, cross_color)
image.pixels[:] = pixels
image.update()
abs_path = bpy.path.abspath(output_path)
image.filepath_raw = abs_path
image.file_format = 'PNG'
image.save()
print(f"Checker texture saved to {abs_path}")
return image
# ----------------------------------------------------------------------
# Mode 3: Circular Dot Markers
# ----------------------------------------------------------------------
def generate_dot_texture(
obj, cam, scene, depsgraph,
dot_diameter_pixels=12,
spacing_pixels=40,
dot_color=(0,0,0,1), # RGBA 0-1
texture_width=2048, texture_height=2048,
output_path="//dot_texture.png",uv_centre=None, cross_length=40, cross_thickness=3, cross_color=(1,0,0,1)
):
render = scene.render
res_x = render.resolution_x
res_y = render.resolution_y
cam_data = cam.data
if cam_data.sensor_fit == 'HORIZONTAL':
sensor_width = cam_data.sensor_width
sensor_height = sensor_width / (res_x / res_y)
else:
sensor_height = cam_data.sensor_height
sensor_width = sensor_height * (res_x / res_y)
focal = cam_data.lens
frame_width = sensor_width / focal
frame_height = sensor_height / focal
margin = spacing_pixels // 2
x_centers = np.arange(margin, res_x, spacing_pixels)
y_centers = np.arange(margin, res_y, spacing_pixels)
offset_pix = max(2, dot_diameter_pixels // 4)
du_ndc = offset_pix / res_x
dv_ndc = offset_pix / res_y
dots = [] # (cx_uv, cy_uv, a_uv, b_uv, angle_uv)
print(f"Raycasting {len(x_centers)*len(y_centers)} dot centers...")
for y_pix in y_centers:
v_ndc = y_pix / res_y
for x_pix in x_centers:
u_ndc = x_pix / res_x
center_loc, face_idx = raycast_from_camera(scene, depsgraph, cam, obj,
u_ndc, v_ndc, frame_width, frame_height)
if center_loc is None: continue
center_uv = get_uv_from_face(obj, face_idx, center_loc)
if center_uv is None: continue
offsets = [
( u_ndc + du_ndc, v_ndc),
( u_ndc - du_ndc, v_ndc),
( u_ndc, v_ndc + dv_ndc),
( u_ndc, v_ndc - dv_ndc)
]
uv_offsets = []
for (u_off, v_off) in offsets:
uv_off = raycast_uv_from_camera(
scene, depsgraph, cam, obj,
u_off, v_off, frame_width, frame_height
)
uv_offsets.append(uv_off)
scale_u = 1.0 / res_x
scale_v = 1.0 / res_y
r_pix = dot_diameter_pixels / 2.0
right_uv, left_uv, up_uv, down_uv = uv_offsets
if right_uv is not None and left_uv is not None:
du_dx = (right_uv - left_uv) / (2 * offset_pix * scale_u)
a_uv = r_pix * du_dx.length * scale_u
elif right_uv is not None:
du_dx = (right_uv - center_uv) / (offset_pix * scale_u)
a_uv = r_pix * du_dx.length * scale_u
elif left_uv is not None:
du_dx = (center_uv - left_uv) / (offset_pix * scale_u)
a_uv = r_pix * du_dx.length * scale_u
else:
du_dx = Vector((1.0, 0.0))
a_uv = r_pix / texture_width
if up_uv is not None and down_uv is not None:
dv_dx = (up_uv - down_uv) / (2 * offset_pix * scale_v)
b_uv = r_pix * dv_dx.length * scale_v
elif up_uv is not None:
dv_dx = (up_uv - center_uv) / (offset_pix * scale_v)
b_uv = r_pix * dv_dx.length * scale_v
elif down_uv is not None:
dv_dx = (center_uv - down_uv) / (offset_pix * scale_v)
b_uv = r_pix * dv_dx.length * scale_v
else:
b_uv = r_pix / texture_height
a_uv = max(a_uv, 1.0 / texture_width)
b_uv = max(b_uv, 1.0 / texture_height)
angle = np.arctan2(du_dx.y, du_dx.x) if du_dx.length > 0 else 0.0
dots.append((center_uv.x, center_uv.y, a_uv, b_uv, angle))
print(f"Found {len(dots)} valid dots.")
# Create image
img_name = "DotTexture"
if img_name in bpy.data.images:
bpy.data.images.remove(bpy.data.images[img_name])
image = bpy.data.images.new(img_name, width=texture_width, height=texture_height, alpha=True)
pixels = list(image.pixels)
# Black background (no mesh)
fill_background(pixels, texture_width, texture_height, (0,0,0,1))
# White for all mesh areas
rasterize_mesh_uvs(obj, pixels, texture_width, texture_height, color=(1,1,1,1))
# Draw dots
for (cx_uv, cy_uv, a_uv, b_uv, angle_uv) in dots:
cx_px = cx_uv * (texture_width - 1)
cy_px = cy_uv * (texture_height - 1)
a_px = a_uv * texture_width
b_px = b_uv * texture_height
draw_ellipse(pixels, texture_width, texture_height,
cx_px, cy_px, a_px, b_px, angle_uv, dot_color)
if uv_centre is not None:
draw_cross(pixels, texture_width, texture_height,
uv_centre[0], uv_centre[1],
cross_length, cross_thickness, cross_color)
image.pixels[:] = pixels
image.update()
abs_path = bpy.path.abspath(output_path)
image.filepath_raw = abs_path
image.file_format = 'PNG'
image.save()
print(f"Dot texture saved to {abs_path}")
return image
# ----------------------------------------------------------------------
# Main
# ----------------------------------------------------------------------
def main():
# ==================== CONFIGURATION ====================
MODE = "patches" # "patches", "checker", or "dots"
OBJ_NAME = "FingerMesh0525.001"
CAM_NAME = "Camera0607"
TEXTURE_WIDTH = 1640*2
TEXTURE_HEIGHT = 1600*2
OUTPUT_NAMES = {
"patches": "patch_texture.png",
"checker": "checker_texture.png",
"dots": "dot_texture.png",
}
OUTPUT_PATH = "//" + OUTPUT_NAMES.get(MODE, "texture.png")
# Patch mode
PATCH_NX = 80; PATCH_NY = 45; COLOR_R = 0.05
EDGE_SAMPLE_GRID = 5
# Checker mode
CHECKER_NX = 40; CHECKER_NY = 40; CHECKER_GAP_FILL_PX = 6
# Dot mode
DOT_DIAMETER_PIX = 12; DOT_SPACING_PIX = 40; DOT_COLOR = (0,0,0,1)
# Cross
CROSS_LENGTH_PX = 40
CROSS_THICKNESS_PX = 3
CROSS_COLOR = (1, 0, 0, 1) # red cross
# ======================================================
start_time = time.time()
scene = bpy.context.scene
depsgraph = bpy.context.evaluated_depsgraph_get()
if OBJ_NAME:
obj = bpy.data.objects.get(OBJ_NAME)
else:
obj = bpy.context.active_object
if obj is None or obj.type != 'MESH':
raise ValueError("No valid mesh object found.")
cam = bpy.data.objects.get(CAM_NAME) if CAM_NAME else scene.camera
if cam is None or cam.type != 'CAMERA':
raise ValueError("No camera found.")
if not obj.data.uv_layers.active:
obj.data.uv_layers.new(name="UVMap")
# Compute UV centre
uv_layer = obj.data.uv_layers.active
uv_coords = np.array([uv_layer.data[i].uv.to_tuple() for i in range(len(uv_layer.data))])
uv_min = uv_coords.min(axis=0)
uv_max = uv_coords.max(axis=0)
uv_centre = (uv_min + uv_max) / 2.0
print(f"UV centre: ({uv_centre[0]:.4f}, {uv_centre[1]:.4f})")
# Compute UV bounding box centre
me = obj.data
uv_layer = me.uv_layers.active
uv_coords = np.array([uv_layer.data[i].uv.to_tuple() for i in range(len(uv_layer.data))])
uv_min = uv_coords.min(axis=0)
uv_max = uv_coords.max(axis=0)
uv_centre = (uv_min + uv_max) / 2.0
print(f"UV centre: {uv_centre[0]:.4f}, {uv_centre[1]:.4f}")
if MODE == "patches":
generate_patch_texture(
obj, cam, scene, depsgraph,
nx=PATCH_NX, ny=PATCH_NY, r=COLOR_R,
texture_width=TEXTURE_WIDTH, texture_height=TEXTURE_HEIGHT,
output_path=OUTPUT_PATH,
uv_centre=uv_centre,
cross_length=CROSS_LENGTH_PX,
cross_thickness=CROSS_THICKNESS_PX,
cross_color=CROSS_COLOR,
edge_sample_grid=EDGE_SAMPLE_GRID
)
elif MODE == "checker":
generate_checker_texture(
obj, cam, scene, depsgraph,
nx=CHECKER_NX, ny=CHECKER_NY,
texture_width=TEXTURE_WIDTH, texture_height=TEXTURE_HEIGHT,
output_path=OUTPUT_PATH,
uv_centre=uv_centre,
cross_length=CROSS_LENGTH_PX,
cross_thickness=CROSS_THICKNESS_PX,
cross_color=CROSS_COLOR,
edge_sample_grid=EDGE_SAMPLE_GRID,
gap_fill_px=CHECKER_GAP_FILL_PX
)
elif MODE == "dots":
generate_dot_texture(
obj, cam, scene, depsgraph,
dot_diameter_pixels=DOT_DIAMETER_PIX,
spacing_pixels=DOT_SPACING_PIX,
dot_color=DOT_COLOR,
texture_width=TEXTURE_WIDTH, texture_height=TEXTURE_HEIGHT,
output_path=OUTPUT_PATH,
uv_centre=uv_centre,
cross_length=CROSS_LENGTH_PX,
cross_thickness=CROSS_THICKNESS_PX,
cross_color=CROSS_COLOR
)
else:
raise ValueError(f"Unknown MODE: {MODE!r}. Use 'patches', 'checker', or 'dots'.")
print(f"Completed in {time.time()-start_time:.2f} seconds.")
if __name__ == "__main__":
main()
File diff suppressed because it is too large Load Diff
Binary file not shown.

Before

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 1.3 MiB

@@ -1,720 +0,0 @@
#!/usr/bin/env python3
"""
Real-time UV-based perspective correction for the exported finger mesh.
This version replaces the previous cylindrical approximation with a remap built
from the mesh UV parameterization. The remap is precomputed once and then
applied to each live frame with cv2.remap().
"""
import argparse
import sys
import time
import cv2
import numpy as np
import orisys
from orisys.configs import load_config
from orisys.frame_processors import FrameUnwarper
from orisys.unroll import (
barycentric_grid,
build_or_load_domain,
compute_uv_distortion_metrics,
face_uvs_from_loop_data,
load_calib_json,
load_exported_data,
recommend_edge_expand,
resolve_uv_runtime_sizes,
scale_camera_geometry,
scaled_map_cache_path,
uv_triangle_to_image,
)
class StepTimer:
"""Accumulate per-step wall times for rolling runtime averages."""
def __init__(self):
self._starts = {}
self.totals = {}
self.counts = {}
def start(self, name):
self._starts[name] = time.perf_counter()
def stop(self, name):
dt = time.perf_counter() - self._starts.pop(name)
self.totals[name] = self.totals.get(name, 0.0) + dt
self.counts[name] = self.counts.get(name, 0) + 1
return dt
def avg_ms(self, name):
if self.counts.get(name, 0) == 0:
return 0.0
return 1000.0 * self.totals[name] / self.counts[name]
def benchmark_step(label, func, *args, **kwargs):
"""Run func once and print its wall time."""
t0 = time.perf_counter()
result = func(*args, **kwargs)
ms = (time.perf_counter() - t0) * 1000.0
print(f"[benchmark] {label}: {ms:.1f} ms")
return result
def print_runtime_benchmark(timer):
"""Format rolling average timings for the live loop."""
parts = [
f"capture={timer.avg_ms('capture'):.1f}ms",
f"uv_prep={timer.avg_ms('uv_prepare'):.1f}ms",
f"uv_remap={timer.avg_ms('uv_remap'):.1f}ms",
f"get_img={timer.avg_ms('get_img'):.1f}ms",
f"flow={timer.avg_ms('flow'):.1f}ms",
f"hhd={timer.avg_ms('hhd'):.1f}ms",
f"contact={timer.avg_ms('contact'):.1f}ms",
f"viz={timer.avg_ms('viz'):.1f}ms",
]
return "bench[" + ", ".join(parts) + "]"
def resize_to_height(image, target_height):
"""Resize an image to a target height while preserving aspect ratio."""
if image.shape[0] == target_height:
return image
scale = target_height / float(image.shape[0])
target_width = max(1, int(round(image.shape[1] * scale)))
return cv2.resize(image, (target_width, target_height), interpolation=cv2.INTER_AREA)
def resize_to_fit(image, target_width, target_height):
"""Resize with aspect ratio preserved and pad to exact tile size."""
h, w = image.shape[:2]
if h <= 0 or w <= 0:
return np.zeros((target_height, target_width, 3), dtype=np.uint8)
scale = min(target_width / float(w), target_height / float(h))
new_w = max(1, int(round(w * scale)))
new_h = max(1, int(round(h * scale)))
resized = cv2.resize(image, (new_w, new_h), interpolation=cv2.INTER_AREA)
if resized.ndim == 2:
resized = cv2.cvtColor(resized, cv2.COLOR_GRAY2BGR)
elif resized.ndim == 3 and resized.shape[2] == 1:
resized = cv2.cvtColor(resized[:, :, 0], cv2.COLOR_GRAY2BGR)
canvas = np.zeros((target_height, target_width, 3), dtype=np.uint8)
y0 = (target_height - new_h) // 2
x0 = (target_width - new_w) // 2
canvas[y0:y0 + new_h, x0:x0 + new_w] = resized
return canvas
def add_panel_title(image, title):
"""Overlay a simple title onto a display tile."""
out = image.copy()
cv2.rectangle(out, (0, 0), (out.shape[1], 36), (0, 0, 0), thickness=-1)
cv2.putText(out, title, (12, 24), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (255, 255, 255), 2, cv2.LINE_AA)
return out
def apply_display_mask(image, valid_mask, border_value=0):
"""Zero invalid UV pixels in visualization panels."""
if image is None or valid_mask is None:
return image
out = image.copy()
mask = valid_mask.astype(np.uint8)
if mask.shape != out.shape[:2]:
mask = cv2.resize(mask, (out.shape[1], out.shape[0]), interpolation=cv2.INTER_NEAREST)
invalid = ~mask.astype(bool)
out[invalid] = border_value
return out
def contact_inside_mask(centroid, valid_mask):
if centroid is None or valid_mask is None:
return None
if len(centroid) < 2:
return None
cx, cy = int(centroid[0]), int(centroid[1])
if cy < 0 or cy >= valid_mask.shape[0] or cx < 0 or cx >= valid_mask.shape[1]:
return False
return bool(valid_mask[cy, cx])
def overlay_runtime_status(image, sensor, centroid=None):
"""Annotate a panel with mask-aware runtime diagnostics."""
out = image.copy()
stats = getattr(sensor, "mask_stats", {}) or {}
coverage = 100.0 * float(stats.get("coverage", 0.0))
flow_valid = float(stats.get("flow_mean_valid", 0.0))
flow_invalid = float(stats.get("flow_mean_invalid", 0.0))
centroid_ok = contact_inside_mask(centroid, getattr(sensor, "valid_mask", None))
centroid_label = "n/a" if centroid_ok is None else ("inside" if centroid_ok else "outside")
lines = [
f"mask {coverage:.1f}%",
f"flow valid {flow_valid:.3f}",
f"flow invalid {flow_invalid:.3f}",
f"contact {centroid_label}",
]
y = 54
for line in lines:
cv2.putText(out, line, (12, y), cv2.FONT_HERSHEY_SIMPLEX, 0.55, (255, 255, 255), 2, cv2.LINE_AA)
y += 24
return out
def make_tiled_dashboard(panels, tile_width=640, tile_height=360):
"""Arrange four panels into a 2x2 dashboard."""
tiles = [
add_panel_title(resize_to_fit(image, tile_width, tile_height), title)
for title, image in panels
]
top_row = np.hstack(tiles[:2])
bottom_row = np.hstack(tiles[2:4])
return np.vstack([top_row, bottom_row])
def make_uv_debug_image(data, out_width, out_height, flip_v=True, edge_band=0.12, edge_expand=1.0):
"""Create a UV-domain reference image to verify orientation of the remap."""
face_uvs = face_uvs_from_loop_data(data["faces"], data["uv_data"], data["loop_to_vertex"])
uv_img = np.zeros((out_height, out_width, 3), dtype=np.uint8)
label_color = (255, 255, 255)
for uv_tri in face_uvs:
tri = uv_triangle_to_image(
uv_tri,
out_width,
out_height,
flip_v=flip_v,
edge_band=edge_band,
edge_expand=edge_expand,
)
pts = np.round(tri).astype(np.int32)
mask = np.zeros((out_height, out_width), dtype=np.uint8)
cv2.fillConvexPoly(mask, pts, 255)
x0 = max(0, int(np.floor(np.min(tri[:, 0]))))
x1 = min(out_width - 1, int(np.ceil(np.max(tri[:, 0]))))
y0 = max(0, int(np.floor(np.min(tri[:, 1]))))
y1 = min(out_height - 1, int(np.ceil(np.max(tri[:, 1]))))
if x1 < x0 or y1 < y0:
continue
yy, xx = np.mgrid[y0:y1 + 1, x0:x1 + 1]
inside = mask[y0:y1 + 1, x0:x1 + 1] > 0
if not np.any(inside):
continue
bary = barycentric_grid(xx + 0.5, yy + 0.5, tri)
if bary is None:
continue
w0, w1, w2 = bary
valid = inside & (w0 >= -1e-6) & (w1 >= -1e-6) & (w2 >= -1e-6)
if not np.any(valid):
continue
uv = (
w0[valid, None] * uv_tri[0]
+ w1[valid, None] * uv_tri[1]
+ w2[valid, None] * uv_tri[2]
)
colors = np.zeros((uv.shape[0], 3), dtype=np.uint8)
colors[:, 0] = np.clip(uv[:, 0] * 255.0, 0, 255).astype(np.uint8)
colors[:, 1] = np.clip(uv[:, 1] * 255.0, 0, 255).astype(np.uint8)
colors[:, 2] = np.clip((1.0 - 0.5 * (uv[:, 0] + uv[:, 1])) * 255.0, 0, 255).astype(np.uint8)
uv_img[yy[valid], xx[valid]] = colors
cv2.polylines(uv_img, [pts.reshape(-1, 1, 2)], isClosed=True, color=(255, 255, 255), thickness=1)
cv2.arrowedLine(uv_img, (20, out_height - 20), (min(out_width - 20, 140), out_height - 20), (0, 0, 255), 2, tipLength=0.15)
cv2.arrowedLine(uv_img, (20, out_height - 20), (20, max(20, out_height - 140)), (0, 255, 0), 2, tipLength=0.15)
cv2.putText(uv_img, "U", (min(out_width - 30, 145), out_height - 28), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 0, 255), 2, cv2.LINE_AA)
cv2.putText(uv_img, "V", (28, max(30, out_height - 145)), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0, 255, 0), 2, cv2.LINE_AA)
cv2.putText(uv_img, "UV Map", (20, 35), cv2.FONT_HERSHEY_SIMPLEX, 0.9, label_color, 2, cv2.LINE_AA)
if edge_expand > 1.0:
cv2.putText(uv_img, f"Edge expand {edge_expand:.2f}x", (20, 64), cv2.FONT_HERSHEY_SIMPLEX, 0.55, label_color, 2, cv2.LINE_AA)
return uv_img
def make_uv_distortion_image(data, distortion, out_width, out_height, flip_v=True, edge_band=0.12, edge_expand=1.0):
"""Render UV distortion heatmap; red means more 3D surface per UV area."""
face_uvs = face_uvs_from_loop_data(data["faces"], data["uv_data"], data["loop_to_vertex"])
uv_img = np.zeros((out_height, out_width, 3), dtype=np.uint8)
rel = distortion["relative_scale"]
edge_mask = distortion["edge_mask"]
valid_rel = np.isfinite(rel) & (rel > 0)
if not np.any(valid_rel):
return uv_img
log_rel = np.log2(rel[valid_rel])
clip = max(0.25, float(np.percentile(np.abs(log_rel), 95)))
for idx, uv_tri in enumerate(face_uvs):
tri = uv_triangle_to_image(
uv_tri,
out_width,
out_height,
flip_v=flip_v,
edge_band=edge_band,
edge_expand=edge_expand,
)
pts = np.round(tri).astype(np.int32)
if not np.isfinite(rel[idx]) or rel[idx] <= 0:
color = (64, 64, 64)
else:
value = np.clip((np.log2(rel[idx]) + clip) / (2.0 * clip), 0.0, 1.0)
color_idx = np.uint8(round(value * 255.0))
color = tuple(int(c) for c in cv2.applyColorMap(np.array([[color_idx]], dtype=np.uint8), cv2.COLORMAP_JET)[0, 0])
cv2.fillConvexPoly(uv_img, pts, color)
edge_color = (255, 255, 255) if edge_mask[idx] else (180, 180, 180)
edge_thickness = 2 if edge_mask[idx] else 1
cv2.polylines(uv_img, [pts.reshape(-1, 1, 2)], isClosed=True, color=edge_color, thickness=edge_thickness)
cv2.putText(uv_img, "UV Distortion", (20, 35), cv2.FONT_HERSHEY_SIMPLEX, 0.9, (255, 255, 255), 2, cv2.LINE_AA)
cv2.putText(uv_img, "Blue: denser UV Red: sparser UV", (20, 64), cv2.FONT_HERSHEY_SIMPLEX, 0.55, (220, 220, 220), 2, cv2.LINE_AA)
cv2.putText(uv_img, "White outlines: edge-band faces", (20, 92), cv2.FONT_HERSHEY_SIMPLEX, 0.55, (220, 220, 220), 2, cv2.LINE_AA)
middle_stats = distortion["middle_stats"]
edge_stats = distortion["edge_stats"]
cv2.putText(
uv_img,
f"middle p05/p95: {middle_stats['relative_scale_p05']:.2f}x / {middle_stats['relative_scale_p95']:.2f}x",
(20, 120),
cv2.FONT_HERSHEY_SIMPLEX,
0.55,
(220, 220, 220),
2,
cv2.LINE_AA,
)
cv2.putText(
uv_img,
f"edge p05/p95: {edge_stats['relative_scale_p05']:.2f}x / {edge_stats['relative_scale_p95']:.2f}x",
(20, 148),
cv2.FONT_HERSHEY_SIMPLEX,
0.55,
(220, 220, 220),
2,
cv2.LINE_AA,
)
return uv_img
def open_source(video_path, camera_id):
return cv2.VideoCapture(video_path) if video_path else cv2.VideoCapture(camera_id)
PERF_TUNING_EPILOG = """
Performance tuning (config JSON + CLI):
camera.resolution_stitch [280, 300] UV output + SDK flow/HHD grid
camera.source_scale 0.4 shrink source before uv_prep/remap
tracker.dis_preset medium -> fast faster DIS optical flow
tracker.hhd_downsample 5 -> 4 smaller HHD Poisson grid
contact.denoise_sigma 12 -> 6 or 0 faster FoE contact step
--source_scale 0.4 override config source scale
--masked_poisson_solver sparse vs capfft benchmark HHD solve time
--no-viz skip dashboard ~20 ms/frame savings
Recommended dd02-ov fast profile:
resolution_stitch [280, 300], source_scale 0.4, hhd_downsample 5
Delete or rename stale unroll_map_uv.npz after changing output scale.
"""
def run_deformation_timed(sensor, timer, check_motion=True, threshold=0, tail_frames=10):
"""Run flow + gated HHD with separate benchmark buckets."""
timer.start("flow")
sensor.compute_flow()
timer.stop("flow")
timer.start("hhd")
run_hhd = True
if check_motion:
is_motion = sensor.detect_motion(threshold=threshold)
try:
tail_frames = max(0, int(tail_frames))
except (TypeError, ValueError):
tail_frames = 0
if is_motion:
sensor.tail_frames_count = 0
elif sensor.tail_frames_count < tail_frames:
sensor.tail_frames_count += 1
else:
run_hhd = False
sensor.tail_frames_count = tail_frames
if run_hhd:
sensor.compute_hhd()
timer.stop("hhd")
def configure_masked_poisson(sensor, args):
"""Apply masked-HHD solver settings to the externally fed SDK sensor."""
sensor.masked_poisson_solver = args.masked_poisson_solver
sensor.masked_poisson_max_boundary = args.masked_poisson_max_boundary
sensor.masked_poisson_boundary_stride = args.masked_poisson_boundary_stride
sensor.masked_poisson_max_iters = args.masked_poisson_max_iters
sensor.masked_poisson_tol = args.masked_poisson_tol
sensor.masked_poisson_reg = args.masked_poisson_reg
def resolve_runtime_geometry(args):
"""Apply config defaults and optional source-scale to UV/camera geometry."""
sdk_config = load_config(args.sdk_config)
camera_config = sdk_config.get("camera", {})
out_width, out_height, source_scale = resolve_uv_runtime_sizes(
camera_config,
out_width=args.out_width,
out_height=args.out_height,
source_scale=args.source_scale,
)
stitch = camera_config.get("resolution_stitch")
if stitch and len(stitch) >= 2 and (int(stitch[0]) != out_width or int(stitch[1]) != out_height):
print(
"Warning: UV output "
f"{out_width}x{out_height} differs from resolution_stitch "
f"{int(stitch[0])}x{int(stitch[1])} in {args.sdk_config}; "
"Sensor.get_img will resize again."
)
args.out_width = out_width
args.out_height = out_height
args.source_scale = source_scale
args.map_cache = scaled_map_cache_path(args.map_cache, source_scale)
return sdk_config, camera_config
def main():
parser = argparse.ArgumentParser(
description="Real-time UV-based finger unrolling",
epilog=PERF_TUNING_EPILOG,
formatter_class=argparse.RawDescriptionHelpFormatter,
)
parser.add_argument("--mesh", required=True, help="Exported mesh.npz path")
parser.add_argument("--K", dest="intrinsics_path", required=True, help="Exported camera_intrinsics.npz path")
parser.add_argument("--pose", required=True, help="Exported world_to_cam.npy path")
parser.add_argument("--calib", default=None, help="OpenCV calibration JSON for the live camera")
# With a 12% UV edge band on each side, only about 76% of each axis is the
# lower-distortion middle region. Use a slightly larger full remap so the
# useful center keeps roughly the same sampling density as the prior 506x550.
parser.add_argument("--out_width", type=int, default=None, help="UV flatten width (default: camera.resolution_stitch in sdk config)")
parser.add_argument("--out_height", type=int, default=None, help="UV flatten height (default: camera.resolution_stitch in sdk config)")
parser.add_argument(
"--source_scale",
type=float,
default=None,
help="Resize live source to this fraction before uv_prep/uv_remap; scales K accordingly (default: camera.source_scale or work4[3] in sdk config, else 1.0)",
)
parser.add_argument("--camera_id", type=int, default=0)
parser.add_argument("--video", help="Video file path (overrides camera)")
parser.add_argument("--map_cache", default="unroll_map_uv.npz")
parser.add_argument("--no_flip_v", action="store_true", help="Keep UV V increasing downward")
parser.add_argument("--no_flip_source_vertical", action="store_true", help="Do not vertically flip the live source before UV remap")
parser.add_argument("--edge_band", type=float, default=0.12, help="Normalized UV band near the boundary treated as high-distortion edge")
parser.add_argument("--edge_expand", type=float, default=0.0, help="Extra UV allocation for edge band; 0 chooses automatically from distortion stats")
parser.add_argument("--sdk_config", default="./config/dd02-ov.json", help="Orisys SDK config path/name")
parser.add_argument("--sdk_source", default=None, help="Dummy source used to initialize Sensor; defaults to --video or --camera_id")
parser.add_argument("--masked_poisson_solver", choices=["sparse", "capfft"], default="sparse", help="Masked HHD Poisson solver used inside the SDK")
parser.add_argument("--masked_poisson_max_boundary", type=int, default=96, help="Maximum number of UV-boundary samples used by the capfft solver")
parser.add_argument("--masked_poisson_boundary_stride", type=int, default=1, help="Subsample stride applied before capfft boundary sampling")
parser.add_argument("--masked_poisson_max_iters", type=int, default=3, help="Maximum capfft correction iterations per Poisson solve")
parser.add_argument("--masked_poisson_tol", type=float, default=1e-4, help="Residual tolerance used by the capfft correction loop")
parser.add_argument("--masked_poisson_reg", type=float, default=1e-8, help="Diagonal regularization applied to the capacitance matrix")
parser.add_argument("--verbose", action="store_true", help="Enable verbose SDK logs")
parser.add_argument(
"--no-viz",
action="store_true",
help="Skip OpenCV dashboard; print forces and benchmark timings only",
)
args = parser.parse_args()
_, camera_config = resolve_runtime_geometry(args)
print("[benchmark] Startup timing:")
data = benchmark_step("load_exported_data", load_exported_data, args.mesh, args.intrinsics_path, args.pose)
flip_v = not args.no_flip_v
flip_source_vertical = not args.no_flip_source_vertical
distortion = benchmark_step(
"compute_uv_distortion_metrics",
compute_uv_distortion_metrics,
data,
edge_band=args.edge_band,
)
distortion_stats = distortion["stats"]
middle_stats = distortion["middle_stats"]
edge_stats = distortion["edge_stats"]
edge_expand = args.edge_expand if args.edge_expand > 1.0 else recommend_edge_expand(distortion)
print(
"UV distortion stats: "
f"median mm/uv={distortion_stats['median_mm_per_uv']:.3f}, "
f"relative scale p05/p95={distortion_stats['relative_scale_p05']:.2f}x/"
f"{distortion_stats['relative_scale_p95']:.2f}x, "
f"anisotropy p50/p95={distortion_stats['anisotropy_p50']:.2f}x/"
f"{distortion_stats['anisotropy_p95']:.2f}x"
)
print(f"Using UV edge warp: band={args.edge_band:.2f}, expand={edge_expand:.2f}x")
print(
"UV middle-region stats: "
f"faces={middle_stats['face_count']}, "
f"relative scale p05/p95={middle_stats['relative_scale_p05']:.2f}x/"
f"{middle_stats['relative_scale_p95']:.2f}x, "
f"anisotropy p50/p95={middle_stats['anisotropy_p50']:.2f}x/"
f"{middle_stats['anisotropy_p95']:.2f}x"
)
print(
"UV edge-band stats: "
f"faces={edge_stats['face_count']}, "
f"relative scale p05/p95={edge_stats['relative_scale_p05']:.2f}x/"
f"{edge_stats['relative_scale_p95']:.2f}x, "
f"anisotropy p50/p95={edge_stats['anisotropy_p50']:.2f}x/"
f"{edge_stats['anisotropy_p95']:.2f}x"
)
if args.calib:
source_K, dist_coeffs, img_size = benchmark_step("load_calib_json", load_calib_json, args.calib)
native_width, native_height = img_size
print(f"Calibration loaded. Live image size: {(native_width, native_height)}")
else:
source_K = data["K"]
dist_coeffs = None
native_width = data["width"]
native_height = data["height"]
print("No calibration file provided; using exported intrinsics for projection.")
source_K, source_width, source_height = scale_camera_geometry(
source_K,
native_width,
native_height,
args.source_scale,
)
print(
"UV runtime geometry: "
f"native={native_width}x{native_height}, "
f"source={source_width}x{source_height} (scale={args.source_scale:g}), "
f"output={args.out_width}x{args.out_height}, "
f"map_cache={args.map_cache}"
)
domain = benchmark_step(
"build_or_load_domain",
build_or_load_domain,
args.mesh,
args.intrinsics_path,
args.pose,
args.map_cache,
source_K,
source_width,
source_height,
args.out_width,
args.out_height,
flip_v=flip_v,
edge_band=args.edge_band,
edge_expand=edge_expand,
include_surface=True,
source_scale=args.source_scale,
)
remap = domain.remap
remap.ensure_fast_maps()
edge_expand = float(remap.metadata.get("edge_expand", edge_expand))
uv_debug = benchmark_step(
"make_uv_debug_image",
make_uv_debug_image,
data,
out_width=args.out_width,
out_height=args.out_height,
flip_v=flip_v,
edge_band=args.edge_band,
edge_expand=edge_expand,
)
uv_distortion = benchmark_step(
"make_uv_distortion_image",
make_uv_distortion_image,
data,
distortion,
out_width=args.out_width,
out_height=args.out_height,
flip_v=flip_v,
edge_band=args.edge_band,
edge_expand=edge_expand,
)
uv_reference = benchmark_step(
"uv_reference_compose",
lambda: np.hstack(
[
resize_to_height(uv_debug, args.out_height),
resize_to_height(uv_distortion, args.out_height),
]
),
)
cap = open_source(args.video, args.camera_id)
if not cap.isOpened():
print("Cannot open camera/video source")
sys.exit(1)
cap.set(cv2.CAP_PROP_FRAME_WIDTH, native_width)
cap.set(cv2.CAP_PROP_FRAME_HEIGHT, native_height)
sensor = benchmark_step(
"sensor_init",
lambda: orisys.Sensor(
"external",
config_name=args.sdk_config,
verbose=args.verbose,
mirror=False,
),
)
configure_masked_poisson(sensor, args)
print(
"Masked HHD solver: "
f"{sensor.masked_poisson_solver} "
f"(max_boundary={sensor.masked_poisson_max_boundary}, "
f"stride={sensor.masked_poisson_boundary_stride}, "
f"iters={sensor.masked_poisson_max_iters}, "
f"tol={sensor.masked_poisson_tol:.1e})"
)
frame_processor = FrameUnwarper(
domain.remap,
camera_matrix=source_K if dist_coeffs is not None else None,
dist_coeffs=dist_coeffs,
source_width=source_width,
source_height=source_height,
flip_source_vertical=flip_source_vertical,
border_value=(0, 0, 0),
)
ok, warmup_frame = cap.read()
if ok:
warmup_prepared = frame_processor.prepare_source(warmup_frame)
warmup_unrolled = frame_processor.flatten_prepared(warmup_prepared)
sensor.get_img(
frame=warmup_unrolled,
already_processed=True,
valid_mask=frame_processor.valid_mask,
)
benchmark_step("hhd_warmup", sensor.warmup_deformation)
else:
print("[benchmark] hhd_warmup: skipped (no frame from source)")
if args.no_viz:
print("Running with --no-viz (dashboard disabled).")
else:
print("Press 'q' to quit.")
runtime_timer = StepTimer()
try:
while True:
runtime_timer.start("capture")
ok, frame = cap.read()
runtime_timer.stop("capture")
if not ok:
break
runtime_timer.start("uv_prepare")
frame_source = frame_processor.prepare_source(frame)
runtime_timer.stop("uv_prepare")
runtime_timer.start("uv_remap")
unrolled = frame_processor.flatten_prepared(frame_source)
runtime_timer.stop("uv_remap")
runtime_timer.start("get_img")
sensor.get_img(
frame=unrolled,
already_processed=True,
valid_mask=frame_processor.valid_mask,
)
runtime_timer.stop("get_img")
run_deformation_timed(sensor, runtime_timer, check_motion=True, threshold=0)
runtime_timer.start("contact")
sensor.compute_contact()
runtime_timer.stop("contact")
fps = sensor._fps
fnormal = sensor.fnormal
fshearx = sensor.fshearx
fsheary = sensor.fsheary
contact_ok = contact_inside_mask(sensor.centroid, sensor.valid_mask)
if not args.no_viz:
runtime_timer.start("viz")
fps, flow, vnormal, vshear, img, contour, centroid, depth_map, fnormal, fshearx, fsheary = sensor.read_info(
sensor.info.FPS,
sensor.info.VRAW,
sensor.info.VNORMAL,
sensor.info.VSHEAR,
sensor.info.IMG,
sensor.info.CONTOUR,
sensor.info.CENTROID,
sensor.info.DEPTH,
sensor.info.FNORMAL,
sensor.info.FSHEARX,
sensor.info.FSHEARY,
)
arrows = orisys.util.draw_arrows(
img,
flow,
threshold=0.4,
grid_spacing=10,
arrow_scale=1.0,
)
contact_img = orisys.util.draw_contact(img, contour, centroid)
depth_norm = depth_map
if depth_norm.size > 0 and depth_norm.max() > depth_norm.min():
depth_u8 = ((depth_norm - depth_norm.min()) / (depth_norm.max() - depth_norm.min()) * 255).astype(np.uint8)
else:
depth_u8 = np.zeros_like(depth_norm, dtype=np.uint8)
depth_colored = cv2.applyColorMap(depth_u8, cv2.COLORMAP_JET)
display_mask = sensor.valid_mask
arrows = apply_display_mask(arrows, display_mask, border_value=(0, 0, 0))
contact_img = apply_display_mask(contact_img, display_mask, border_value=(0, 0, 0))
depth_colored = apply_display_mask(depth_colored, display_mask, border_value=(0, 0, 0))
arrows = overlay_runtime_status(arrows, sensor, centroid=centroid)
contact_ok = contact_inside_mask(centroid, display_mask)
dashboard = make_tiled_dashboard(
[
("Original", frame_source),
("UV Map | Distortion", uv_reference),
("UV Unrolled", unrolled),
("Flow", arrows),
],
tile_width=640,
tile_height=360,
)
cv2.imshow("UV Unroll Dashboard", dashboard)
runtime_timer.stop("viz")
print(
f"FPS={fps:.2f}, 法向力={fnormal:.4f}, 切向力X={fshearx:.4f}, "
f"切向力Y={fsheary:.4f}, 展开图尺寸={unrolled.shape[:2]}, "
f"solver={sensor.masked_poisson_solver}, "
f"mask={sensor.mask_stats['coverage']:.3f}, "
f"flow(valid/invalid)=({sensor.mask_stats['flow_mean_valid']:.3f}/"
f"{sensor.mask_stats['flow_mean_invalid']:.3f}), "
f"contact_in_mask={contact_ok}, "
f"{print_runtime_benchmark(runtime_timer)}",
end="\r",
)
key = cv2.waitKey(1) & 0xFF
if key == ord("q"):
print("\n正在退出程序...")
break
if key == ord("r"):
sensor.reset()
print("\n追踪器已重置。")
finally:
cap.release()
sensor.disconnect()
cv2.destroyAllWindows()
if __name__ == "__main__":
main()
Binary file not shown.

Before

Width:  |  Height:  |  Size: 934 KiB

Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
-188
View File
@@ -1,188 +0,0 @@
import numpy as np
import cv2
import os
def load_exported_data(data_dir):
"""Load all exported files."""
# Mesh
mesh = np.load(os.path.join(data_dir, "mesh.npz"))
vertices = mesh['vertices']
faces = mesh['faces']
uv_data = mesh['uv_data']
loop_to_vertex = mesh['loop_to_vertex']
# Camera intrinsics
cam_data = np.load(os.path.join(data_dir, "camera_intrinsics.npz"))
K = cam_data['K']
width = int(cam_data['width'])
height = int(cam_data['height'])
# Camera pose (world-to-camera)
world_to_cam = np.load(os.path.join(data_dir, "world_to_cam.npy"))
# Precomputed pixel-to-3D map (optional)
xyz_map_path = os.path.join(data_dir, "pixel_to_xyz.npz")
if os.path.exists(xyz_map_path):
xyz_map = np.load(xyz_map_path)['xyz_map']
else:
xyz_map = None
return {
'vertices': vertices,
'faces': faces,
'uv_data': uv_data,
'loop_to_vertex': loop_to_vertex,
'K': K,
'width': width,
'height': height,
'world_to_cam': world_to_cam,
'xyz_map': xyz_map
}
def ray_intersect_triangle(ray_origin, ray_dir, v0, v1, v2):
"""Möller–Trumbore intersection algorithm. Returns t (distance) or None."""
epsilon = 1e-8
edge1 = v1 - v0
edge2 = v2 - v0
h = np.cross(ray_dir, edge2)
a = np.dot(edge1, h)
if abs(a) < epsilon:
return None
f = 1.0 / a
s = ray_origin - v0
u = f * np.dot(s, h)
if u < 0.0 or u > 1.0:
return None
q = np.cross(s, edge1)
v = f * np.dot(ray_dir, q)
if v < 0.0 or u + v > 1.0:
return None
t = f * np.dot(edge2, q)
if t > epsilon:
return t
return None
def pixel_to_3d_ray(u, v, K, world_to_cam):
"""
Given pixel coordinates (u,v) and camera intrinsics K,
return a ray in world coordinates: origin and direction.
"""
# Invert K to get normalized image coordinates
fx = K[0,0]; fy = K[1,1]; cx = K[0,2]; cy = K[1,2]
x_norm = (u - cx) / fx
y_norm = (v - cy) / fy
# Ray in camera space: (x_norm, y_norm, 1.0)
ray_cam = np.array([x_norm, y_norm, 1.0], dtype=np.float64)
# Transform to world
cam_to_world = np.linalg.inv(world_to_cam)
R = cam_to_world[:3, :3]
t = cam_to_world[:3, 3]
ray_world_dir = R @ ray_cam
ray_world_dir /= np.linalg.norm(ray_world_dir)
ray_world_origin = t
return ray_world_origin, ray_world_dir
def map_image_to_3d(image_path, data, use_precomputed=True, flip_vertical=True):
"""
Map each pixel of the image to a 3D world point.
If use_precomputed and xyz_map exists, return that directly.
If `flip_vertical` is True, flip the loaded image top-to-bottom so its
row order matches the exported Blender camera/pixel mapping.
Otherwise, perform raycasting on the mesh.
Returns:
points: (H, W, 3) array of 3D points (NaN for misses).
colors: (H, W, 3) array of RGB colors from image.
"""
img = cv2.imread(image_path)
if img is None:
raise FileNotFoundError(f"Could not read image: {image_path}")
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB) # to RGB
if flip_vertical:
img = np.flipud(img).copy()
h, w = img.shape[:2]
if use_precomputed and data['xyz_map'] is not None:
# Use precomputed map (must match resolution)
if h != data['height'] or w != data['width']:
raise ValueError("Image resolution does not match precomputed map.")
points = data['xyz_map'].copy()
colors = img.astype(np.float32) / 255.0
return points, colors
# Otherwise, perform raycasting on the mesh
vertices = data['vertices']
faces = data['faces']
K = data['K']
world_to_cam = data['world_to_cam']
points = np.full((h, w, 3), np.nan, dtype=np.float32)
colors = img.astype(np.float32) / 255.0
# For speed, we can use a bounding volume hierarchy (BVH).
# Here we implement a simple brute-force over all faces (slow for large meshes).
# In practice, precomputing the xyz_map is strongly recommended.
for y in range(h):
for x in range(w):
origin, direction = pixel_to_3d_ray(x, y, K, world_to_cam)
best_t = np.inf
best_hit = None
for face in faces:
v0 = vertices[face[0]]
v1 = vertices[face[1]]
v2 = vertices[face[2]]
t = ray_intersect_triangle(origin, direction, v0, v1, v2)
if t is not None and t < best_t:
best_t = t
best_hit = origin + direction * t
if best_hit is not None:
points[y, x] = best_hit
return points, colors
def save_colored_pointcloud(points, colors, output_path):
"""Save as PLY file (colored point cloud)."""
# Flatten valid points
mask = ~np.isnan(points[:,:,0])
valid_points = points[mask].reshape(-1, 3)
valid_colors = (colors[mask].reshape(-1, 3) * 255).astype(np.uint8)
with open(output_path, 'w') as f:
f.write("ply\n")
f.write("format ascii 1.0\n")
f.write(f"element vertex {len(valid_points)}\n")
f.write("property float x\n")
f.write("property float y\n")
f.write("property float z\n")
f.write("property uchar red\n")
f.write("property uchar green\n")
f.write("property uchar blue\n")
f.write("end_header\n")
for pt, col in zip(valid_points, valid_colors):
f.write(f"{pt[0]} {pt[1]} {pt[2]} {col[0]} {col[1]} {col[2]}\n")
print(f"Saved point cloud to {output_path}")
if __name__ == "__main__":
data_dir = "./exported_data" # Path to exported files
image_path = "captured_image_1.jpg" # Input camera image
output_ply = "mapped_points.ply"
data = load_exported_data(data_dir)
points, colors = map_image_to_3d(
image_path,
data,
use_precomputed=True,
flip_vertical=True,
)
save_colored_pointcloud(points, colors, output_ply)
# Optional: Visualize using Open3D
import open3d as o3d
pcd = o3d.geometry.PointCloud()
mask = ~np.isnan(points[:,:,0])
pcd.points = o3d.utility.Vector3dVector(points[mask].reshape(-1,3))
pcd.colors = o3d.utility.Vector3dVector(colors[mask].reshape(-1,3))
o3d.visualization.draw_geometries([pcd])
-629
View File
@@ -1,629 +0,0 @@
import argparse
import os
import sys
import time
import cv2
import numpy as np
# UV-unroll mesh viewer (grey-white surface, UV-space field mapping):
# from examples.viz.mesh_viewer import UVMeshViewer, build_uv_surface_mesh
# # or: from orisys.finger import UVMeshViewer, build_uv_surface_mesh
# See examples/finger_pipeline_3dview.py
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
REPO_ROOT = os.path.abspath(os.path.join(SCRIPT_DIR, "..", ".."))
SRC_DIR = os.path.join(REPO_ROOT, "src")
if SRC_DIR not in sys.path:
sys.path.insert(0, SRC_DIR)
import orisys
def load_exported_data(data_dir):
mesh = np.load(os.path.join(data_dir, "mesh.npz"))
cam_data = np.load(os.path.join(data_dir, "camera_intrinsics.npz"))
xyz_map_path = os.path.join(data_dir, "pixel_to_xyz.npz")
if not os.path.exists(xyz_map_path):
raise FileNotFoundError(f"Missing precomputed 3D map: {xyz_map_path}")
xyz_map = np.load(xyz_map_path)["xyz_map"]
width = int(cam_data["width"])
height = int(cam_data["height"])
if xyz_map.shape[:2] != (height, width):
raise ValueError(
f"pixel_to_xyz shape {xyz_map.shape[:2]} does not match camera size {(height, width)}"
)
return {
"vertices": mesh["vertices"],
"faces": mesh["faces"],
"uv_data": mesh["uv_data"],
"loop_to_vertex": mesh["loop_to_vertex"],
"K": cam_data["K"],
"width": width,
"height": height,
"world_to_cam": np.load(os.path.join(data_dir, "world_to_cam.npy")),
"xyz_map": xyz_map,
}
def parse_source(value):
value = str(value).strip()
for prefix in ("camera:", "cam:", "webcam:"):
if value.lower().startswith(prefix):
value = value[len(prefix):].strip()
break
try:
return int(value)
except ValueError:
return value
def prepare_static_points(data, display_stride):
xyz_map = data["xyz_map"]
height, width = xyz_map.shape[:2]
valid_mask = ~np.isnan(xyz_map[:, :, 0])
valid_points = xyz_map[valid_mask].reshape(-1, 3)
if len(valid_points) == 0:
raise ValueError("pixel_to_xyz.npz contains no valid mesh hit points.")
display_indices = np.arange(len(valid_points))[::max(1, int(display_stride))]
flat_indices = np.flatnonzero(valid_mask.ravel())[display_indices]
display_ys = flat_indices // width
display_xs = flat_indices % width
return (
valid_mask,
display_indices,
valid_points[display_indices],
display_ys,
display_xs,
)
def mesh_pixel_to_world_scale(xyz_map, valid_mask):
valid_points = xyz_map[valid_mask]
if len(valid_points) == 0:
return 1e-3
extent = np.nanmax(valid_points, axis=0) - np.nanmin(valid_points, axis=0)
mesh_span = float(np.max(extent))
height, width = xyz_map.shape[:2]
return mesh_span / max(height, width, 1)
def sample_field_at_display_points(field, display_ys, display_xs):
if field.ndim == 2:
return field[display_ys, display_xs]
return field[display_ys, display_xs, :]
def deform_point_cloud(rest_points, displacement, pixel_to_world, deform_scale):
return rest_points + displacement * pixel_to_world * float(deform_scale)
def lift_flow_to_3d(flow, tangent_x, tangent_y, normals):
"""Lift raw 2D optical flow to tangent-plane 3D displacement."""
tx_norm = np.linalg.norm(tangent_x, axis=2, keepdims=True)
ty_norm = np.linalg.norm(tangent_y, axis=2, keepdims=True)
t1 = np.divide(tangent_x, tx_norm, out=np.zeros_like(tangent_x), where=tx_norm > 1e-8)
t2 = np.divide(tangent_y, ty_norm, out=np.zeros_like(tangent_y), where=ty_norm > 1e-8)
flow_3d = flow[:, :, 0:1] * t1 + flow[:, :, 1:2] * t2
normal_component = np.sum(flow_3d * normals, axis=2, keepdims=True)
return flow_3d - normal_component * normals
def align_image_to_map(img_bgr, data, valid_mask, display_indices, flip_vertical=True):
target_w = int(data["width"])
target_h = int(data["height"])
h, w = img_bgr.shape[:2]
if (w, h) != (target_w, target_h):
img_bgr = cv2.resize(img_bgr, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
img_rgb = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB)
if flip_vertical:
img_rgb = np.flipud(img_rgb).copy()
colors = img_rgb.astype(np.float64) / 255.0
return colors[valid_mask].reshape(-1, 3)[display_indices]
def arrow_map_reference_colors(count, color=(0.04, 0.04, 0.04)):
return np.tile(np.asarray(color, dtype=np.float64), (count, 1))
def magnitude_map_colors(flow, valid_mask, display_indices, max_magnitude=None):
mag = np.linalg.norm(flow, axis=2)
valid_mag = mag[valid_mask].reshape(-1)
display_mag = valid_mag[display_indices]
if max_magnitude is None or max_magnitude <= 0:
max_magnitude = max(float(np.percentile(display_mag, 95)), 1e-6)
value = np.clip(display_mag / max_magnitude, 0.0, 1.0)
return np.column_stack(
(
value,
np.clip(1.0 - np.abs(value - 0.5) * 2.0, 0.0, 1.0),
1.0 - value,
)
).astype(np.float64)
def align_scalar_field_to_map(field, data, flip_vertical=True):
if field.ndim == 3:
field = np.linalg.norm(field, axis=2)
target_w = int(data["width"])
target_h = int(data["height"])
h, w = field.shape[:2]
if (w, h) != (target_w, target_h):
field = cv2.resize(field, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
else:
field = field.copy()
if flip_vertical:
field = np.flipud(field).copy()
return field
def scalar_map_colors(field, valid_mask, display_indices, value_min=None, value_max=None):
valid_values = field[valid_mask].reshape(-1)
display_values = valid_values[display_indices]
finite_values = display_values[np.isfinite(display_values)]
if len(finite_values) == 0:
return arrow_map_reference_colors(len(display_indices))
if value_min is None:
value_min = float(np.percentile(finite_values, 5))
if value_max is None:
value_max = float(np.percentile(finite_values, 95))
if value_max <= value_min:
value_max = value_min + 1e-6
value = np.clip((display_values - value_min) / (value_max - value_min), 0.0, 1.0)
return np.column_stack(
(
value,
np.clip(1.0 - np.abs(value - 0.5) * 2.0, 0.0, 1.0),
1.0 - value,
)
).astype(np.float64)
def align_vector_field_to_map(field, data, flip_vertical=True, scale_as_flow=True):
target_w = int(data["width"])
target_h = int(data["height"])
h, w = field.shape[:2]
if (w, h) != (target_w, target_h):
sx = target_w / max(w, 1)
sy = target_h / max(h, 1)
field = cv2.resize(field, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
if scale_as_flow:
field[:, :, 0] *= sx
field[:, :, 1] *= sy
else:
field = field.copy()
if flip_vertical:
field = np.flipud(field).copy()
field[:, :, 1] *= -1.0
return field
def estimate_surface_normals(xyz_map):
valid = ~np.isnan(xyz_map[:, :, 0])
filled = xyz_map.copy()
filled[~valid] = 0
d_y = np.gradient(filled, axis=0)
d_x = np.gradient(filled, axis=1)
normals = np.cross(d_x, d_y)
norm = np.linalg.norm(normals, axis=2, keepdims=True)
normals = np.divide(normals, norm, out=np.zeros_like(normals), where=norm > 1e-8)
normals[~valid] = 0
return normals
def compute_surface_frames(xyz_map):
valid = ~np.isnan(xyz_map[:, :, 0])
filled = xyz_map.copy()
filled[~valid] = 0
d_y = np.gradient(filled, axis=0)
d_x = np.gradient(filled, axis=1)
normals = np.cross(d_x, d_y)
norm = np.linalg.norm(normals, axis=2, keepdims=True)
normals = np.divide(normals, norm, out=np.zeros_like(normals), where=norm > 1e-8)
d_x[~valid] = 0
d_y[~valid] = 0
normals[~valid] = 0
return d_x, d_y, normals, valid
def lift_deformation_to_3d(vshear, vnormal, tangent_x, tangent_y, normals, normal_gain=1.0):
# Unit tangents: map 2D pixel shear to 3D direction without world-scale Jacobian blow-up.
tx_norm = np.linalg.norm(tangent_x, axis=2, keepdims=True)
ty_norm = np.linalg.norm(tangent_y, axis=2, keepdims=True)
t1 = np.divide(tangent_x, tx_norm, out=np.zeros_like(tangent_x), where=tx_norm > 1e-8)
t2 = np.divide(tangent_y, ty_norm, out=np.zeros_like(tangent_y), where=ty_norm > 1e-8)
shear_3d = vshear[:, :, 0:1] * t1 + vshear[:, :, 1:2] * t2
normal_mag = np.linalg.norm(vnormal, axis=2, keepdims=True) * float(normal_gain)
normal_3d = normal_mag * normals
normal_component = np.sum(shear_3d * normals, axis=2, keepdims=True)
shear_tangent = shear_3d - normal_component * normals
return shear_tangent + normal_3d
def line_colors(mag, arrow_map_style):
if len(mag) == 0:
return np.empty((0, 3), dtype=np.float64)
if arrow_map_style:
return np.tile(np.array([0.0, 0.9, 1.0], dtype=np.float64), (len(mag), 1))
norm_mag = mag / max(float(np.percentile(mag, 95)), 1e-6)
norm_mag = np.clip(norm_mag, 0.0, 1.0)
return np.column_stack((norm_mag, 1.0 - norm_mag, np.zeros_like(norm_mag))).astype(np.float64)
def lines_from_vectors(starts, vectors, mag, flow_scale, arrow_map_style, unit_direction=False):
vec_norm = np.linalg.norm(vectors, axis=1)
good = np.isfinite(vec_norm) & (vec_norm > 1e-8)
starts, vectors, mag = starts[good], vectors[good], mag[good]
if len(starts) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
if unit_direction:
vec_norm = vec_norm[:, None]
vectors = np.divide(vectors, vec_norm, out=np.zeros_like(vectors), where=vec_norm > 1e-8)
points = np.empty((len(starts) * 2, 3), dtype=np.float64)
points[0::2] = starts
points[1::2] = starts + vectors * flow_scale
lines = np.column_stack((np.arange(0, len(points), 2), np.arange(1, len(points), 2))).astype(np.int32)
return points, lines, line_colors(mag, arrow_map_style)
def build_flow_lines(
xyz_map,
valid_mask,
flow,
arrow_stride=32,
flow_scale=1.0,
min_magnitude=0.5,
arrow_map_style=False,
):
h, w = flow.shape[:2]
ys, xs = np.mgrid[0:h:arrow_stride, 0:w:arrow_stride]
ys = ys.reshape(-1)
xs = xs.reshape(-1)
sample_flow = flow[ys, xs]
mag = np.linalg.norm(sample_flow, axis=1)
keep = valid_mask[ys, xs] & (mag >= min_magnitude)
xs, ys, sample_flow, mag = xs[keep], ys[keep], sample_flow[keep], mag[keep]
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
tx = np.rint(xs + sample_flow[:, 0]).astype(np.int32)
ty = np.rint(ys + sample_flow[:, 1]).astype(np.int32)
inside = (0 <= tx) & (tx < w) & (0 <= ty) & (ty < h)
xs, ys, tx, ty, mag = xs[inside], ys[inside], tx[inside], ty[inside], mag[inside]
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
valid_target = valid_mask[ty, tx]
xs, ys, tx, ty, mag = xs[valid_target], ys[valid_target], tx[valid_target], ty[valid_target], mag[valid_target]
starts = xyz_map[ys, xs]
vectors = xyz_map[ty, tx] - starts
return lines_from_vectors(starts, vectors, mag, flow_scale, arrow_map_style)
def build_resultant_force_lines(
xyz_map,
valid_mask,
resultant_force,
arrow_stride=32,
flow_scale=1.0,
min_magnitude=0.01,
arrow_map_style=False,
):
h, w = resultant_force.shape[:2]
ys, xs = np.mgrid[0:h:arrow_stride, 0:w:arrow_stride]
ys = ys.reshape(-1)
xs = xs.reshape(-1)
vectors = resultant_force[ys, xs]
mag = np.linalg.norm(vectors, axis=1)
keep = valid_mask[ys, xs] & (mag >= min_magnitude)
xs, ys, vectors, mag = (
xs[keep],
ys[keep],
vectors[keep],
mag[keep],
)
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
starts = xyz_map[ys, xs]
return lines_from_vectors(
starts, vectors, mag, flow_scale, arrow_map_style, unit_direction=True
)
def init_viewer(
cloud_points,
cloud_colors,
flow_points,
flow_lines,
flow_colors,
arrow_map_style=False,
show_arrows=True,
):
import open3d as o3d
vis = o3d.visualization.Visualizer()
vis.create_window(window_name="SDK 3D Flow", width=1280, height=720)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(cloud_points.astype(np.float64))
pcd.colors = o3d.utility.Vector3dVector(cloud_colors.astype(np.float64))
vis.add_geometry(pcd)
line_set = None
if show_arrows:
line_set = o3d.geometry.LineSet()
line_set.points = o3d.utility.Vector3dVector(flow_points.astype(np.float64))
line_set.lines = o3d.utility.Vector2iVector(flow_lines)
line_set.colors = o3d.utility.Vector3dVector(flow_colors.astype(np.float64))
vis.add_geometry(line_set)
render_option = vis.get_render_option()
render_option.point_size = 1.0 if arrow_map_style else 2.0
render_option.line_width = 2.0
render_option.background_color = np.array([0.0, 0.0, 0.0])
vis.poll_events()
vis.update_renderer()
return vis, pcd, line_set, o3d
def run_sdk_3d_flow(args):
data = load_exported_data(args.data_dir)
valid_mask, display_indices, rest_cloud_points, display_ys, display_xs = prepare_static_points(
data, args.display_stride
)
tangent_x, tangent_y, surface_normals, valid_mask = compute_surface_frames(data["xyz_map"])
pixel_to_world = mesh_pixel_to_world_scale(data["xyz_map"], valid_mask)
show_arrows = args.display_mode == "arrows"
deform_cloud = args.display_mode == "deform"
sensor = orisys.Sensor(
parse_source(args.source),
config_name=args.config,
isstitch=False,
cuda=args.cuda,
verbose=args.verbose,
mirror=not args.no_mirror,
)
vis = pcd = line_set = o3d = None
frame_count = 0
t0 = time.perf_counter()
try:
while True:
if vis is not None and not vis.poll_events():
break
if sensor.get_img() is None:
break
needs_decomposition = args.vector_mode in ("shear-normal", "depth-map")
sensor.compute_deformation(decompose=needs_decomposition, check_motion=False)
if args.vector_mode == "depth-map":
sensor.compute_contact()
flow, img_bgr, vshear, vnormal, depth_map = sensor.read_info(
sensor.info.VRAW,
sensor.info.IMG,
sensor.info.VSHEAR,
sensor.info.VNORMAL,
sensor.info.DEPTH,
)
flow_map = align_vector_field_to_map(flow, data, flip_vertical=not args.no_flip_vertical)
if args.vector_mode == "vraw-magnitude":
cloud_colors = magnitude_map_colors(flow_map, valid_mask, display_indices, args.magnitude_max)
elif args.vector_mode == "depth-map":
depth_map = align_scalar_field_to_map(depth_map, data, flip_vertical=not args.no_flip_vertical)
cloud_colors = scalar_map_colors(
depth_map,
valid_mask,
display_indices,
value_min=args.depth_min,
value_max=args.depth_max,
)
elif args.arrow_map_style:
cloud_colors = arrow_map_reference_colors(len(rest_cloud_points))
else:
cloud_colors = align_image_to_map(
img_bgr,
data,
valid_mask,
display_indices,
flip_vertical=not args.no_flip_vertical,
)
flip_vertical = not args.no_flip_vertical
displacement_map = None
if args.vector_mode == "shear-normal":
vshear_map = align_vector_field_to_map(vshear, data, flip_vertical=flip_vertical)
vnormal_map = align_vector_field_to_map(vnormal, data, flip_vertical=flip_vertical)
displacement_map = lift_deformation_to_3d(
vshear_map,
vnormal_map,
tangent_x,
tangent_y,
surface_normals,
normal_gain=args.normal_gain,
)
elif args.vector_mode == "depth-map" and deform_cloud:
vnormal_map = align_vector_field_to_map(vnormal, data, flip_vertical=flip_vertical)
zero_shear = np.zeros_like(vnormal_map)
displacement_map = lift_deformation_to_3d(
zero_shear,
vnormal_map,
tangent_x,
tangent_y,
surface_normals,
normal_gain=args.normal_gain,
)
elif args.vector_mode in ("vraw", "vraw-magnitude"):
displacement_map = lift_flow_to_3d(flow_map, tangent_x, tangent_y, surface_normals)
if deform_cloud and args.deform_color == "magnitude" and displacement_map is not None:
cloud_colors = magnitude_map_colors(
displacement_map, valid_mask, display_indices, args.magnitude_max
)
if deform_cloud:
if displacement_map is None:
raise ValueError(
f"display-mode deform is not supported with vector-mode {args.vector_mode!r}"
)
point_displacement = sample_field_at_display_points(
displacement_map, display_ys, display_xs
)
cloud_points = deform_point_cloud(
rest_cloud_points,
point_displacement,
pixel_to_world,
args.deform_scale,
)
flow_points = np.empty((0, 3), dtype=np.float64)
flow_lines = np.empty((0, 2), dtype=np.int32)
flow_colors = np.empty((0, 3), dtype=np.float64)
else:
cloud_points = rest_cloud_points
if args.vector_mode == "shear-normal":
flow_points, flow_lines, flow_colors = build_resultant_force_lines(
data["xyz_map"],
valid_mask,
displacement_map,
arrow_stride=args.arrow_stride,
flow_scale=args.flow_scale,
min_magnitude=args.min_flow,
arrow_map_style=args.arrow_map_style,
)
elif args.vector_mode in ("vraw-magnitude", "depth-map"):
flow_points = np.empty((0, 3), dtype=np.float64)
flow_lines = np.empty((0, 2), dtype=np.int32)
flow_colors = np.empty((0, 3), dtype=np.float64)
else:
flow_points, flow_lines, flow_colors = build_flow_lines(
data["xyz_map"],
valid_mask,
flow_map,
arrow_stride=args.arrow_stride,
flow_scale=args.flow_scale,
min_magnitude=args.min_flow,
arrow_map_style=args.arrow_map_style,
)
if vis is None:
vis, pcd, line_set, o3d = init_viewer(
cloud_points,
cloud_colors,
flow_points,
flow_lines,
flow_colors,
arrow_map_style=args.arrow_map_style,
show_arrows=show_arrows,
)
mode_label = "deformed cloud" if deform_cloud else "arrows"
print(f"SDK 3D flow viewer started ({mode_label}). Close the Open3D window to stop.")
else:
pcd.points = o3d.utility.Vector3dVector(cloud_points.astype(np.float64))
pcd.colors = o3d.utility.Vector3dVector(cloud_colors)
if line_set is not None:
line_set.points = o3d.utility.Vector3dVector(flow_points.astype(np.float64))
line_set.lines = o3d.utility.Vector2iVector(flow_lines)
line_set.colors = o3d.utility.Vector3dVector(flow_colors.astype(np.float64))
vis.update_geometry(line_set)
vis.update_geometry(pcd)
vis.update_renderer()
frame_count += 1
if frame_count % 30 == 0:
elapsed = max(time.perf_counter() - t0, 1e-6)
arrow_count = len(flow_lines) if show_arrows else 0
print(f"frames={frame_count}, viewer FPS={frame_count / elapsed:.1f}, arrows={arrow_count}")
finally:
sensor.disconnect()
if vis is not None:
vis.destroy_window()
def parse_args():
parser = argparse.ArgumentParser(description="Map SDK flow fields to a 3D point cloud.")
parser.add_argument("--data-dir", default=os.path.join(SCRIPT_DIR, "exported_data"))
parser.add_argument(
"--source",
default=os.path.join(SCRIPT_DIR, "20260602-094916.mp4"),
help="Video path or live camera ID, for example 0, 1, camera:0, or cam:1.",
)
parser.add_argument("--config", default=os.path.join(REPO_ROOT, "config", "dd02-ov.json"))
parser.add_argument("--display-stride", type=int, default=4)
parser.add_argument("--arrow-stride", type=int, default=8)
parser.add_argument(
"--flow-scale",
type=float,
default=3.0,
help="Arrow length in 3D viewer units (shear-normal uses unit direction; magnitude shown by color).",
)
parser.add_argument(
"--min-flow",
type=float,
default=0.5,
help="Minimum 2D-equivalent resultant magnitude before drawing an arrow.",
)
parser.add_argument("--normal-gain", type=float, default=1.0)
parser.add_argument(
"--display-mode",
choices=["arrows", "deform"],
default="arrows",
help="arrows: line vectors; deform: displace the point cloud by the 3D deformation field.",
)
parser.add_argument(
"--deform-scale",
type=float,
default=1.0,
help="Gain on cloud displacement (multiplied by mesh pixel-to-world scale).",
)
parser.add_argument(
"--deform-color",
choices=["texture", "magnitude"],
default="texture",
help="Point colors in deform mode: camera texture or deformation magnitude heat-map.",
)
parser.add_argument("--magnitude-max", type=float, default=None,
help="Upper VRAW magnitude for heat-map color scaling; default uses frame percentile.")
parser.add_argument("--depth-min", type=float, default=None,
help="Minimum depth value for heat-map color scaling; default uses frame percentile.")
parser.add_argument("--depth-max", type=float, default=None,
help="Maximum depth value for heat-map color scaling; default uses frame percentile.")
parser.add_argument(
"--vector-mode",
choices=["vraw", "vraw-magnitude", "shear-normal", "depth-map"],
default="vraw",
help="Vector display mode. shear-normal shows 3D resultant vectors (surface-tangent shear + surface-normal component).",
)
parser.add_argument(
"--arrow-map-style",
action="store_true",
help="Hide real texture and use a dark Qt-viewer-like arrow map style.",
)
parser.add_argument("--no-flip-vertical", action="store_true")
parser.add_argument("--no-mirror", action="store_true")
parser.add_argument("--cuda", action="store_true")
parser.add_argument("--verbose", action="store_true")
return parser.parse_args()
if __name__ == "__main__":
run_sdk_3d_flow(parse_args())
-564
View File
@@ -1,564 +0,0 @@
import argparse
import os
import sys
import time
import cv2
import numpy as np
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
REPO_ROOT = os.path.abspath(os.path.join(SCRIPT_DIR, "..", ".."))
SRC_DIR = os.path.join(REPO_ROOT, "src")
if SRC_DIR not in sys.path:
sys.path.insert(0, SRC_DIR)
import orisys
from orisys.unroll import UVDomain
def load_flatten_remap(path):
"""Load flatten cache; returns dict compatible with legacy viewer code."""
domain = UVDomain.load(path)
result = {
"map_x": domain.remap.map_x,
"map_y": domain.remap.map_y,
"valid_mask": domain.valid_mask,
}
if domain.surface is not None:
result["flat_xyz"] = domain.surface.flat_xyz
return result
def flip_sensor_frames_vertical(sensor):
"""Flip captured BGR/gray frames to match Blender-exported map orientation."""
if sensor._img_bgr_raw is not None:
sensor._img_bgr_raw = cv2.flip(sensor._img_bgr_raw, 0)
if sensor.img is not None:
cv2.flip(sensor.img, 0, dst=sensor.img)
def flatten_frame(frame, remap, border_value=0):
flat = cv2.remap(
frame,
remap["map_x"],
remap["map_y"],
interpolation=cv2.INTER_LINEAR,
borderMode=cv2.BORDER_CONSTANT,
borderValue=border_value,
)
flat = flat.copy()
flat[~remap["valid_mask"]] = border_value
return flat
def load_exported_data(data_dir):
mesh = np.load(os.path.join(data_dir, "mesh.npz"))
cam_data = np.load(os.path.join(data_dir, "camera_intrinsics.npz"))
xyz_map_path = os.path.join(data_dir, "pixel_to_xyz.npz")
if not os.path.exists(xyz_map_path):
raise FileNotFoundError(f"Missing precomputed 3D map: {xyz_map_path}")
xyz_map = np.load(xyz_map_path)["xyz_map"]
width = int(cam_data["width"])
height = int(cam_data["height"])
if xyz_map.shape[:2] != (height, width):
raise ValueError(
f"pixel_to_xyz shape {xyz_map.shape[:2]} does not match camera size {(height, width)}"
)
return {
"vertices": mesh["vertices"],
"faces": mesh["faces"],
"uv_data": mesh["uv_data"],
"loop_to_vertex": mesh["loop_to_vertex"],
"K": cam_data["K"],
"width": width,
"height": height,
"world_to_cam": np.load(os.path.join(data_dir, "world_to_cam.npy")),
"xyz_map": xyz_map,
}
def parse_source(value):
value = str(value).strip()
for prefix in ("camera:", "cam:", "webcam:"):
if value.lower().startswith(prefix):
value = value[len(prefix):].strip()
break
try:
return int(value)
except ValueError:
return value
def prepare_static_points(data, display_stride):
xyz_map = data["xyz_map"]
valid_mask = ~np.isnan(xyz_map[:, :, 0])
valid_points = xyz_map[valid_mask].reshape(-1, 3)
if len(valid_points) == 0:
raise ValueError("pixel_to_xyz.npz contains no valid mesh hit points.")
display_indices = np.arange(len(valid_points))[::max(1, int(display_stride))]
return valid_mask, display_indices, valid_points[display_indices]
def align_image_to_map(img_bgr, data, valid_mask, display_indices, flip_vertical=True, flip_horizontal=True):
target_w = int(data["width"])
target_h = int(data["height"])
h, w = img_bgr.shape[:2]
if (w, h) != (target_w, target_h):
img_bgr = cv2.resize(img_bgr, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
img_rgb = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB)
if flip_vertical:
img_rgb = np.flipud(img_rgb).copy()
if flip_horizontal:
img_rgb = np.fliplr(img_rgb).copy()
colors = img_rgb.astype(np.float64) / 255.0
return colors[valid_mask].reshape(-1, 3)[display_indices]
def arrow_map_reference_colors(count, color=(0.04, 0.04, 0.04)):
return np.tile(np.asarray(color, dtype=np.float64), (count, 1))
def magnitude_map_colors(flow, valid_mask, display_indices, max_magnitude=None):
mag = np.linalg.norm(flow, axis=2)
valid_mag = mag[valid_mask].reshape(-1)
display_mag = valid_mag[display_indices]
if max_magnitude is None or max_magnitude <= 0:
max_magnitude = max(float(np.percentile(display_mag, 95)), 1e-6)
value = np.clip(display_mag / max_magnitude, 0.0, 1.0)
return np.column_stack(
(
value,
np.clip(1.0 - np.abs(value - 0.5) * 2.0, 0.0, 1.0),
1.0 - value,
)
).astype(np.float64)
def align_scalar_field_to_map(field, data, flip_vertical=True):
if field.ndim == 3:
field = np.linalg.norm(field, axis=2)
target_w = int(data["width"])
target_h = int(data["height"])
h, w = field.shape[:2]
if (w, h) != (target_w, target_h):
field = cv2.resize(field, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
else:
field = field.copy()
if flip_vertical:
field = np.flipud(field).copy()
return field
def scalar_map_colors(field, valid_mask, display_indices, value_min=None, value_max=None):
valid_values = field[valid_mask].reshape(-1)
display_values = valid_values[display_indices]
finite_values = display_values[np.isfinite(display_values)]
if len(finite_values) == 0:
return arrow_map_reference_colors(len(display_indices))
if value_min is None:
value_min = float(np.percentile(finite_values, 5))
if value_max is None:
value_max = float(np.percentile(finite_values, 95))
if value_max <= value_min:
value_max = value_min + 1e-6
value = np.clip((display_values - value_min) / (value_max - value_min), 0.0, 1.0)
return np.column_stack(
(
value,
np.clip(1.0 - np.abs(value - 0.5) * 2.0, 0.0, 1.0),
1.0 - value,
)
).astype(np.float64)
def align_vector_field_to_map(field, data, flip_vertical=True, scale_as_flow=True):
target_w = int(data["width"])
target_h = int(data["height"])
h, w = field.shape[:2]
if (w, h) != (target_w, target_h):
sx = target_w / max(w, 1)
sy = target_h / max(h, 1)
field = cv2.resize(field, (target_w, target_h), interpolation=cv2.INTER_LINEAR)
if scale_as_flow:
field[:, :, 0] *= sx
field[:, :, 1] *= sy
else:
field = field.copy()
if flip_vertical:
field = np.flipud(field).copy()
field[:, :, 1] *= -1.0
return field
def estimate_surface_normals(xyz_map):
valid = ~np.isnan(xyz_map[:, :, 0])
filled = xyz_map.copy()
filled[~valid] = 0
d_y = np.gradient(filled, axis=0)
d_x = np.gradient(filled, axis=1)
normals = np.cross(d_x, d_y)
norm = np.linalg.norm(normals, axis=2, keepdims=True)
normals = np.divide(normals, norm, out=np.zeros_like(normals), where=norm > 1e-8)
normals[~valid] = 0
return normals
def line_colors(mag, arrow_map_style):
if len(mag) == 0:
return np.empty((0, 3), dtype=np.float64)
if arrow_map_style:
return np.tile(np.array([0.0, 0.9, 1.0], dtype=np.float64), (len(mag), 1))
norm_mag = mag / max(float(np.percentile(mag, 95)), 1e-6)
norm_mag = np.clip(norm_mag, 0.0, 1.0)
return np.column_stack((norm_mag, 1.0 - norm_mag, np.zeros_like(norm_mag))).astype(np.float64)
def lines_from_vectors(starts, vectors, mag, flow_scale, arrow_map_style):
vec_norm = np.linalg.norm(vectors, axis=1)
good = np.isfinite(vec_norm) & (vec_norm > 1e-8)
starts, vectors, mag = starts[good], vectors[good], mag[good]
if len(starts) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
points = np.empty((len(starts) * 2, 3), dtype=np.float64)
points[0::2] = starts
points[1::2] = starts + vectors * flow_scale
lines = np.column_stack((np.arange(0, len(points), 2), np.arange(1, len(points), 2))).astype(np.int32)
return points, lines, line_colors(mag, arrow_map_style)
def build_flow_lines(
xyz_map,
valid_mask,
flow,
arrow_stride=32,
flow_scale=1.0,
min_magnitude=0.5,
arrow_map_style=False,
):
h, w = flow.shape[:2]
ys, xs = np.mgrid[0:h:arrow_stride, 0:w:arrow_stride]
ys = ys.reshape(-1)
xs = xs.reshape(-1)
sample_flow = flow[ys, xs]
mag = np.linalg.norm(sample_flow, axis=1)
keep = valid_mask[ys, xs] & (mag >= min_magnitude)
xs, ys, sample_flow, mag = xs[keep], ys[keep], sample_flow[keep], mag[keep]
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
tx = np.rint(xs + sample_flow[:, 0]).astype(np.int32)
ty = np.rint(ys + sample_flow[:, 1]).astype(np.int32)
inside = (0 <= tx) & (tx < w) & (0 <= ty) & (ty < h)
xs, ys, tx, ty, mag = xs[inside], ys[inside], tx[inside], ty[inside], mag[inside]
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
valid_target = valid_mask[ty, tx]
xs, ys, tx, ty, mag = xs[valid_target], ys[valid_target], tx[valid_target], ty[valid_target], mag[valid_target]
starts = xyz_map[ys, xs]
vectors = xyz_map[ty, tx] - starts
return lines_from_vectors(starts, vectors, mag, flow_scale, arrow_map_style)
def build_shear_normal_lines(
xyz_map,
valid_mask,
vshear,
vnormal,
normals,
arrow_stride=32,
flow_scale=1.0,
min_magnitude=0.01,
normal_gain=1.0,
arrow_map_style=False,
):
h, w = vshear.shape[:2]
ys, xs = np.mgrid[0:h:arrow_stride, 0:w:arrow_stride]
ys = ys.reshape(-1)
xs = xs.reshape(-1)
shear = vshear[ys, xs]
normal_field = vnormal[ys, xs]
shear_mag = np.linalg.norm(shear, axis=1)
normal_mag = np.linalg.norm(normal_field, axis=1) * normal_gain
combined_mag = shear_mag + normal_mag
keep = valid_mask[ys, xs] & (combined_mag >= min_magnitude)
xs, ys, shear, shear_mag, normal_mag, combined_mag = (
xs[keep],
ys[keep],
shear[keep],
shear_mag[keep],
normal_mag[keep],
combined_mag[keep],
)
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
tx = np.rint(xs + shear[:, 0]).astype(np.int32)
ty = np.rint(ys + shear[:, 1]).astype(np.int32)
inside = (0 <= tx) & (tx < w) & (0 <= ty) & (ty < h)
xs, ys, tx, ty, shear_mag, normal_mag, combined_mag = (
xs[inside],
ys[inside],
tx[inside],
ty[inside],
shear_mag[inside],
normal_mag[inside],
combined_mag[inside],
)
if len(xs) == 0:
return np.empty((0, 3)), np.empty((0, 2), dtype=np.int32), np.empty((0, 3))
valid_target = valid_mask[ty, tx]
xs, ys, tx, ty, shear_mag, normal_mag, combined_mag = (
xs[valid_target],
ys[valid_target],
tx[valid_target],
ty[valid_target],
shear_mag[valid_target],
normal_mag[valid_target],
combined_mag[valid_target],
)
starts = xyz_map[ys, xs]
tangent = xyz_map[ty, tx] - starts
tangent_norm = np.linalg.norm(tangent, axis=1, keepdims=True)
tangent_unit = np.divide(tangent, tangent_norm, out=np.zeros_like(tangent), where=tangent_norm > 1e-8)
normal_unit = normals[ys, xs]
shear_ratio = shear_mag / np.maximum(combined_mag, 1e-8)
normal_ratio = normal_mag / np.maximum(combined_mag, 1e-8)
vectors = tangent_unit * shear_ratio[:, None] + normal_unit * normal_ratio[:, None]
vector_norm = np.linalg.norm(vectors, axis=1, keepdims=True)
vectors = np.divide(vectors, vector_norm, out=np.zeros_like(vectors), where=vector_norm > 1e-8)
return lines_from_vectors(starts, vectors, combined_mag, flow_scale, arrow_map_style)
def init_viewer(cloud_points, cloud_colors, flow_points, flow_lines, flow_colors, arrow_map_style=False):
import open3d as o3d
vis = o3d.visualization.Visualizer()
vis.create_window(window_name="SDK 3D Flow", width=1280, height=720)
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(cloud_points.astype(np.float64))
pcd.colors = o3d.utility.Vector3dVector(cloud_colors.astype(np.float64))
vis.add_geometry(pcd)
line_set = o3d.geometry.LineSet()
line_set.points = o3d.utility.Vector3dVector(flow_points.astype(np.float64))
line_set.lines = o3d.utility.Vector2iVector(flow_lines)
line_set.colors = o3d.utility.Vector3dVector(flow_colors.astype(np.float64))
vis.add_geometry(line_set)
render_option = vis.get_render_option()
render_option.point_size = 1.0 if arrow_map_style else 2.0
render_option.line_width = 2.0
render_option.background_color = np.array([0.0, 0.0, 0.0])
vis.poll_events()
vis.update_renderer()
return vis, pcd, line_set, o3d
def run_sdk_3d_flow(args):
data = load_exported_data(args.data_dir)
flatten_remap = load_flatten_remap(args.flatten_remap) if args.flatten_remap else None
view_data = data
if flatten_remap is not None:
if "flat_xyz" not in flatten_remap:
raise ValueError("Flatten remap must contain flat_xyz for 3D visualization.")
view_data = dict(data)
view_data["xyz_map"] = flatten_remap["flat_xyz"]
view_data["height"], view_data["width"] = view_data["xyz_map"].shape[:2]
valid_mask, display_indices, cloud_points = prepare_static_points(view_data, args.display_stride)
surface_normals = estimate_surface_normals(view_data["xyz_map"])
flip_vertical = not args.no_flip_vertical
sensor = orisys.Sensor(
parse_source(args.source),
config_name=args.config,
isstitch=False,
cuda=args.cuda,
verbose=args.verbose,
mirror=not args.no_mirror,
)
vis = pcd = line_set = o3d = None
frame_count = 0
t0 = time.perf_counter()
try:
while True:
if vis is not None and not vis.poll_events():
break
if sensor.get_img() is None:
break
if flip_vertical:
flip_sensor_frames_vertical(sensor)
flat_bgr = None
if flatten_remap is not None:
flat_bgr = flatten_frame(sensor.read_info(sensor.info.IMG)[0], flatten_remap)
flat_gray = cv2.cvtColor(flat_bgr, cv2.COLOR_BGR2GRAY)
sensor.img = flat_gray
sensor._img_bgr_raw = flat_bgr
sensor.img_size_raw = flat_bgr.shape
needs_decomposition = args.vector_mode in ("shear-normal", "depth-map")
sensor.compute_deformation(decompose=needs_decomposition, check_motion=False)
if args.vector_mode == "depth-map":
sensor.compute_contact()
flow, img_bgr, vshear, vnormal, depth_map = sensor.read_info(
sensor.info.VRAW,
sensor.info.IMG,
sensor.info.VSHEAR,
sensor.info.VNORMAL,
sensor.info.DEPTH,
)
if flat_bgr is not None:
img_bgr = flat_bgr
active_data = view_data if flatten_remap is not None else data
# Input is already flipped above; flattened domain uses flat_xyz directly.
flip_for_map = False if (flatten_remap is not None or flip_vertical) else True
flow_map = align_vector_field_to_map(flow, active_data, flip_vertical=flip_for_map)
if args.vector_mode in ("vraw-magnitude", "vraw-magnitude-arrows"):
cloud_colors = magnitude_map_colors(flow_map, valid_mask, display_indices, args.magnitude_max)
elif args.vector_mode == "depth-map":
depth_map = align_scalar_field_to_map(depth_map, active_data, flip_vertical=flip_for_map)
cloud_colors = scalar_map_colors(
depth_map,
valid_mask,
display_indices,
value_min=args.depth_min,
value_max=args.depth_max,
)
elif args.arrow_map_style:
cloud_colors = arrow_map_reference_colors(len(cloud_points))
else:
cloud_colors = align_image_to_map(
img_bgr,
active_data,
valid_mask,
display_indices,
flip_vertical=flip_for_map,
)
if args.vector_mode == "shear-normal":
vshear_map = align_vector_field_to_map(vshear, active_data, flip_vertical=flip_for_map)
vnormal_map = align_vector_field_to_map(vnormal, active_data, flip_vertical=flip_for_map)
flow_points, flow_lines, flow_colors = build_shear_normal_lines(
active_data["xyz_map"],
valid_mask,
vshear_map,
vnormal_map,
surface_normals,
arrow_stride=args.arrow_stride,
flow_scale=args.flow_scale,
min_magnitude=args.min_flow,
normal_gain=args.normal_gain,
arrow_map_style=args.arrow_map_style,
)
elif args.vector_mode in ("vraw-magnitude", "depth-map"):
flow_points = np.empty((0, 3), dtype=np.float64)
flow_lines = np.empty((0, 2), dtype=np.int32)
flow_colors = np.empty((0, 3), dtype=np.float64)
else:
flow_points, flow_lines, flow_colors = build_flow_lines(
active_data["xyz_map"],
valid_mask,
flow_map,
arrow_stride=args.arrow_stride,
flow_scale=args.flow_scale,
min_magnitude=args.min_flow,
arrow_map_style=args.arrow_map_style,
)
if vis is None:
vis, pcd, line_set, o3d = init_viewer(
cloud_points,
cloud_colors,
flow_points,
flow_lines,
flow_colors,
arrow_map_style=args.arrow_map_style,
)
print("SDK 3D flow viewer started. Close the Open3D window to stop.")
else:
pcd.colors = o3d.utility.Vector3dVector(cloud_colors)
line_set.points = o3d.utility.Vector3dVector(flow_points.astype(np.float64))
line_set.lines = o3d.utility.Vector2iVector(flow_lines)
line_set.colors = o3d.utility.Vector3dVector(flow_colors.astype(np.float64))
vis.update_geometry(pcd)
vis.update_geometry(line_set)
vis.update_renderer()
frame_count += 1
if frame_count % 30 == 0:
elapsed = max(time.perf_counter() - t0, 1e-6)
print(f"frames={frame_count}, viewer FPS={frame_count / elapsed:.1f}, arrows={len(flow_lines)}")
finally:
sensor.disconnect()
if vis is not None:
vis.destroy_window()
def parse_args():
parser = argparse.ArgumentParser(description="Map SDK flow fields to a 3D point cloud.")
parser.add_argument("--data-dir", default=os.path.join(SCRIPT_DIR, "exported_data"))
parser.add_argument(
"--source",
default=os.path.join(SCRIPT_DIR, "20260602-094916.mp4"),
help="Video path or live camera ID, for example 0, 1, camera:0, or cam:1.",
)
parser.add_argument("--config", default=os.path.join(REPO_ROOT, "config", "dd02-ov.json"))
parser.add_argument("--display-stride", type=int, default=4)
parser.add_argument("--arrow-stride", type=int, default=32)
parser.add_argument("--flow-scale", type=float, default=3.0)
parser.add_argument("--min-flow", type=float, default=0.5)
parser.add_argument("--normal-gain", type=float, default=1.0)
parser.add_argument("--magnitude-max", type=float, default=None,
help="Upper VRAW magnitude for heat-map color scaling; default uses frame percentile.")
parser.add_argument("--depth-min", type=float, default=None,
help="Minimum depth value for heat-map color scaling; default uses frame percentile.")
parser.add_argument("--depth-max", type=float, default=None,
help="Maximum depth value for heat-map color scaling; default uses frame percentile.")
parser.add_argument("--flatten-remap", default=None,
help="Optional flatten_remap.npz to run optical flow in UV-flattened image space.")
parser.add_argument(
"--vector-mode",
choices=["vraw", "vraw-magnitude", "vraw-magnitude-arrows", "shear-normal", "depth-map"],
default="vraw",
)
parser.add_argument(
"--arrow-map-style",
action="store_true",
help="Hide real texture and use a dark Qt-viewer-like arrow map style.",
)
parser.add_argument(
"--no-flip-vertical",
action="store_true",
help="Do not flip captured video vertically before processing.",
)
parser.add_argument("--no-mirror", action="store_true")
parser.add_argument("--cuda", action="store_true")
parser.add_argument("--verbose", action="store_true")
return parser.parse_args()
if __name__ == "__main__":
run_sdk_3d_flow(parse_args())
Binary file not shown.

Before

Width:  |  Height:  |  Size: 1.2 MiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 795 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 898 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 1002 KiB

-4
View File
@@ -1,4 +0,0 @@
"""Backward-compatible re-exports; prefer ``examples/viz/mesh_viewer.py`` or ``orisys.finger``."""
from orisys.finger.mesh import * # noqa: F403
from orisys.finger.viewer import * # noqa: F403
File diff suppressed because one or more lines are too long

Before

Width:  |  Height:  |  Size: 74 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 173 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 929 KiB

Binary file not shown.
Binary file not shown.

Before

Width:  |  Height:  |  Size: 626 KiB

File diff suppressed because one or more lines are too long

Before

Width:  |  Height:  |  Size: 2.6 MiB

File diff suppressed because one or more lines are too long

Before

Width:  |  Height:  |  Size: 247 KiB

@@ -1,53 +0,0 @@
#%%
import numpy as np
import cv2
def generate_color_pattern(h, w, d, r):
# Initialize the pattern with zeros
pattern = np.zeros((h, w, 3))
# Calculate the number of patches in each dimension
num_patches_x = h // d
num_patches_y = w // d
# Maximum number of attempts to find a suitable color
max_attempts = 1000
# Iterate over the patches
for i in range(num_patches_x):
for j in range(num_patches_y):
# Get the neighbors
neighbors = []
if i > 0:
neighbors.append(pattern[(i-1)*d:i*d, j*d:(j+1)*d])
if j > 0:
neighbors.append(pattern[i*d:(i+1)*d, (j-1)*d:j*d])
# Generate a new color
for attempt in range(max_attempts):
new_color = np.random.rand(3)
if all(np.min((new_color - neighbor.mean(axis=(0, 1)))**2) >= r for neighbor in neighbors):
break
if attempt == max_attempts - 1:
print(f"Warning: Could not find a suitable color for patch ({i}, {j}) after {max_attempts} attempts")
# Apply the new color to the patch
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
print('Patch ({}, {})'.format(i, j))
return pattern
# Test the function
pattern = generate_color_pattern(162, 162, 3, 0.75)
# Convert the pattern to a format suitable for display with OpenCV
pattern_display = (pattern * 255).astype(np.uint8)
# # Display the pattern
# cv2.imshow('Pattern', pattern_display)
# cv2.waitKey(0)
# cv2.destroyAllWindows()
# #%%
cv2.imwrite('color_pattern.png', pattern_display)
@@ -1,773 +0,0 @@
#%%
import numpy as np
import cv2
import matplotlib.pyplot as plt
from PIL import Image
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def generate_color_pattern(h, w, d, r):
"""
Generate a dense color pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
"""
pattern = np.zeros((h, w, 3), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch colors (using mean because patches have uniform color).
neighbors = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
neighbors.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
neighbors.append(np.mean(patch_left, axis=(0, 1)))
if len(neighbors) == 0:
# For the first patch use a random color.
new_color = np.random.rand(3)
else:
new_color = generate_valid_color(neighbors, r)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
# print('Patch ({}, {})'.format(i, j))
return pattern
def convert_to_grayscale(color_pattern):
"""
Convert color pattern to grayscale using standard luminance weights.
Parameters:
color_pattern: Input color image (H x W x 3) with values in range [0, 1]
Returns:
Grayscale image (H x W) with values in range [0, 1]
"""
# Use standard RGB to grayscale conversion: 0.299*R + 0.587*G + 0.114*B
grayscale = 0.299 * color_pattern[:, :, 2] + 0.587 * color_pattern[:, :, 1] + 0.114 * color_pattern[:, :, 0]
return grayscale
def optimize_grayscale_patches(grayscale, d, min_diff=0.1, max_iterations=100):
"""
Optimize grayscale image to improve value difference between neighboring patches.
Parameters:
grayscale: Input grayscale image (H x W) with values in range [0, 1]
d: patch size in pixels (each patch is d x d)
min_diff: minimum absolute difference required between neighboring patches
max_iterations: maximum number of optimization iterations
Returns:
Optimized grayscale image with improved patch differences
"""
h, w = grayscale.shape
num_patches_x = h // d
num_patches_y = w // d
# Create optimized version
optimized = grayscale.copy()
for iteration in range(max_iterations):
improved = False
for i in range(num_patches_x):
for j in range(num_patches_y):
# Get current patch value (mean of patch)
current_patch = optimized[i*d:(i+1)*d, j*d:(j+1)*d]
current_value = np.mean(current_patch)
# Collect neighbor patch values
neighbor_values = []
if i > 0:
patch_above = optimized[(i-1)*d:i*d, j*d:(j+1)*d]
neighbor_values.append(np.mean(patch_above))
if j > 0:
patch_left = optimized[i*d:(i+1)*d, (j-1)*d:j*d]
neighbor_values.append(np.mean(patch_left))
if i < num_patches_x - 1:
patch_below = optimized[(i+1)*d:(i+2)*d, j*d:(j+1)*d]
neighbor_values.append(np.mean(patch_below))
if j < num_patches_y - 1:
patch_right = optimized[i*d:(i+1)*d, (j+1)*d:(j+2)*d]
neighbor_values.append(np.mean(patch_right))
if len(neighbor_values) == 0:
continue
# Check if current value is too close to any neighbor
min_neighbor_diff = min([abs(current_value - nv) for nv in neighbor_values])
if min_neighbor_diff < min_diff:
# Need to adjust this patch
# Find valid range that maintains min_diff from all neighbors
forbidden = []
for nv in neighbor_values:
forbidden.append((max(0, nv - min_diff), min(1, nv + min_diff)))
# Merge forbidden intervals
merged_forbidden = merge_intervals(forbidden)
# Find allowed intervals
allowed = []
if merged_forbidden:
if merged_forbidden[0][0] > 0:
allowed.append((0, merged_forbidden[0][0]))
for k in range(1, len(merged_forbidden)):
gap = merged_forbidden[k][0] - merged_forbidden[k-1][1]
if gap > 0:
allowed.append((merged_forbidden[k-1][1], merged_forbidden[k][0]))
if merged_forbidden[-1][1] < 1:
allowed.append((merged_forbidden[-1][1], 1))
else:
allowed.append((0, 1))
# Sample new value from allowed intervals
if allowed:
new_value = sample_from_allowed(allowed)
if new_value is not None:
optimized[i*d:(i+1)*d, j*d:(j+1)*d] = new_value
improved = True
if not improved:
print(f"Converged after {iteration + 1} iterations")
break
return optimized
def analyze_patch_differences(image, d):
"""
Analyze and report the minimum differences between neighboring patches.
Parameters:
image: Grayscale image (H x W) with values in range [0, 1]
d: patch size in pixels
Returns:
Dictionary with statistics about patch differences
"""
h, w = image.shape
num_patches_x = h // d
num_patches_y = w // d
horizontal_diffs = []
vertical_diffs = []
for i in range(num_patches_x):
for j in range(num_patches_y):
current_patch = image[i*d:(i+1)*d, j*d:(j+1)*d]
current_value = np.mean(current_patch)
# Check horizontal neighbor (right)
if j < num_patches_y - 1:
patch_right = image[i*d:(i+1)*d, (j+1)*d:(j+2)*d]
right_value = np.mean(patch_right)
horizontal_diffs.append(abs(current_value - right_value))
# Check vertical neighbor (below)
if i < num_patches_x - 1:
patch_below = image[(i+1)*d:(i+2)*d, j*d:(j+1)*d]
below_value = np.mean(patch_below)
vertical_diffs.append(abs(current_value - below_value))
stats = {
'min_diff': min(horizontal_diffs + vertical_diffs) if (horizontal_diffs + vertical_diffs) else 0,
'max_diff': max(horizontal_diffs + vertical_diffs) if (horizontal_diffs + vertical_diffs) else 0,
'mean_diff': np.mean(horizontal_diffs + vertical_diffs) if (horizontal_diffs + vertical_diffs) else 0,
'median_diff': np.median(horizontal_diffs + vertical_diffs) if (horizontal_diffs + vertical_diffs) else 0,
}
return stats
def rgb_to_cmyk(rgb_image):
"""
Convert RGB image to CMYK color space for printing.
Parameters:
rgb_image: RGB image (H x W x 3) with values in range [0, 1], BGR order (OpenCV format)
Returns:
CMYK image (H x W x 4) with values in range [0, 1]
"""
# OpenCV uses BGR, so we need to handle it as B, G, R in channels 0, 1, 2
b = rgb_image[:, :, 0]
g = rgb_image[:, :, 1]
r = rgb_image[:, :, 2]
# Calculate K (black) channel
k = 1 - np.maximum(np.maximum(r, g), b)
# Avoid division by zero
k_inv = 1 - k
k_inv = np.where(k_inv == 0, 1e-10, k_inv)
# Calculate CMY channels
c = (1 - r - k) / k_inv
m = (1 - g - k) / k_inv
y = (1 - b - k) / k_inv
# Stack into CMYK image
cmyk = np.stack([c, m, y, k], axis=2)
cmyk = np.clip(cmyk, 0, 1)
return cmyk
def cmyk_to_rgb(cmyk_image):
"""
Convert CMYK image back to RGB color space.
Parameters:
cmyk_image: CMYK image (H x W x 4) with values in range [0, 1]
Returns:
RGB image (H x W x 3) with values in range [0, 1], BGR order (OpenCV format)
"""
c = cmyk_image[:, :, 0]
m = cmyk_image[:, :, 1]
y = cmyk_image[:, :, 2]
k = cmyk_image[:, :, 3]
# Convert back to RGB
r = (1 - c) * (1 - k)
g = (1 - m) * (1 - k)
b = (1 - y) * (1 - k)
# Stack in BGR order for OpenCV
rgb = np.stack([b, g, r], axis=2)
rgb = np.clip(rgb, 0, 1)
return rgb
def optimize_cmyk_patches(cmyk_image, d, min_diff=0.1, max_iterations=100):
"""
Optimize CMYK image to improve color difference between neighboring patches.
Uses Euclidean distance in CMYK space.
Parameters:
cmyk_image: Input CMYK image (H x W x 4) with values in range [0, 1]
d: patch size in pixels (each patch is d x d)
min_diff: minimum Euclidean distance required between neighboring patches in CMYK space
max_iterations: maximum number of optimization iterations
Returns:
Optimized CMYK image with improved patch differences
"""
h, w = cmyk_image.shape[:2]
num_patches_x = h // d
num_patches_y = w // d
optimized = cmyk_image.copy()
for iteration in range(max_iterations):
improved = False
for i in range(num_patches_x):
for j in range(num_patches_y):
# Get current patch CMYK values (mean of patch)
current_patch = optimized[i*d:(i+1)*d, j*d:(j+1)*d]
current_cmyk = np.mean(current_patch, axis=(0, 1))
# Collect neighbor patch CMYK values
neighbor_cmyk = []
if i > 0:
patch_above = optimized[(i-1)*d:i*d, j*d:(j+1)*d]
neighbor_cmyk.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = optimized[i*d:(i+1)*d, (j-1)*d:j*d]
neighbor_cmyk.append(np.mean(patch_left, axis=(0, 1)))
if i < num_patches_x - 1:
patch_below = optimized[(i+1)*d:(i+2)*d, j*d:(j+1)*d]
neighbor_cmyk.append(np.mean(patch_below, axis=(0, 1)))
if j < num_patches_y - 1:
patch_right = optimized[i*d:(i+1)*d, (j+1)*d:(j+2)*d]
neighbor_cmyk.append(np.mean(patch_right, axis=(0, 1)))
if len(neighbor_cmyk) == 0:
continue
# Calculate minimum Euclidean distance to neighbors in CMYK space
distances = [np.linalg.norm(current_cmyk - n) for n in neighbor_cmyk]
min_distance = min(distances)
if min_distance < min_diff:
# Need to adjust this patch - try to move it away from neighbors
# Find the direction away from the closest neighbor
closest_idx = np.argmin(distances)
closest_neighbor = neighbor_cmyk[closest_idx]
# Move away from closest neighbor
direction = current_cmyk - closest_neighbor
direction_norm = np.linalg.norm(direction)
if direction_norm > 0:
direction = direction / direction_norm
# Move by the deficit amount
move_amount = (min_diff - min_distance) * 1.2 # 1.2 factor for faster convergence
new_cmyk = current_cmyk + direction * move_amount
new_cmyk = np.clip(new_cmyk, 0, 1)
# Check if new value is better
new_distances = [np.linalg.norm(new_cmyk - n) for n in neighbor_cmyk]
new_min_distance = min(new_distances)
if new_min_distance > min_distance:
optimized[i*d:(i+1)*d, j*d:(j+1)*d] = new_cmyk
improved = True
if not improved:
print(f"CMYK optimization converged after {iteration + 1} iterations")
break
return optimized
def analyze_cmyk_patch_differences(cmyk_image, d):
"""
Analyze CMYK patch differences using Euclidean distance in CMYK space.
Parameters:
cmyk_image: CMYK image (H x W x 4) with values in range [0, 1]
d: patch size in pixels
Returns:
Dictionary with statistics about CMYK patch differences
"""
h, w = cmyk_image.shape[:2]
num_patches_x = h // d
num_patches_y = w // d
horizontal_diffs = []
vertical_diffs = []
for i in range(num_patches_x):
for j in range(num_patches_y):
current_patch = cmyk_image[i*d:(i+1)*d, j*d:(j+1)*d]
current_cmyk = np.mean(current_patch, axis=(0, 1))
# Check horizontal neighbor (right)
if j < num_patches_y - 1:
patch_right = cmyk_image[i*d:(i+1)*d, (j+1)*d:(j+2)*d]
right_cmyk = np.mean(patch_right, axis=(0, 1))
horizontal_diffs.append(np.linalg.norm(current_cmyk - right_cmyk))
# Check vertical neighbor (below)
if i < num_patches_x - 1:
patch_below = cmyk_image[(i+1)*d:(i+2)*d, j*d:(j+1)*d]
below_cmyk = np.mean(patch_below, axis=(0, 1))
vertical_diffs.append(np.linalg.norm(current_cmyk - below_cmyk))
all_diffs = horizontal_diffs + vertical_diffs
stats = {
'min_diff': min(all_diffs) if all_diffs else 0,
'max_diff': max(all_diffs) if all_diffs else 0,
'mean_diff': np.mean(all_diffs) if all_diffs else 0,
'median_diff': np.median(all_diffs) if all_diffs else 0,
}
return stats
def grayscale_to_cmyk_optimized(grayscale_target, d, min_cmyk_diff=0.2, max_iterations=100, color_diversity=0.3):
"""
Generate a CMYK image that when converted to grayscale matches the target grayscale,
while maintaining good color separation between patches.
Parameters:
grayscale_target: Target grayscale image (H x W) with values in range [0, 1]
d: patch size in pixels
min_cmyk_diff: minimum Euclidean distance required between neighboring patches in CMYK space
max_iterations: maximum number of optimization iterations
color_diversity: amount of color diversity to inject (0.0 to 1.0)
Returns:
CMYK image (H x W x 4) optimized for printing
"""
h, w = grayscale_target.shape
num_patches_x = h // d
num_patches_y = w // d
# Initialize CMYK image
cmyk_image = np.zeros((h, w, 4), dtype=np.float32)
# Generate initial CMYK values for each patch
for i in range(num_patches_x):
for j in range(num_patches_y):
# Get target grayscale value for this patch
target_gray = np.mean(grayscale_target[i*d:(i+1)*d, j*d:(j+1)*d])
# Generate random CMY values with diversity
c = np.random.rand() * color_diversity
m = np.random.rand() * color_diversity
y = np.random.rand() * color_diversity
# Calculate K based on target grayscale
# Grayscale from CMYK: gray ≈ 1 - K when CMY are small
# For more accurate: gray = (1-C)*(1-K)*0.114 + (1-M)*(1-K)*0.587 + (1-Y)*(1-K)*0.299
k = 1 - target_gray # Initial approximation
# Adjust to match target grayscale
# Create temporary CMYK and test
temp_cmyk = np.array([c, m, y, k])
temp_rgb = np.array([(1-c)*(1-k), (1-m)*(1-k), (1-y)*(1-k)]) # R, G, B
temp_gray = 0.299 * temp_rgb[0] + 0.587 * temp_rgb[1] + 0.114 * temp_rgb[2]
# Adjust K to match target
k_adjustment = target_gray - temp_gray
k = np.clip(k - k_adjustment, 0, 1)
cmyk_image[i*d:(i+1)*d, j*d:(j+1)*d] = [c, m, y, k]
# Optimize CMYK to improve patch differences while maintaining grayscale target
print("Optimizing CMYK for color separation while matching grayscale...")
for iteration in range(max_iterations):
improved = False
for i in range(num_patches_x):
for j in range(num_patches_y):
# Get current patch
current_patch = cmyk_image[i*d:(i+1)*d, j*d:(j+1)*d]
current_cmyk = np.mean(current_patch, axis=(0, 1))
# Get target grayscale for this patch
target_gray = np.mean(grayscale_target[i*d:(i+1)*d, j*d:(j+1)*d])
# Collect neighbor CMYK values
neighbor_cmyk = []
if i > 0:
neighbor_cmyk.append(np.mean(cmyk_image[(i-1)*d:i*d, j*d:(j+1)*d], axis=(0, 1)))
if j > 0:
neighbor_cmyk.append(np.mean(cmyk_image[i*d:(i+1)*d, (j-1)*d:j*d], axis=(0, 1)))
if i < num_patches_x - 1:
neighbor_cmyk.append(np.mean(cmyk_image[(i+1)*d:(i+2)*d, j*d:(j+1)*d], axis=(0, 1)))
if j < num_patches_y - 1:
neighbor_cmyk.append(np.mean(cmyk_image[i*d:(i+1)*d, (j+1)*d:(j+2)*d], axis=(0, 1)))
if len(neighbor_cmyk) == 0:
continue
# Check CMYK distance to neighbors
distances = [np.linalg.norm(current_cmyk - n) for n in neighbor_cmyk]
min_distance = min(distances)
if min_distance < min_cmyk_diff:
# Try to improve color separation
# Vary CMY channels while keeping grayscale constant
for attempt in range(5):
# Generate new CMY values
new_c = np.clip(current_cmyk[0] + (np.random.rand() - 0.5) * 0.2, 0, 1)
new_m = np.clip(current_cmyk[1] + (np.random.rand() - 0.5) * 0.2, 0, 1)
new_y = np.clip(current_cmyk[2] + (np.random.rand() - 0.5) * 0.2, 0, 1)
# Calculate K to maintain grayscale
# gray = (1-C)*(1-K)*0.114 + (1-M)*(1-K)*0.587 + (1-Y)*(1-K)*0.299
# gray = (1-K) * [0.114*(1-C) + 0.587*(1-M) + 0.299*(1-Y)]
rgb_factor = 0.114*(1-new_c) + 0.587*(1-new_m) + 0.299*(1-new_y)
if rgb_factor > 0:
new_k = 1 - (target_gray / rgb_factor)
new_k = np.clip(new_k, 0, 1)
else:
new_k = current_cmyk[3]
new_cmyk = np.array([new_c, new_m, new_y, new_k])
# Verify grayscale match
test_rgb = np.array([(1-new_c)*(1-new_k), (1-new_m)*(1-new_k), (1-new_y)*(1-new_k)])
test_gray = 0.299 * test_rgb[0] + 0.587 * test_rgb[1] + 0.114 * test_rgb[2]
# Check if grayscale is close enough and CMYK distance improved
if abs(test_gray - target_gray) < 0.05: # Grayscale tolerance
new_distances = [np.linalg.norm(new_cmyk - n) for n in neighbor_cmyk]
new_min_distance = min(new_distances)
if new_min_distance > min_distance:
cmyk_image[i*d:(i+1)*d, j*d:(j+1)*d] = new_cmyk
improved = True
break
if not improved:
print(f"Grayscale-to-CMYK optimization converged after {iteration + 1} iterations")
break
return cmyk_image
def cmyk_to_grayscale(cmyk_image):
"""
Convert CMYK image to grayscale.
Parameters:
cmyk_image: CMYK image (H x W x 4) with values in range [0, 1]
Returns:
Grayscale image (H x W) with values in range [0, 1]
"""
# Convert to RGB first, then to grayscale
rgb = cmyk_to_rgb(cmyk_image)
grayscale = 0.299 * rgb[:, :, 2] + 0.587 * rgb[:, :, 1] + 0.114 * rgb[:, :, 0]
return grayscale
def main():
# Parameters for the pattern.
# Image dimensions in millimeters
width_mm = 50 # Width in millimeters
height_mm = 50 # Height in millimeters
dpi = 600 # Dots per inch for printing
patch_size = 1 # 3 for 300dpi, 6 for 600dpi
# Convert mm to pixels: pixels = (mm / 25.4) * dpi
width = int((width_mm / 25.4) * dpi)
height = int((height_mm / 25.4) * dpi)
print(f"Image dimensions: {width_mm}mm x {height_mm}mm at {dpi} DPI")
print(f"Pixel dimensions: {width} x {height} pixels")
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1
# Step 1: Generate initial color pattern
print("Generating color pattern...")
pattern = generate_color_pattern(height, width, patch_size, r)
pattern_display = (pattern * 255).astype(np.uint8)
# Step 2: Convert to grayscale
print("\nConverting to grayscale...")
grayscale = convert_to_grayscale(pattern)
grayscale_display = (grayscale * 255).astype(np.uint8)
# Step 3: Analyze original grayscale differences
print("\nAnalyzing original grayscale patch differences...")
original_stats = analyze_patch_differences(grayscale, patch_size)
print(f"Original - Min: {original_stats['min_diff']:.4f}, Max: {original_stats['max_diff']:.4f}, "
f"Mean: {original_stats['mean_diff']:.4f}, Median: {original_stats['median_diff']:.4f}")
# Step 4: Optimize grayscale to improve patch differences
print("\nOptimizing grayscale patch differences...")
min_diff_threshold = 0.15 # Adjust this value to control minimum difference
optimized_grayscale = optimize_grayscale_patches(grayscale, patch_size, min_diff=min_diff_threshold, max_iterations=50)
optimized_display = (optimized_grayscale * 255).astype(np.uint8)
# Step 5: Analyze optimized grayscale differences
print("\nAnalyzing optimized grayscale patch differences...")
optimized_stats = analyze_patch_differences(optimized_grayscale, patch_size)
print(f"Optimized - Min: {optimized_stats['min_diff']:.4f}, Max: {optimized_stats['max_diff']:.4f}, "
f"Mean: {optimized_stats['mean_diff']:.4f}, Median: {optimized_stats['median_diff']:.4f}")
# Step 6: Generate CMYK image from optimized grayscale
print("\n" + "="*60)
print("GENERATING CMYK FROM OPTIMIZED GRAYSCALE")
print("="*60)
print("\nGenerating CMYK image to match optimized grayscale...")
cmyk_from_grayscale = grayscale_to_cmyk_optimized(
optimized_grayscale,
patch_size,
min_cmyk_diff=0.2,
max_iterations=50,
color_diversity=0.5
)
# Analyze CMYK differences
print("\nAnalyzing CMYK patch differences...")
cmyk_stats = analyze_cmyk_patch_differences(cmyk_from_grayscale, patch_size)
print(f"CMYK - Min: {cmyk_stats['min_diff']:.4f}, Max: {cmyk_stats['max_diff']:.4f}, "
f"Mean: {cmyk_stats['mean_diff']:.4f}, Median: {cmyk_stats['median_diff']:.4f}")
# Convert CMYK back to RGB for display
rgb_from_cmyk = cmyk_to_rgb(cmyk_from_grayscale)
rgb_from_cmyk_display = (rgb_from_cmyk * 255).astype(np.uint8)
# Convert CMYK back to grayscale to verify
grayscale_from_cmyk = cmyk_to_grayscale(cmyk_from_grayscale)
grayscale_from_cmyk_display = (grayscale_from_cmyk * 255).astype(np.uint8)
# Analyze how well grayscale was preserved
grayscale_from_cmyk_stats = analyze_patch_differences(grayscale_from_cmyk, patch_size)
print(f"Grayscale from CMYK - Min: {grayscale_from_cmyk_stats['min_diff']:.4f}")
# Calculate grayscale reconstruction error
gray_error = np.mean(np.abs(optimized_grayscale - grayscale_from_cmyk))
print(f"Grayscale reconstruction error: {gray_error:.4f}")
# Display all images with matplotlib
fig, axes = plt.subplots(2, 3, figsize=(18, 12), dpi=300)
# Row 1: Original grayscale optimization
# Color pattern (convert BGR to RGB for matplotlib)
pattern_rgb = cv2.cvtColor(pattern_display, cv2.COLOR_BGR2RGB)
axes[0, 0].imshow(pattern_rgb)
axes[0, 0].set_title(f'1. Original Color Pattern\nSize: {height}x{width} pixels', fontsize=10)
axes[0, 0].axis('on')
axes[0, 0].grid(True, alpha=0.3)
# Grayscale from color
axes[0, 1].imshow(grayscale_display, cmap='gray', vmin=0, vmax=255)
axes[0, 1].set_title(f'2. Grayscale (from color)\nSize: {height}x{width} pixels\nMin diff: {original_stats["min_diff"]:.4f}', fontsize=10)
axes[0, 1].axis('on')
axes[0, 1].grid(True, alpha=0.3)
# Optimized grayscale
axes[0, 2].imshow(optimized_display, cmap='gray', vmin=0, vmax=255)
axes[0, 2].set_title(f'3. Optimized Grayscale (Target)\nSize: {height}x{width} pixels\nMin diff: {optimized_stats["min_diff"]:.4f}', fontsize=10, fontweight='bold')
axes[0, 2].axis('on')
axes[0, 2].grid(True, alpha=0.3)
# Row 2: CMYK generation from optimized grayscale
# RGB from CMYK
rgb_from_cmyk_viz = cv2.cvtColor(rgb_from_cmyk_display, cv2.COLOR_BGR2RGB)
axes[1, 0].imshow(rgb_from_cmyk_viz)
axes[1, 0].set_title(f'4. RGB from Generated CMYK\nCMYK Min diff: {cmyk_stats["min_diff"]:.4f}\n(Best for printing)', fontsize=10, fontweight='bold')
axes[1, 0].axis('on')
axes[1, 0].grid(True, alpha=0.3)
# Grayscale from CMYK (verification)
axes[1, 1].imshow(grayscale_from_cmyk_display, cmap='gray', vmin=0, vmax=255)
axes[1, 1].set_title(f'5. Grayscale from CMYK\nMin diff: {grayscale_from_cmyk_stats["min_diff"]:.4f}\nError: {gray_error:.4f}', fontsize=10)
axes[1, 1].axis('on')
axes[1, 1].grid(True, alpha=0.3)
# Show CMYK channels visualization
cmyk_viz = np.zeros((height, width, 3), dtype=np.uint8)
cmyk_viz[:, :, 0] = (cmyk_from_grayscale[:, :, 0] * 255).astype(np.uint8) # Cyan in B channel
cmyk_viz[:, :, 1] = (cmyk_from_grayscale[:, :, 1] * 255).astype(np.uint8) # Magenta in G channel
cmyk_viz[:, :, 2] = (cmyk_from_grayscale[:, :, 2] * 255).astype(np.uint8) # Yellow in R channel
axes[1, 2].imshow(cv2.cvtColor(cmyk_viz, cv2.COLOR_BGR2RGB))
axes[1, 2].set_title(f'6. CMYK Channels (CMY only)\nK channel separate', fontsize=10)
axes[1, 2].axis('on')
axes[1, 2].grid(True, alpha=0.3)
plt.tight_layout()
print(f"\nDisplay info: Figure DPI={fig.dpi}, Size={fig.get_size_inches()} inches")
plt.show()
# Save all images
cv2.imwrite('color_pattern.png', pattern_display)
cv2.imwrite('grayscale_from_color.png', grayscale_display)
cv2.imwrite('optimized_grayscale.png', optimized_display)
cv2.imwrite('cmyk_generated_rgb.png', rgb_from_cmyk_display)
cv2.imwrite('grayscale_from_cmyk.png', grayscale_from_cmyk_display)
# Save CMYK data
np.save('cmyk_from_grayscale.npy', cmyk_from_grayscale)
# Save CMYK as TIFF (proper format for printing)
print("\nSaving CMYK as TIFF...")
try:
# Convert CMYK to uint8 (0-255 range)
cmyk_uint8 = (cmyk_from_grayscale * 255).astype(np.uint8)
# For TIFF CMYK, we need to ensure proper channel order
# PIL CMYK expects channels in C, M, Y, K order
# Split channels
c_channel = cmyk_uint8[:, :, 0]
m_channel = cmyk_uint8[:, :, 1]
y_channel = cmyk_uint8[:, :, 2]
k_channel = cmyk_uint8[:, :, 3]
# Create PIL images for each channel and merge
c_img = Image.fromarray(c_channel, mode='L')
m_img = Image.fromarray(m_channel, mode='L')
y_img = Image.fromarray(y_channel, mode='L')
k_img = Image.fromarray(k_channel, mode='L')
# Merge into CMYK image
cmyk_pil = Image.merge('CMYK', (c_img, m_img, y_img, k_img))
# Save as TIFF with CMYK color space at 300 DPI for professional printing
cmyk_pil.save('cmyk_pattern.tiff', format='TIFF', dpi=(300, 300))
# Verify the TIFF was saved by checking file size
import os
if os.path.exists('cmyk_pattern.tiff'):
file_size = os.path.getsize('cmyk_pattern.tiff')
print(f"CMYK TIFF saved successfully: {file_size} bytes at 300 DPI")
except Exception as e:
print(f"Error saving CMYK TIFF: {e}")
print("\nImages saved successfully!")
print("CMYK-generated RGB: cmyk_generated_rgb.png")
print("Grayscale verification: grayscale_from_cmyk.png")
print("CMYK data: cmyk_from_grayscale.npy")
print("CMYK TIFF (300 DPI): cmyk_pattern.tiff (professional printing format)")
if __name__ == '__main__':
main()
@@ -1,140 +0,0 @@
#%%
import numpy as np
import cv2
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def generate_color_pattern(h, w, d, r):
"""
Generate a dense color pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
"""
pattern = np.zeros((h, w, 3), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch colors (using mean because patches have uniform color).
neighbors = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
neighbors.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
neighbors.append(np.mean(patch_left, axis=(0, 1)))
if len(neighbors) == 0:
# For the first patch use a random color.
new_color = np.random.rand(3)
else:
new_color = generate_valid_color(neighbors, r)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
print('Patch ({}, {})'.format(i, j))
return pattern
def main():
# Parameters for the pattern.
# For example, using 100x100 pixels image with patches of 3x3 pixels.
# 11111111
height = 1200
width = 1000
patch_size = 12
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1000
pattern = generate_color_pattern(height, width, patch_size, r)
# Scale to 0-255 for display and saving.
pattern_display = (pattern * 255).astype(np.uint8)
# Display the pattern with OpenCV.
cv2.imshow('Pattern', pattern_display)
cv2.waitKey(0)
cv2.destroyAllWindows()
# Save the pattern as an image file.
cv2.imwrite('color_pattern.png', pattern_display)
if __name__ == '__main__':
main()
@@ -1,46 +0,0 @@
#%%
import numpy as np
import cv2
def generate_checker_pattern(height, width, patch_size):
"""
Generate a black-and-white checkerboard pattern.
Each `patch_size x patch_size` block alternates between black and white.
Partial blocks at the image borders are also filled correctly.
"""
pattern = np.zeros((height, width, 3), dtype=np.uint8)
num_patches_x = (height + patch_size - 1) // patch_size
num_patches_y = (width + patch_size - 1) // patch_size
for i in range(num_patches_x):
for j in range(num_patches_y):
color = 255 if (i + j) % 2 == 0 else 0
y0 = i * patch_size
y1 = min((i + 1) * patch_size, height)
x0 = j * patch_size
x1 = min((j + 1) * patch_size, width)
pattern[y0:y1, x0:x1] = color
print("Patch ({}, {}) -> {}".format(i, j, "white" if color == 255 else "black"))
return pattern
def main():
height = 1200
width = 1000
patch_size = 12
pattern = generate_checker_pattern(height, width, patch_size)
cv2.imshow("Checker Pattern", pattern)
cv2.waitKey(0)
cv2.destroyAllWindows()
cv2.imwrite("checker_pattern.png", pattern)
if __name__ == '__main__':
main()
-139
View File
@@ -1,139 +0,0 @@
#%%
import numpy as np
import cv2
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def generate_color_pattern(h, w, d, r):
"""
Generate a dense color pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
"""
pattern = np.zeros((h, w, 3), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch colors (using mean because patches have uniform color).
neighbors = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
neighbors.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
neighbors.append(np.mean(patch_left, axis=(0, 1)))
if len(neighbors) == 0:
# For the first patch use a random color.
new_color = np.random.rand(3)
else:
new_color = generate_valid_color(neighbors, r)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
print('Patch ({}, {})'.format(i, j))
return pattern
def main():
# Parameters for the pattern.
# For example, using 100x100 pixels image with patches of 3x3 pixels.
height = 1200
width = 1000
patch_size = 12
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1000
pattern = generate_color_pattern(height, width, patch_size, r)
# Scale to 0-255 for display and saving.
pattern_display = (pattern * 255).astype(np.uint8)
# Display the pattern with OpenCV.
cv2.imshow('Pattern', pattern_display)
cv2.waitKey(0)
cv2.destroyAllWindows()
# Save the pattern as an image file.
cv2.imwrite('color_pattern.png', pattern_display)
if __name__ == '__main__':
main()
-141
View File
@@ -1,141 +0,0 @@
#%%
import numpy as np
import cv2
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def generate_color_pattern(h, w, d, r):
"""
Generate a dense color pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
"""
pattern = np.zeros((h, w, 3), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch colors (using mean because patches have uniform color).
neighbors = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
neighbors.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
neighbors.append(np.mean(patch_left, axis=(0, 1)))
if len(neighbors) == 0:
# For the first patch use a random color.
new_color = np.random.rand(3)
else:
new_color = generate_valid_color(neighbors, r)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
print('Patch ({}, {})'.format(i, j))
return pattern
def main():
# Parameters for the pattern.
# For example, using 100x100 pixels image with patches of 3x3 pixels.
# height = 2750
# width = 1950
height = 1200
width = 1000
patch_size = 12
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1000
pattern = generate_color_pattern(height, width, patch_size, r)
# Scale to 0-255 for display and saving.
pattern_display = (pattern * 255).astype(np.uint8)
# Display the pattern with OpenCV.
cv2.imshow('Pattern', pattern_display)
cv2.waitKey(0)
cv2.destroyAllWindows()
# Save the pattern as an image file.
cv2.imwrite('color_pattern.png', pattern_display)
if __name__ == '__main__':
main()
-141
View File
@@ -1,141 +0,0 @@
#%%
import numpy as np
import cv2
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def generate_color_pattern(h, w, d, r):
"""
Generate a dense color pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
"""
pattern = np.zeros((h, w, 3), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch colors (using mean because patches have uniform color).
neighbors = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
neighbors.append(np.mean(patch_above, axis=(0, 1)))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
neighbors.append(np.mean(patch_left, axis=(0, 1)))
if len(neighbors) == 0:
# For the first patch use a random color.
new_color = np.random.rand(3)
else:
new_color = generate_valid_color(neighbors, r)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = new_color
print('Patch ({}, {})'.format(i, j))
return pattern
def main():
# Parameters for the pattern.
# For example, using 100x100 pixels image with patches of 3x3 pixels.
# height = 2750
# width = 1950
height = 1340
width = 1360
patch_size = 12
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1000
pattern = generate_color_pattern(height, width, patch_size, r)
# Scale to 0-255 for display and saving.
pattern_display = (pattern * 255).astype(np.uint8)
# Display the pattern with OpenCV.
cv2.imshow('Pattern', pattern_display)
cv2.waitKey(0)
cv2.destroyAllWindows()
# Save the pattern as an image file.
cv2.imwrite('color_pattern.png', pattern_display)
if __name__ == '__main__':
main()
@@ -1,233 +0,0 @@
#%%
import numpy as np
import cv2
def merge_intervals(intervals):
# Merge overlapping intervals.
if not intervals:
return []
intervals.sort(key=lambda x: x[0])
merged = [intervals[0]]
for current in intervals[1:]:
last = merged[-1]
if current[0] <= last[1]:
merged[-1] = (last[0], max(last[1], current[1]))
else:
merged.append(current)
return merged
def get_allowed_intervals(neighbor_values, delta):
# Compute forbidden intervals for this channel from neighbor values.
forbidden = []
for c in neighbor_values:
low = max(0, c - delta)
high = min(1, c + delta)
forbidden.append((low, high))
# Merge overlapping forbidden intervals.
merged = merge_intervals(forbidden)
allowed = []
# Allowed region before the first forbidden interval.
if merged:
if merged[0][0] > 0:
allowed.append((0, merged[0][0]))
# Allowed regions between forbidden intervals.
for i in range(1, len(merged)):
allowed.append((merged[i-1][1], merged[i][0]))
# Allowed region after the last forbidden interval.
if merged[-1][1] < 1:
allowed.append((merged[-1][1], 1))
else:
allowed.append((0, 1))
return allowed
def sample_from_allowed(allowed):
# Choose a random value from the allowed intervals with probability proportional to interval length.
lengths = [end - start for start, end in allowed]
total_length = sum(lengths)
if total_length <= 0:
return None
point = np.random.rand() * total_length
cumulative = 0
for (start, end), length in zip(allowed, lengths):
cumulative += length
if point <= cumulative:
return np.random.uniform(start, end)
return allowed[-1][0]
def sample_channel(neighbor_values, delta):
# Compute allowed intervals and sample a value.
allowed = get_allowed_intervals(neighbor_values, delta)
if not allowed:
return None
return sample_from_allowed(allowed)
def clip_allowed_to_min(allowed, min_value):
clipped = []
for start, end in allowed:
start = max(start, min_value)
if end > start:
clipped.append((start, end))
return clipped
def generate_valid_gray(neighbor_grays, r, black_threshold):
delta = np.sqrt(r)
min_gray = min_gray_normalized(black_threshold)
allowed = clip_allowed_to_min(get_allowed_intervals(neighbor_grays, delta), min_gray)
if not allowed:
return np.random.uniform(min_gray, 1.0)
sampled = sample_from_allowed(allowed)
if sampled is None:
return np.random.uniform(min_gray, 1.0)
return enforce_min_gray(sampled, black_threshold)
def rgb_to_gray(color):
# ITU-R BT.601 luma, color channels in [0, 1].
r, g, b = color
return 0.299 * r + 0.587 * g + 0.114 * b
def min_gray_normalized(black_threshold):
# Strictly above black_threshold in uint8 space.
return (float(black_threshold) + 1.0) / 255.0
def gray_to_rgb(gray):
gray = float(gray)
return np.array([gray, gray, gray], dtype=np.float32)
def enforce_min_gray(gray, black_threshold):
return max(float(gray), min_gray_normalized(black_threshold))
def generate_valid_color(neighbors, r):
"""
Generate a random color that satisfies, for each neighbor,
|R - R_neighbor| >= sqrt(r), and similarly for G and B.
neighbors: list of neighbor colors (each is an array of shape (3,))
r: threshold for squared difference.
"""
delta = np.sqrt(r)
color = np.zeros(3)
for channel in range(3):
neighbor_channel = [n[channel] for n in neighbors]
sampled = sample_channel(neighbor_channel, delta)
# If the allowed intervals are empty, use a fallback.
if sampled is None:
sampled = np.random.rand()
color[channel] = sampled
return color
def fill_pattern_remainder(pattern, h, w, d, black_threshold):
"""Fill pixels not covered by full dxd patches (image size may not be divisible by d)."""
min_gray = min_gray_normalized(black_threshold)
filled_h = (h // d) * d
filled_w = (w // d) * d
if filled_h < h and filled_w > 0:
remainder_h = h - filled_h
source = pattern[filled_h - d:filled_h, :filled_w]
pattern[filled_h:h, :filled_w] = np.tile(
source, ((remainder_h + d - 1) // d, 1)
)[:remainder_h, :]
if filled_w < w and filled_h > 0:
remainder_w = w - filled_w
source = pattern[:filled_h, filled_w - d:filled_w]
pattern[:filled_h, filled_w:w] = np.tile(source, (1, (remainder_w + d - 1) // d))[:, :remainder_w]
if filled_h < h and filled_w < w:
remainder_h = h - filled_h
remainder_w = w - filled_w
source = pattern[filled_h - d:filled_h, filled_w - d:filled_w]
pattern[filled_h:h, filled_w:w] = np.tile(
source,
((remainder_h + d - 1) // d, (remainder_w + d - 1) // d),
)[:remainder_h, :remainder_w]
pattern[:] = np.maximum(pattern, min_gray)
return pattern
def generate_color_pattern(h, w, d, r, black_threshold):
"""
Generate a dense grayscale pattern.
h, w: height and width of the output image in pixels.
d: patch size in pixels (each patch is d x d).
r: randomness regulation factor. For each neighbor the condition is:
(|channel_new - channel_neighbor|)^2 >= r.
black_threshold: minimum grayscale value in [0, 255]; patches cannot be black.
"""
pattern = np.zeros((h, w), dtype=np.float32)
num_patches_x = h // d
num_patches_y = w // d
min_gray = min_gray_normalized(black_threshold)
for i in range(num_patches_x):
for j in range(num_patches_y):
# Collect neighbor patch grays as RGB triples for the existing sampler.
neighbors = []
neighbor_grays = []
if i > 0:
patch_above = pattern[(i-1)*d:i*d, j*d:(j+1)*d]
gray_above = float(np.mean(patch_above))
neighbor_grays.append(gray_above)
neighbors.append(gray_to_rgb(gray_above))
if j > 0:
patch_left = pattern[i*d:(i+1)*d, (j-1)*d:j*d]
gray_left = float(np.mean(patch_left))
neighbor_grays.append(gray_left)
neighbors.append(gray_to_rgb(gray_left))
if len(neighbors) == 0:
# First patch: sample directly above the black threshold.
gray = np.random.uniform(min_gray, 1.0)
else:
new_color = generate_valid_color(neighbors, r)
gray = rgb_to_gray(new_color)
if gray < min_gray:
gray = generate_valid_gray(neighbor_grays, r, black_threshold)
else:
gray = enforce_min_gray(gray, black_threshold)
pattern[i*d:(i+1)*d, j*d:(j+1)*d] = gray
print('Patch ({}, {}) gray={:.3f}'.format(i, j, gray))
return fill_pattern_remainder(pattern, h, w, d, black_threshold)
def main():
# Parameters for the pattern.
# For example, using 100x100 pixels image with patches of 3x3 pixels.
# height = 2750
# width = 1950
height = 1340
width = 1360
patch_size = 12
black_threshold = 125
# Set r to a value that gives sufficient difference.
# Note: r is the threshold for squared difference per channel, and the condition is
# |c - c_neighbor| >= sqrt(r). Very high r may make it impossible to satisfy.
r = 1000
pattern = generate_color_pattern(height, width, patch_size, r, black_threshold)
# Scale to 0-255 for display and saving.
pattern_display = (pattern * 255).astype(np.uint8)
if pattern_display.min() <= black_threshold:
raise RuntimeError(
f"Generated pattern contains values <= black_threshold ({black_threshold}); "
f"min={int(pattern_display.min())}"
)
# Display the pattern with OpenCV.
cv2.imshow('Pattern', pattern_display)
cv2.waitKey(0)
cv2.destroyAllWindows()
# Save the pattern as an image file.
cv2.imwrite('gray_pattern.png', pattern_display)
if __name__ == '__main__':
main()
Binary file not shown.
Binary file not shown.
Binary file not shown.